The Tibet AS$\gamma$ Collaboration

# Gamma-ray Observation of the Cygnus Region in the 100 TeV Energy Region

M. Amenomori  Department of Physics, Hirosaki University, Hirosaki 036-8561,
Japan  Y. W. Bao  School of Astronomy and Space Science, Nanjing University,
Nanjing 210093, China  X. J. Bi  Key Laboratory of Particle Astrophysics,
Institute of High Energy Physics, Chinese Academy of Sciences, Beijing 100049,
China  D. Chen<EMAIL_ADDRESS>National Astronomical Observatories,
Chinese Academy of Sciences, Beijing 100012, China  T. L. Chen  Department of
Mathematics and Physics, Tibet University, Lhasa 850000, China  W. Y. Chen
Key Laboratory of Particle Astrophysics, Institute of High Energy Physics,
Chinese Academy of Sciences, Beijing 100049, China  Xu Chen  Key Laboratory of
Particle Astrophysics, Institute of High Energy Physics, Chinese Academy of
Sciences, Beijing 100049, China  Y. Chen  School of Astronomy and Space
Science, Nanjing University, Nanjing 210093, China  Cirennima  Department of
Mathematics and Physics, Tibet University, Lhasa 850000, China  S. W. Cui
Department of Physics, Hebei Normal University, Shijiazhuang 050016, China
Danzengluobu  Department of Mathematics and Physics, Tibet University, Lhasa
850000, China  L. K. Ding  Key Laboratory of Particle Astrophysics, Institute
of High Energy Physics, Chinese Academy of Sciences, Beijing 100049, China  J.
H. Fang  Key Laboratory of Particle Astrophysics, Institute of High Energy
Physics, Chinese Academy of Sciences, Beijing 100049, China  University of
Chinese Academy of Sciences, Beijing 100049, China  K. Fang  Key Laboratory of
Particle Astrophysics, Institute of High Energy Physics, Chinese Academy of
Sciences, Beijing 100049, China  C. F. Feng  Institute of Frontier and
Interdisciplinary Science and Key Laboratory of Particle Physics and Particle
Irradiation (MOE), Shandong University, Qingdao 266237, China  Zhaoyang Feng
Key Laboratory of Particle Astrophysics, Institute of High Energy Physics,
Chinese Academy of Sciences, Beijing 100049, China  Z. Y. Feng  Institute of
Modern Physics, SouthWest Jiaotong University, Chengdu 610031, China  Qi Gao
Department of Mathematics and Physics, Tibet University, Lhasa 850000, China
A. Gomi  Faculty of Engineering, Yokohama National University, Yokohama
240-8501, Japan  Q. B. Gou  Key Laboratory of Particle Astrophysics, Institute
of High Energy Physics, Chinese Academy of Sciences, Beijing 100049, China  Y.
Q. Guo  Key Laboratory of Particle Astrophysics, Institute of High Energy
Physics, Chinese Academy of Sciences, Beijing 100049, China  Y. Y. Guo  Key
Laboratory of Particle Astrophysics, Institute of High Energy Physics, Chinese
Academy of Sciences, Beijing 100049, China  H. H. He  Key Laboratory of
Particle Astrophysics, Institute of High Energy Physics, Chinese Academy of
Sciences, Beijing 100049, China  Z. T. He  Department of Physics, Hebei Normal
University, Shijiazhuang 050016, China  K. Hibino  Faculty of Engineering,
Kanagawa University, Yokohama 221-8686, Japan  N. Hotta  Faculty of Education,
Utsunomiya University, Utsunomiya 321-8505, Japan  Haibing Hu  Department of
Mathematics and Physics, Tibet University, Lhasa 850000, China  H. B. Hu  Key
Laboratory of Particle Astrophysics, Institute of High Energy Physics, Chinese
Academy of Sciences, Beijing 100049, China  J. Huang<EMAIL_ADDRESS>Key Laboratory of Particle Astrophysics, Institute of High Energy Physics,
Chinese Academy of Sciences, Beijing 100049, China  H. Y. Jia  Institute of
Modern Physics, SouthWest Jiaotong University, Chengdu 610031, China  L. Jiang
Key Laboratory of Particle Astrophysics, Institute of High Energy Physics,
Chinese Academy of Sciences, Beijing 100049, China  P. Jiang  National
Astronomical Observatories, Chinese Academy of Sciences, Beijing 100012, China
H. B. Jin  National Astronomical Observatories, Chinese Academy of Sciences,
Beijing 100012, China  K. Kasahara  Faculty of Systems Engineering, Shibaura
Institute of Technology, Omiya 330-8570, Japan  Y. Katayose  Faculty of
Engineering, Yokohama National University, Yokohama 240-8501, Japan  C. Kato
Department of Physics, Shinshu University, Matsumoto 390-8621, Japan  S. Kato
Institute for Cosmic Ray Research, University of Tokyo, Kashiwa 277-8582,
Japan  K. Kawata  Institute for Cosmic Ray Research, University of Tokyo,
Kashiwa 277-8582, Japan  M. Kozai  Institute of Space and Astronautical
Science, Japan Aerospace Exploration Agency (ISAS/JAXA), Sagamihara 252-5210,
Japan  D. Kurashige  Faculty of Engineering, Yokohama National University,
Yokohama 240-8501, Japan  Labaciren  Department of Mathematics and Physics,
Tibet University, Lhasa 850000, China  G. M. Le  National Center for Space
Weather, China Meteorological Administration, Beijing 100081, China  A. F. Li
School of Information Science and Engineering, Shandong Agriculture
University, Taian 271018, China  Institute of Frontier and Interdisciplinary
Science and Key Laboratory of Particle Physics and Particle Irradiation (MOE),
Shandong University, Qingdao 266237, China  Key Laboratory of Particle
Astrophysics, Institute of High Energy Physics, Chinese Academy of Sciences,
Beijing 100049, China  H. J. Li  Department of Mathematics and Physics, Tibet
University, Lhasa 850000, China  W. J. Li  Key Laboratory of Particle
Astrophysics, Institute of High Energy Physics, Chinese Academy of Sciences,
Beijing 100049, China  Institute of Modern Physics, SouthWest Jiaotong
University, Chengdu 610031, China  Y. Li  National Astronomical Observatories,
Chinese Academy of Sciences, Beijing 100012, China  Y. H. Lin  Key Laboratory
of Particle Astrophysics, Institute of High Energy Physics, Chinese Academy of
Sciences, Beijing 100049, China  University of Chinese Academy of Sciences,
Beijing 100049, China  B. Liu  Department of Astronomy, School of Physical
Sciences, University of Science and Technology of China, Hefei, Anhui 230026,
China  C. Liu  Key Laboratory of Particle Astrophysics, Institute of High
Energy Physics, Chinese Academy of Sciences, Beijing 100049, China  J. S. Liu
Key Laboratory of Particle Astrophysics, Institute of High Energy Physics,
Chinese Academy of Sciences, Beijing 100049, China  L. Y. Liu  National
Astronomical Observatories, Chinese Academy of Sciences, Beijing 100012, China
M. Y. Liu  Department of Mathematics and Physics, Tibet University, Lhasa
850000, China  W. Liu  Key Laboratory of Particle Astrophysics, Institute of
High Energy Physics, Chinese Academy of Sciences, Beijing 100049, China  X. L.
Liu  National Astronomical Observatories, Chinese Academy of Sciences, Beijing
100012, China  Y.-Q. Lou  Department of Physics and Tsinghua Centre for
Astrophysics (THCA), Tsinghua University, Beijing 100084, China  Tsinghua
University-National Astronomical Observatories of China (NAOC) Joint Research
Center for Astrophysics, Tsinghua University, Beijing 100084, China
Department of Astronomy, Tsinghua University, Beijing 100084, China  H. Lu
Key Laboratory of Particle Astrophysics, Institute of High Energy Physics,
Chinese Academy of Sciences, Beijing 100049, China  X. R. Meng  Department of
Mathematics and Physics, Tibet University, Lhasa 850000, China  K. Munakata
Department of Physics, Shinshu University, Matsumoto 390-8621, Japan  H.
Nakada  Faculty of Engineering, Yokohama National University, Yokohama
240-8501, Japan  Y. Nakamura  Key Laboratory of Particle Astrophysics,
Institute of High Energy Physics, Chinese Academy of Sciences, Beijing 100049,
China  Institute for Cosmic Ray Research, University of Tokyo, Kashiwa
277-8582, Japan  Y. Nakazawa  College of Industrial Technology, Nihon
University, Narashino 275-8575, Japan  H. Nanjo  Department of Physics,
Hirosaki University, Hirosaki 036-8561, Japan  C. C. Ning  Department of
Mathematics and Physics, Tibet University, Lhasa 850000, China  M. Nishizawa
National Institute of Informatics, Tokyo 101-8430, Japan  M. Ohnishi
Institute for Cosmic Ray Research, University of Tokyo, Kashiwa 277-8582,
Japan  T. Ohura  Faculty of Engineering, Yokohama National University,
Yokohama 240-8501, Japan  S. Okukawa  Faculty of Engineering, Yokohama
National University, Yokohama 240-8501, Japan  S. Ozawa  National Institute of
Information and Communications Technology, Tokyo 184-8795, Japan  L. Qian
National Astronomical Observatories, Chinese Academy of Sciences, Beijing
100012, China  X. Qian  National Astronomical Observatories, Chinese Academy
of Sciences, Beijing 100012, China  X. L. Qian  Department of Mechanical and
Electrical Engineering, Shangdong Management University, Jinan 250357, China
X. B. Qu  College of Science, China University of Petroleum, Qingdao 266555,
China  T. Saito  Tokyo Metropolitan College of Industrial Technology, Tokyo
116-8523, Japan  M. Sakata  Department of Physics, Konan University, Kobe
658-8501, Japan  T. Sako  Institute for Cosmic Ray Research, University of
Tokyo, Kashiwa 277-8582, Japan  T. K. Sako<EMAIL_ADDRESS>Institute for Cosmic Ray Research, University of Tokyo, Kashiwa 277-8582,
Japan  J. Shao  Key Laboratory of Particle Astrophysics, Institute of High
Energy Physics, Chinese Academy of Sciences, Beijing 100049, China  Institute
of Frontier and Interdisciplinary Science and Key Laboratory of Particle
Physics and Particle Irradiation (MOE), Shandong University, Qingdao 266237,
China  M. Shibata  Faculty of Engineering, Yokohama National University,
Yokohama 240-8501, Japan  A. Shiomi  College of Industrial Technology, Nihon
University, Narashino 275-8575, Japan  H. Sugimoto  Shonan Institute of
Technology, Fujisawa 251-8511, Japan  W. Takano  Faculty of Engineering,
Kanagawa University, Yokohama 221-8686, Japan  M. Takita
<EMAIL_ADDRESS>Institute for Cosmic Ray Research, University of
Tokyo, Kashiwa 277-8582, Japan  Y. H. Tan  Key Laboratory of Particle
Astrophysics, Institute of High Energy Physics, Chinese Academy of Sciences,
Beijing 100049, China  N. Tateyama  Faculty of Engineering, Kanagawa
University, Yokohama 221-8686, Japan  S. Torii  Research Institute for Science
and Engineering, Waseda University, Tokyo 162-0044, Japan  H. Tsuchiya  Japan
Atomic Energy Agency, Tokai-mura 319-1195, Japan  S. Udo  Faculty of
Engineering, Kanagawa University, Yokohama 221-8686, Japan  H. Wang  Key
Laboratory of Particle Astrophysics, Institute of High Energy Physics, Chinese
Academy of Sciences, Beijing 100049, China  Y. P. Wang  Department of
Mathematics and Physics, Tibet University, Lhasa 850000, China  Wangdui
Department of Mathematics and Physics, Tibet University, Lhasa 850000, China
H. R. Wu  Key Laboratory of Particle Astrophysics, Institute of High Energy
Physics, Chinese Academy of Sciences, Beijing 100049, China  Q. Wu  Department
of Mathematics and Physics, Tibet University, Lhasa 850000, China  J. L. Xu
National Astronomical Observatories, Chinese Academy of Sciences, Beijing
100012, China  L. Xue  Institute of Frontier and Interdisciplinary Science and
Key Laboratory of Particle Physics and Particle Irradiation (MOE), Shandong
University, Qingdao 266237, China  Y. Yamamoto  deceased  Department of
Physics, Konan University, Kobe 658-8501, Japan  Z. Yang  Key Laboratory of
Particle Astrophysics, Institute of High Energy Physics, Chinese Academy of
Sciences, Beijing 100049, China  Y. Q. Yao  National Astronomical
Observatories, Chinese Academy of Sciences, Beijing 100012, China  J. Yin
National Astronomical Observatories, Chinese Academy of Sciences, Beijing
100012, China  Y. Yokoe  Institute for Cosmic Ray Research, University of
Tokyo, Kashiwa 277-8582, Japan  N. P. Yu  National Astronomical Observatories,
Chinese Academy of Sciences, Beijing 100012, China  A. F. Yuan  Department of
Mathematics and Physics, Tibet University, Lhasa 850000, China  L. M. Zhai
National Astronomical Observatories, Chinese Academy of Sciences, Beijing
100012, China  C. P. Zhang  National Astronomical Observatories, Chinese
Academy of Sciences, Beijing 100012, China  H. M. Zhang  Key Laboratory of
Particle Astrophysics, Institute of High Energy Physics, Chinese Academy of
Sciences, Beijing 100049, China  J. L. Zhang  Key Laboratory of Particle
Astrophysics, Institute of High Energy Physics, Chinese Academy of Sciences,
Beijing 100049, China  X. Zhang  School of Astronomy and Space Science,
Nanjing University, Nanjing 210093, China  X. Y. Zhang  Institute of Frontier
and Interdisciplinary Science and Key Laboratory of Particle Physics and
Particle Irradiation (MOE), Shandong University, Qingdao 266237, China  Y.
Zhang  Key Laboratory of Particle Astrophysics, Institute of High Energy
Physics, Chinese Academy of Sciences, Beijing 100049, China  Yi Zhang  Key
Laboratory of Dark Matter and Space Astronomy, Purple Mountain Observatory,
Chinese Academy of Sciences, Nanjing 210034, China  Ying Zhang  Key Laboratory
of Particle Astrophysics, Institute of High Energy Physics, Chinese Academy of
Sciences, Beijing 100049, China  S. P. Zhao  Key Laboratory of Particle
Astrophysics, Institute of High Energy Physics, Chinese Academy of Sciences,
Beijing 100049, China  Zhaxisangzhu  Department of Mathematics and Physics,
Tibet University, Lhasa 850000, China  X. X. Zhou  Institute of Modern
Physics, SouthWest Jiaotong University, Chengdu 610031, China

###### Abstract

We report observations of gamma-ray emissions with energies in the 100 TeV
energy region from the Cygnus region in our Galaxy. Two sources are
significantly detected in the directions of the Cygnus OB1 and OB2
associations. Based on their positional coincidences, we associate one with a
pulsar PSR J2032$+$4127 and the other mainly with a pulsar wind nebula PWN
G75.2+0.1 with the pulsar moving away from its original birthplace situated
around the centroid of the observed gamma-ray emission. This work would
stimulate further studies of particle acceleration mechanisms at these gamma-
ray sources.

##  I Introduction

The Cygnus region is a large active star-forming complex in our Galaxy,
hosting numerous pulsar wind nebulae (PWNe) and supernova remnants (SNRs) as
well as Wolf-Rayet stars, OB associations, open clusters etc. This region is a
natural laboratory to study cosmic-ray acceleration and transport, and has
been observed in various wavebands. Especially interesting in light of TeV
gamma-ray observations are the Cygnus OB1 and OB2 associations.

In the direction of Cygnus OB2, the HEGRA Imaging Air Cherenkov Telescope
(IACT) discovered a gamma-ray source TeV J2032$+$4130 [1, 2], and follow-up
observations have been made by other IACTs including Whipple [3], MAGIC [4]
and VERITAS [5, 6]. They reported gamma-ray fluxes of roughly 3% that of the
Crab nebula above 1 TeV and source extensions of radii 0.1∘–0.2∘. Meanwhile,
air-shower arrays generally reported fluxes much higher than those of IACTs at
multi-TeV energies [7, 8, 9, 10, 11]. For example, ARGO [9] employed a larger
variable-size window with a radius of approximately 1.5∘ at 1 TeV and reported
flux data points 5 –10 times those of IACTs at multi-TeV energies, and the
flux reported at 7 TeV by HAWC with a window of radius 0.7∘ was approximately
twice that based on their point-source analysis [11], which appears to be
consistent with the IACT measurements. A recent HAWC paper [12] reported that
there is a largely extended ‘Cocoon’ region (HAWC J2030$+$409) with a Gaussian
width of $2.13^{\circ}$, counterpart of the GeV Cygnus Cocoon [13], underneath
the gamma-ray emission of HAWC J2031$+$415 with a width of $0.27^{\circ}$,
counterpart of TeV J2032$+$4130\. A plausible explanation of the discrepancy
in the measured fluxes between IACTs and air-shower arrays is that IACTs
extract the emission of HAWC J2031$+$415 by counting the extended emission of
the Cocoon region in estimating the background, while air-shower arrays tend
to integrate part of the gamma-ray emission of the Cocoon region inside their
relatively large window sizes as well as the gamma-ray emission of HAWC
J2031$+$415\. Additionally, HAWC detected eHWC J2030$+$412 above 56 TeV with
an extension of 0.18∘ [14], although the association of eHWC J2030$+$412 with
TeV J2032$+$4130/PSR J2032$+$4127 is a matter of debate, given that eHWC
J2030$+$412 is $0.36^{\circ}$ away from TeV J2032$+$4130 and $0.33^{\circ}$
away from PSR J2032$+$4127.

Located 0.18∘ southeast of the centroid of TeV J2032$+$4130 is the radio
pulsar PSR J2032$+$4127, with the distance from the Earth of 1.4–1.7 kpc [15,
16, 17], a characteristic age of 180 kyr, and a spin-down luminosity of
$1.7\times 10^{35}$ erg/s [18]. PSR J2032$+$4127 forms a binary system with a
massive B0Ve star MT91 213, and a flux enhancement of this binary system in
the 2017 autumn periastron period was observed at the 0.2 –2 TeV energy region
by VERITAS and MAGIC and at the 0.3 –10 keV energy region by Swift-X-Ray
Telescope and NuSTAR [19, 20, 21].

In the direction of Cygnus OB1, a gamma-ray source MGRO J2019$+$37 was
discovered by Milagro at the median energy of 12 TeV with a source extension
of $0.3^{\circ}$ [22], and this source was also recently observed by HAWC [14,
23]. VERITAS observed the same region above 0.6 TeV and separated the gamma-
ray emissions into two sources: VER J2019$+$368 and VER J2016$+$371 [24].
VERITAS also reported that the morphology of VER J2019$+$368 is asymmetrical
with
$\sigma_{\mathrm{major}}(\sigma_{\mathrm{minor}})~{}\approx~{}0.34^{\circ}~{}(0.13^{\circ})$.
The flux of VER J2019$+$368 measured by VERITAS in 2014 [24] is a few times
higher than that in 2018 [6] due to different collection areas of photon
integration windows.

Located 0.36∘ east of the centroid of VER J2019$+$368 is the radio pulsar PSR
J2021$+$3651, having a characteristic age of 17 kyr and a spin-down luminosity
of $3.4\times 10^{36}$ erg/s [25]. PSR J2021$+$3651 is one of the relatively
few pulsars which have gamma-ray pulsations observed by Fermi-LAT [26]. The
distance $D$ to the pulsar has remained a matter of debate. $D\sim 12$ kpc was
suggested by the dispersion measure [25] and $D\geq 5$ kpc by the pulsar
polarization rotation measure [26]. Comparing the HI column density along the
pulsar line of sight with the hydrogen absorbing column density implied $D\sim
10$ kpc [27] or $3-4$ kpc [28], and $D~{}=~{}1.8^{+1.7}_{-1.4}$ kpc was
obtained from the absorption-distance relation using red-clump stars in the
direction of the pulsar [29]. This value is consistent with 1.58 kpc, which is
the estimated distance to Cygnus OB1 [30, 31]. PSR J2021$+$3651 appears to
form a pulsar wind nebula referred to as PWN G75.2$+$0.1 [27] and the radio
and X-ray morphologies of this nebula feature a bright bow-shaped tail
extending westwards from the pulsar, indicating that the pulsar is moving
eastwards with its birthplace as far west as the apparent end of the tail at
$\approx 0.2^{\circ}$ west of the current pulsar position [24, 32]. In
addition, a number of hard X-ray sources were identified along the western
edge of VER J2019$+$368 by $NuSTAR$ space observations including a young
massive stellar cluster [33].

##  II Experiment and Data Analysis

The Tibet air-shower (AS) array has been observing cosmic rays and gamma rays
above TeV energies since 1990 at Yangbajing (90.522∘E, 30.102∘N; 4300 m above
sea level) in Tibet, China [34]. In this work, we use data obtained by the
Tibet air shower array combined with the muon detector array during 719 live
days from 2014 February to 2017 May. Our data analysis method and data
selection criteria as well as the array configuration are the same as
described in our previous papers [35, 36]. We use the Equi-Zenith-Angle method
[37] to estimate the gamma-ray excess count and the number of background
events. Twenty OFF regions are taken, and the radius of the analysis window is
variable depending on the recorded air-shower size from approximately 0.7∘ at
lower energies around 10 TeV to a lower limit of 0.5∘.

In the following section, we discuss two gamma-ray sources detected
significantly above 10 TeV in the directions of Cygnus OB1 and OB2
respectively.

##  III Results and Discussions

### Cygnus OB2

Figure 1(a) shows a detection significance map around the gamma-ray source
detected by this work with photon energies above 10 TeV in the direction of
Cygnus OB2. The sky is gridded in $0.1^{\circ}\times 0.1^{\circ}$ pixels and
the significance value of each pixel calculated according to [43] is smoothed
by a circular search window of radius $R_{\mathrm{w}}$ centered at the pixel.
Assuming a symmetrical 2D Gaussian distribution for the gamma-ray excess, we
fit the events within the $4^{\circ}\times 4^{\circ}$ region around the source
using the unbinned maximum likelihood method. The centroid of gamma-ray
emissions detected at the pre-trial (post-trial) detection significance of
5.3$\sigma$ (4.7$\sigma$) above 10 TeV is estimated at (R.A., Dec.) =
($308.04^{\circ}\pm 0.08^{\circ},41.46^{\circ}\pm 0.06^{\circ}$). We name this
source TASG J2032$+414$. The location of TASG J2032$+414$ is in good agreement
with that of the pulsar PSR J2032$+$4127 and consistent with that of HAWC
J2031$+$415 [12] at the 1.7$\sigma$ level, while it appears to deviate from
that of TeV J2032$+$4130 reported in [20] at the 2.8$\sigma$ level. We also
find that most of the gamma-ray emission detected above 10 TeV is confined
inside a void where the radio (frequency 1420 MHz) [39] and infrared
(wavelength 24 $\mu$m) [40, 41] emissions are very weak. This morphology was
also seen by VERITAS [6].

Figure 2(a) shows the distribution of the number of observed events above 10
TeV as a function of the square of the opening angle $\phi$ between the
estimated arrival direction and the TASG J2032$+414$ centroid. To estimate a
possible source extension, we perform the $\chi^{2}$ fitting of the data with
the function
$A\exp[-\phi^{2}/2(\sigma_{\mathrm{PSF}}^{2}+\sigma_{\mathrm{EXT}}^{2})]+N_{\mathrm{BG}}$
where $A$ and $\sigma_{\mathrm{EXT}}$ are two fitting parameters and
$\sigma_{\mathrm{PSF}}=0.36^{\circ}$ and $N_{\mathrm{BG}}=224.5$ are the point
spread function (PSF) of our instrument above 10 TeV and the number of
background events estimated from the background cosmic-ray data, we get
$\sigma_{\mathrm{EXT}}=0.00^{\circ}\pm 0.14^{\circ}$, which is consistent with
that obtained from the maximum likelihood fitting described above. The
$\chi^{2}$/ndf of the fitting is 33.8/38. With a large error of 0.14∘, the
$\sigma_{\mathrm{EXT}}$ value above 10 TeV does not indicate whether TASG
J2032$+414$ is extended or not even though it is consistent with the previous
measurements at multi-TeV energies by IACTs, ARGO and HAWC within the
2$\sigma$ level [2, 3, 4, 9, 6, 20, 12].

At higher energies above 40 TeV, our source location is consistent with that
of PSR J2032$+$4127, and the estimated value of
$\sigma_{\mathrm{EXT}}=0.16^{\circ}\pm 0.09^{\circ}$ suggests a slight source
extension at the 1.9$\sigma$ level (refer to Supplemental Material Figure S1
and S2). There seems to be a tension between the centroid of our source above
40 TeV and that of eHWC J2030$+$412 above 56 TeV [14] at a statistical
significance of 3.5$\sigma$, which might result from the complex morphology of
the Cygnus Cocoon region as these source locations are obtained under the
assumption of a Gaussian spatial distribution for gamma-ray signals.

Figure 3(a) shows the differential energy spectrum of TASG J2032$+414$ (red
filled squares and downward arrows). Although there is a discrepancy in flux
at multi-TeV energies as explained in the Introduction, our flux data points
above 10 TeV are consistent with previous measurements of IACTs when the
spill-over of gamma-ray signals outside their integration radius is taken into
account. Our spectrum from 10 TeV to 120 TeV can be expressed by a simple
power-law as $dF/dE=N_{0}(E/40~{}\mathrm{\mbox{TeV}})^{-\Gamma}$ where
$N_{0}=(4.13\pm 0.83)\times 10^{-16}$ TeV-1cm-2s-1 is the differential gamma-
ray flux at 40 TeV and $\Gamma=3.12\pm 0.21$ is the spectral index
($\chi^{2}$/ndf = 1.6/4).

Given that the TASG J2032$+414$ centroid is in good agreement with the
location of the radio pulsar PSR J2032$+$4127, the observed gamma rays would
be produced by relativistic electrons injected by the pulsar through the
inverse-Compton scattering (ICS) with synchrotron and ambient photons.
Unfortunately, the time span of our data does not cover the 2017 autumn
periastron period of the binary system PSR J2032$+$4127/MT91 213. TASG
J2032$+414$, therefore, is considered to be associated with the pulsar itself
rather than the whole binary system. VERITAS and MAGIC detected TeV gamma-ray
emissions coincident with TeV J2032$+$4130, which is considered to be the
pulsar wind nebula of PSR J2032$+$4127 with its centroid 0.2∘ away from the
location of PSR J2032$+$4127\. If the pulsar is moving in the southeast
direction along the elongation of the TeV gamma-ray emission as suggested in
[5], it would be natural to consider that gamma rays above 10 TeV produced by
parent electrons with energies of $\gtrsim$ 100 TeV are confined around the
current pulsar location, while parent electrons of TeV gamma rays were
accelerated earlier when the pulsar was at the centroid of TeV J2032$+$4130.

The other two TeV sources reported earlier in this region, VER
J2019$+$407/3HWC J2020$+$403 and 2HWC J2024$+$417∗ are not significantly
detected in this work. We obtain 99% C.L. integral flux upper limits above 10
TeV (100 TeV) of $1.0\times 10^{-13}$ cm-2s-1 ($1.7\times 10^{-15}$ cm-2s-1)
for VER J2019$+$407/3HWC J2020$+$403, and $0.92\times 10^{-13}$ cm-2s-1
($2.0\times 10^{-15}$ cm-2s-1) for 2HWC J2024$+$417∗.

### Cygnus OB1

Figure 1(b) shows a significance map above 10 TeV in the direction of Cygnus
OB1 obtained by this work. The centroid of gamma-ray emissions is estimated at
(R.A., Dec.) = ($304.99^{\circ}\pm 0.11^{\circ},36.84^{\circ}\pm
0.08^{\circ}$) with the pre-trial (post-trial) detection significance of
6.7$\sigma$ (6.2$\sigma$). We name this source TASG J2019$+$368\. The centroid
of TASG J2019$+$368 is consistent with that reported by HAWC [23] within the
$1\sigma$ level, and by VERITAS [6] within the 2$\sigma$ level. Possible
particle acceleration sites are also indicated in the figure: NuSTAR X-ray
sources [33], Wolf-Rayet stars [42] and Fermi-LAT sources [38]. The pulsar PSR
J2021$+$3651, located 0.23∘ east of the TASG J2019$+$368 centroid, has a
nebula extending westwards from the pulsar, PWN G75.2$+$0.1, which is
coincident with the location of TASG J2019$+$368.

Figure 2(b) shows the $\phi^{2}$ distribution of the events observed above 10
TeV. The experimental data can be fitted with a Gaussian function with a
source extension of $\sigma_{\mathrm{EXT}}=0.28^{\circ}~{}\pm~{}0.07^{\circ}$
above 10 TeV, consistent with the extension reported by VERITAS [6] (HAWC
[14]) at the 2.1$\sigma$ (0.3$\sigma$) level. The $\chi^{2}$/ndf of the
fitting is 49.1/38.

At higher energies above 40 TeV, our source location is consistent with that
of VERITAS [6] and HAWC [23, 14], and the estimated value of
$\sigma_{\mathrm{EXT}}=0.22^{\circ}\pm 0.05^{\circ}$ suggests that the
extension could become smaller as photon energy increases (refer to
Supplemental Material Figure S1 and S2). This tendency was also seen earlier
by HAWC [14] ($\sigma_{\mathrm{EXT}}=0.30^{\circ}\pm
0.02^{\circ}(0.20^{\circ}\pm 0.05^{\circ})$ above $\sim$1 TeV (56 TeV)).

Figure 3(b) shows the differential gamma-ray energy spectrum of TASG
J2019$+$368, which is in good agreement with the HAWC spectrum and
consistently connects with the VERITAS spectrum reported in 2014. Our spectrum
can be expressed either as $dF/dE=N_{0}(E/40~{}\mathrm{\mbox{TeV}})^{-\Gamma}$
with $N_{0}=(10.6\pm 1.3)\times 10^{-16}$ TeV-1cm-2s-1 and $\Gamma=2.70\pm
0.13$ ($\chi^{2}$/ndf = 10.4/5), or including an exponential cutoff as
$dF/dE=N_{0}(E/40~{}\mathrm{\mbox{TeV}})^{-\Gamma}\exp(-E/E_{\mathrm{cut}})$
with $N_{0}=(3.6\pm 2.0)\times 10^{-15}$ TeV-1cm-2s-1, $\Gamma=1.6\pm 0.5$ and
$E_{\mathrm{cut}}=44\pm 21$ TeV ($\chi^{2}$/ndf = 3.0/4). The latter fitting
is preferred because the best-fit power-law function of the former fitting
conflicts with the HAWC upper limit at 170 TeV.

The current position of the pulsar PSR J2021$+$3651 is 0.23∘ east of our
gamma-ray emission centroid, and the estimated location of the pulsar’s
birthplace resides in our 1$\sigma$ position error circle, indicating that the
gamma-ray emission observed above 10 TeV by this work would be mainly caused
by relativistic electrons produced by PSR J2021$+$3651 around the time when
the pulsar was very young. If the pulsar birthplace is the gamma-ray emission
centroid obtained by this work, the transverse velocity of the pulsar can be
expressed as $v=420~{}(D/1.8~{}\mathrm{kpc})~{}(\tau/17~{}\mathrm{kyr})^{-1}$
km/s, where $\tau$ is the pulsar age. Among known pulsar velocities [44],
$v=420$ km/s is a plausible value. We consider that the observed gamma-ray
emission is mainly associated with PSR J2021$+$3651, although there might be
some contributions from nearby sources such as the $NuSTAR$ sources.

The other TeV source VER J2016$+$371 reported earlier in this region is not
significantly detected by this work. We obtain a 99% C.L. integral flux upper
limit above 10 TeV (100 TeV) of $0.91\times 10^{-13}$ cm-2s-1 ($1.4\times
10^{-15}$ cm-2s-1) for VER J2016+371.

##  IV Conclusions

Using the Tibet air shower array combined with the underground muon detector
array, we have detected two gamma-ray sources in the direction of the Cygnus
region significantly with photon energies above 10 TeV: TASG J2032$+$414 in
the direction of Cygnus OB2 at (R.A., Dec.) = ($308.04^{\circ}\pm
0.08^{\circ},41.46^{\circ}\pm 0.06^{\circ}$), coincident with PSR
J2032$+$4127, and TASG J2019$+$368 in the direction of Cygnus OB1 at (R.A.,
Dec.) = ($304.99^{\circ}\pm 0.11^{\circ},36.84^{\circ}\pm 0.08^{\circ}$),
coincident with PWN G75.2$+$0.1.

We provide for the first time flux data points from 40 TeV to 120 TeV for the
gamma-ray emission region in the direction of Cygnus OB2, and we confirm the
HAWC spectrum up to 100 TeV and provide flux data points from 120 TeV to 200
TeV for the gamma-ray emission region in the direction of Cygnus OB1. For both
cases, the observed gamma-ray emissions likely result from relativistic
electrons injected by the pulsar through the ICS with synchrotron and ambient
photons. PSR J2021$+$3651 could be moving eastwards with a velocity of $\sim
400$ km/s away from its original birthplace situated around the observed
gamma-ray emission centroid. Future observations and theoretical studies would
shed more light on the physical mechanisms of particle acceleration at these
two high-energy gamma-ray sources.

Figure 1: Significance maps around the two gamma-ray emission sources detected
above 10 TeV in the directions of Cygnus OB2 (a) and OB1 (b), smoothed by
search windows (see the text). The point spread function (PSF) is shown in the
inset figure. The red filled star with a position error circle is the centroid
of TASG J2032$+$414/TASG J2019$+$368 obtained by this work, while the magenta
open cross is the centroid of VER J2031$+$415 [20]/VER J2019$+$368 [6] and the
blue asterisk is that of HAWC J2031$+$415 [12]/3HWC J2019$+$367 [23]. The
green filled diamonds show Fermi-LAT sources [38]. (a) The blue open triangle
indicates the centroid of MAGIC J2031$+$4134 [20]. The green filled diamond
coincident with our gamma-ray emission centroid is the pulsar PSR
J2032$+$4127\. The sky-blue contours indicate 1420 MHz radio emissions
provided by the Canadian Galactic Plane Survey [39], and the pink contours
indicate 24 $\mu$m infrared emissions by the Cygnus-X Spitzer Legacy Survey
[40, 41]. (b) The white open circles are NuSTAR X-ray sources [33], and the
gray filled inverted triangles are Wolf-Rayet stars [42]. The green filled
diamond located at 0.23∘ east of our emission centroid is the pulsar PSR
J2021$+$3651\. The magenta open cross located at (R.A., Dec.) = (303.99∘,
37.21∘) is another VERITAS source VER J2016$+$371 [6], which is not detected
significantly in this work.  Figure 2: Number of events observed with photon
energies above 10 TeV as a function of the square of the opening angle between
the estimated arrival direction and the centroid of TASG J2032$+$414 in (a)
and TASG J2019$+$368 in (b). The red filled circles are the experimental data,
with the best-fit Gaussian function indicated by the solid line. The blue
histogram is the distribution of events expected by the MC simulation assuming
a point-like gamma-ray source.  Figure 3: Differential gamma-ray energy
spectra of (a) TASG J2032$+414$ and (b) TASG J2019$+$368 with 95% C.L. upper
limits measured by this work (red filled squares/arrows). In both panels, the
blue filled circles/arrow (sky-blue open circles) show the gamma-ray spectrum
of reported by VERITAS in 2018 (2014) [6, 5, 24], the gray open
triangles/arrow by HAWC [11, 14], and the dark-green filled circles/arrows by
ARGO [9]. Additionally in panel (a), the gold filled diamonds are reported by
Fermi-LAT [6], the green open squares by MAGIC [4], and the pink
pentagons/arrow by HEGRA [2]. The upper limits of Fermi-LAT, HAWC and VERITAS
are at the 95% confidence level, while those of HEGRA in (a) and ARGO in (b)
are at the 99% and 90% confidence levels, respectively.

###### Acknowledgements.

The collaborative experiment of the Tibet Air Shower Arrays has been conducted
under the auspices of the Ministry of Science and Technology of China and the
Ministry of Foreign Affairs of Japan. This work was supported in part by a
Grant-in-Aid for Scientific Research on Priority Areas from the Ministry of
Education, Culture, Sports, Science and Technology, and by Grants-in-Aid for
Science Research from the Japan Society for the Promotion of Science in Japan.
This work is supported by the National Key R&D Program of China (No.
2016YFE0125500), the Grants from the National Natural Science Foundation of
China (Nos. 11533007, 11673041, 11873065, 11773019, 11773014, 11633007,
11803011, and 11851305), and the Key Laboratory of Particle Astrophysics,
Institute of High Energy Physics, CAS. The research presented in this paper
has used data supplied through the Canadian Galactic Plane Survey. This work
is also supported by the joint research program of the Institute for Cosmic
Ray Research (ICRR), the University of Tokyo.

## References

  * [1] F. Aharonian et al., Astron. Astrophys. 393, L37 (2002).
  * [2] F. Aharonian et al., Astron. Astrophys. 431, 197 (2005).
  * [3] A. Konopelko et al., Astrophys. J. 658, 1062 (2007).
  * [4] J. Albert et al., Astrophys. J. Lett. 675, L25 (2008).
  * [5] E. Aliu et al., Astrophys. J. 783, 16 (2014).
  * [6] A. U. Abeysekara et al., Astrophys. J. 861, 134 (2018).
  * [7] A. A. Abdo et al., Astrophys. J. Lett. 664, L91 (2007).
  * [8] A. A. Abdo et al., Astrophys. J. 753, 159 (2012).
  * [9] B. Bartoli et al., Astrophys. J. Lett. 745, L22 (2012).
  * [10] B. Bartoli et al., Astrophys. J. 790, 152 (2014).
  * [11] A. U. Abeysekara et al., Astrophys. J. 843, 40 (2017).
  * [12] A. U. Abeysekara et al., Nat. Astron. Let. (2021) https://doi.org/10.1038/s41550-021-01318-y.
  * [13] M. Ackermann et al., Science 334, 1103 (2011).
  * [14] A. U. Abeysekara et al., Phys. Rev. Lett. 124, 021102 (2020).
  * [15] P. Massey and A. B. Thompson, Astron. J. 101, 1408 (1991).
  * [16] M. M. Hanson, Astrophys. J. 597, 957 (2003).
  * [17] K. L. J. Rygl, A. Brunthaler, A. Sanna, K. M. Menten, M. J. Reid, H. J. van Langevelde, M. Honma, K. J. E. Torstensson, and K. Fujisawa, Astron. Astrophys. 539, 79 (2012).
  * [18] A. G. Lyne, B. W. Stappers, M. J. Keith, P. S. Ray, M. Kerr, and F. Camilo, Mon. Not. R. Astron. Soc. 451, 581 (2015).
  * [19] W. C. G. Ho, C.-Y. Ng, A. G. Lyne, B. W. Stappers, M. J. Coe, J. P. Halpern, T. J. Johnson, and I. A. Steele, Mon. Not. R. Astron. Soc. 464, 1211 (2017).
  * [20] A. U. Abeysekara et al., Astrophys. J. Lett. 867, L19 (2018).
  * [21] M. Chernyakova, D. Malyshev, P. Blay, B. van Soelen, and S. Tsygankov, Mon. Not. R. Astron. Soc. 495, 365 (2020).
  * [22] A. A. Abdo et al., Astrophys. J. Lett. 658, L33 (2007).
  * [23] A. Albert et al., Astrophys. J. 905, 76 (2020).
  * [24] E. Aliu et al., Astrophys. J. 788, 78 (2014).
  * [25] M. S. E. Roberts, J. W. T. Hessels, S. M. Ransom, V. M. Kaspi, P. C. C. Freire, F. Crawford, and D. R. Lorimer, Astrophys. J. Lett. 577, L19 (2002).
  * [26] A. A. Abdo et al., Astrophys. J. 700, 1059 (2009).
  * [27] J. W. T. Hessels, M. S. E. Roberts, S. M. Ransom, V. M. Kaspi, R. W. Romani, C.-Y. Ng, P. C. C. Freire, and B. M. Gaensler, Astrophys. J. 612, 389 (2004).
  * [28] A. van Etten, R. W. Romani, and C.-Y. Ng, Astrophys. J. 680, 1417 (2008).
  * [29] A. Kirichenko, A. Danilenko, P. Shternin, Y. Shibanov, E. Ryspaeva, D. Zyuzin, M. Durant, O. Kargaltsev, G. Pavlov, and A. Cabrera-Lavers, Astrophys. J. 802, 17 (2015).
  * [30] E. M. Levesque, P. Massey, K. A. G. Olsen, B. Plez, E. Josselin, A. Maeder, and G. Meynet, Astrophys. J. 628, 973 (2005).
  * [31] N. Mauron and E. Josselin, Astron. Astrophys. 526, 156 (2011).
  * [32] T. Mizuno, N. Tanaka, H. Takahashi, J. Katsuta, K. Hayashi, and R. Yamazaki, Astrophys. J. 841, 104 (2017).
  * [33] E. V. Gotthelf et al., Astrophys. J. 826, 25 (2016).
  * [34] M. Amenomori et al., Phys. Rev. Lett. 69, 2468 (1992).
  * [35] M. Amenomori et al., Phys. Rev. Lett. 123, 051101 (2019).
  * [36] M. Amenomori et al., Nat. Astron. Let. (2021) https://doi.org/10.1038/s41550-020-01294-9.
  * [37] M. Amenomori et al., Astrophys. J. 633, 1005 (2005).
  * [38] S. Abdollahi et al., Astrophys. J. Suppl. Ser. 247, 33 (2020).
  * [39] A. R. Taylor et al., Astron. J. 125, 3145 (2003).
  * [40] I. M. Beerer et al., Astrophys. J. 720, 679 (2010).
  * [41] K. E. Kraemer et al., Astron. J. 139, 2319 (2010).
  * [42] K. A. van der Hucht, New Astron. Rev. 45, 135 (2001).
  * [43] T.-P. Li and Y.-Q. Ma, Astrophys. J. 272, 317 (1983).
  * [44] G. Hobbs, D. R. Lorimer, A. G. Lyne, and M. Kramer, Mon. Not. R. Astron. Soc. 360, 974 (2005).