# High energy gamma-ray emission powered by a young protostar: the case of
S255 NIRS 3

Emma de Oña Wilhelmi1, Rubén López-Coto2, Yang Su3
1 Deutsches Elektronen-Synchrotron DESY, Platanenallee 6, 15738 Zeuthen,
Germany
2 Instituto de Astrofísica de Andalucía, CSIC, 18080 Granada, Spain.
3Purple Mountain Observatory and Key Laboratory of Radio Astronomy, Chinese
Academy of Sciences, Nanjing 210034, China
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Evidence of efficient acceleration of cosmic rays in massive young stellar
objects has been recently reported. Among these massive protostars, S255 NIRS
3 for which extreme flaring events associated with radio jets have been
detected, is one of the best objects to test this hypothesis. We search for
gamma-ray emission associated with this object in Fermi-LAT data and inspect
the gas content in different molecular lines using the MWISP survey. A GeV
source dubbed 4FGL J0613.1+1749c lies on top of the MYSO region, where two
filamentary $\sim$10 pc CO structures extend along the same direction of the
sub-parsec radio jets. We investigate the spectrum, morphology, and light
curve of the gamma-ray source and compare it with the theoretical emission
expected from hadronic and leptonic populations accelerated in the radio jets.
We argue that the gamma-ray source could be powered by particles accelerated
in the S255 NIRS 3 jets, radiating via Bremsstrahlung or proton-proton
interaction, and with a synchrotron component shinning in radio from primary
or secondary electrons in the case of a leptonic or hadronic population.

###### keywords:

stars: flare, stars: protostars, stars: massive, stars: individual object:
S255 NIRS 3, gamma-rays: stars

††pubyear: ††pagerange: High energy gamma-ray emission powered by a young
protostar: the case of S255 NIRS 3–High energy gamma-ray emission powered by a
young protostar: the case of S255 NIRS 3

##  1 Introduction

The recent discoveries of high energy gamma-ray radiation from stellar winds,
such as novae or stellar clusters have demonstrated the efficient acceleration
of particles in shocks associated with stellar objects (Cheung et al., 2022;
H. E. S. S. Collaboration et al., 2022a; Ackermann et al., 2014; Acciari et
al., 2022; Aharonian et al., 2019; Abdo et al., 2010). Related to them,
massive young stellar objects (MYSOs) have also been hypothesized as cosmic-
ray accelerators, as they power fast outflows associated to accretion
mechanisms or supersonic ejecta (Araudo et al., 2021; Araudo et al., 2007;
Padovani et al., 2015; Bosch-Ramon et al., 2010). The observation of non-
thermal emission from jets associated with some of these objects indicates the
presence of relativistic electrons (Brooks et al., 2003; Garay et al., 2003;
Carrasco-González et al., 2010; Marti et al., 1993). These outflows appear at
linear distances of less than a parsec from the MYSO, where the environment is
strongly affected by its dense IR photon field and high-density matter from
the hosting cloud. The population of relativistic particles might potentially
radiate in gamma rays via leptonic (inverse Compton and Bremsstrahlung) and/or
hadronic (proton-proton interaction) processes.

S255 NIRS 3 (or G192.60–00.04) is a well-studied $\sim$20 M⊙ high-mass young
stellar object (HMYSO) in the S255IR massive star-forming region, located at a
distance of $\sim$1.8 kpc (Burns et al., 2016). It exhibits a disk-like
rotating structure associated with an accretion disk. A molecular bipolar
outflow or jet-like outflow (with a size of $\sim 1$ arcmin) has also been
detected perpendicular to the disk at position angle 67o (Zinchenko et al.,
2015). Several measurements indicate variability of such jet, pointing to
irregular outburst episodes, that occur every few hundred years (Cesaroni et
al., 2018; Burns et al., 2016). This type of massive star is believed to be
formed through disk-mediated accretion. S255 NIRS 3 gained appreciable
attention for its extreme flaring signatures in 6.7 GHz methanol maser, IR,
and radio continuum (Szymczak et al., 2018; Burns et al., 2016; Stecklum et
al., 2016; Caratti o Garatti et al., 2017; Cesaroni et al., 2018). The
detection of a sudden increase in the FIR stellar bolometric luminosity, from
2.9$\times 10^{4}$L⊙ to 1.6$\times 10^{5}$L⊙ at around 30-100 $\mu$m band, was
invoked to demonstrate the irregularity of the accretion flow, rather than a
smooth, steady one. These accretion events are believed to trigger an increase
in the mass ejection rate, which results in an enhancement of the photon field
in the surroundings of the star and eventually a jet ejection.

Motivated by this discovery, we explore the high energy emission using data
from the LAT telescope on board the Fermi satellite. A source dubbed 4FGL
J0613.1+1749c111The c assigned to the name designs a region coincident with
interstellar clump listed in the Fermi-LAT Fourth Source Catalog (4FGL-DR3,
Abdollahi et al. 2020) lies at RAJ2000=93o.28, DECJ2000=17o.822, coincident
with the position of the MYSO source.

In the following, we investigate the connection between the LAT source with
S255 NIRS 3, and discuss the implications of such a connection, under the
hypothesis of a non-thermal jet powered by the protostar.

##  2 Data Analysis and Results

Figure 1: Spectral energy distribution derived from 4FGL J0613.1+1749c. The
1-$\sigma$ error in the best-fit parameters is shown with a shaded blue area.
The dashed and dotted lines show the best-fit models obtained for hadronic and
leptonic emission mechanisms.

###  2.1 Fermi-LAT analysis

To investigate the nature of the gamma-ray source located at the position of
the MYSO source, we obtained Fermi-LAT data from 4th August 2008 to 25th
January 2023. The Fermi-LAT (P8R3, pass 8 release 3) data were analyzed in the
energy range from 100 MeV to 600 GeV, around the region of interest (ROI),
which was defined by a radius of $20^{\circ}$ around the position of S255 NIRS
3. The analysis of the LAT data described below is performed by means of the
fermipy python package (version 0.18.0), based on the Fermi Science Tools
(version 2.0.8, Wood et al. 2017; Fermi Science Support Development Team
2019). The SOURCE class events with a maximum zenith angle of $90^{\circ}$
were selected to eliminate Earth limb events. To derive the global spectrum,
the data are binned in 8 energy bins per decade and spatial bins of
$0.1^{\circ}$ size. The response of the LAT instrument is evaluated with the
IRFs (version P8R3_SOURCE_V3). The gamma-ray emission in the ROI was described
using all LAT sources listed in the Fermi-LAT 4FGL-DR3 catalog in a radius of
$25^{\circ}$ around the position of the MYSO. The contributions from Galactic
and extra-galactic diffuse gamma-rays are modeled using the Galactic
(gll_iem_v07) and isotropic (iso_P8R3_SOURCE_V3_v1) diffuse emission models,
available from the Fermi Science Support Center (FSSC)
222http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html. The
energy dispersion correction is applied to all sources in our model, except
for the isotropic diffuse emission model. Additional gamma-ray emission
regions which might not be accounted for in the 4FGL catalog were investigated
using the findsource method, applying a criterium of Test Statistics TS$>$25,
resulting in one additional source which was added to the model (located at
(RA, DEC)J2000=(90.88,18.34)o). The inclusion of this source in the model had
no effect on the fit results of the source of interest. The spectral
parameters of the sources within $4^{\circ}$ of S255 NIRS 3, sources with TS
larger than 500, and the Galactic and isotropic diffuse emission components
were left free in the likelihood fit. To compute the lightcurve, the data set
was split into 10 bins using the lightcurve method in fermipy.

Figure 2: Gamma-ray light-curve of the observations ranging 14 years. The
dashed vertical line shows the beginning of the maser flare reported by
Fujisawa et al. (2015) whereas the dotted line marks the beginning of the
radio flare reported by Cesaroni et al. (2018). The shaded area shows the
1$\sigma$ error on the photon flux integrated over the analysis energy range.

The spectral and morphological results from the fit for the GeV source
coincident with S255 NIRS 3 are compatible with the ones in the 4FGL catalog:
4FGL J0613.1+1749c is well-described (with test-statistics of TS=113.6) by a
point-like source with a flux well-described by a power-law spectral function
of the form $\phi(E)=\phi_{0}\left(\frac{E}{E_{o}}\right)^{-\alpha}$, with
$\phi_{0}=(6.21\pm 0.73)\times 10^{-13}$ MeV-1cm-2s-1 and $\alpha=2.54\pm
0.07$ (Eo fixed to 1.231 GeV). The spectral energy distribution is shown in
Fig. 1. To test the extension of the source, we selected only events with
energy above 3 GeV to optimize the LAT angular
resolution333https://www.slac.stanford.edu/exp/glast/groups/canda/lat_Performance.htm.
We obtain an upper limit (95 percent CL) on the extension of 0.25o.

The localization of the point-like source results in a best-fit position at
RAJ2000 =(93.27$\pm$0.03)o and DEC2000 = (17.96$\pm$0.04)o, well in agreement
with the position of the MYSO source. The contours above 12$\sigma$ obtained
from the skymap produced for events with energy higher than 1 GeV (selected as
a compromise between enough statistics and angular resolution) is shown in
Fig. 3 in white.

The light curve was derived by dividing the data set into 10 equal bins in
time. We chose 10 bins as a good compromise between time resolution and
detection significance. Figure 2 shows the integral flux for each time bin on
the upper panel, and the TS on the lower one. An upper limit on the flux (95
percent CL) is displayed if the test statistics (TS) is less than 4. The two
vertical lines mark the start of the 6.7 GHz methanol maser reported by
Fujisawa et al. (2015) and the beginning of the radio flare reported by
Cesaroni et al. (2018) (dashed and dotted lines, respectively). The gamma-ray
source is detected well above 20 TS at the moment of the radio outburst,
whereas before 2016 the flux level seems to be very low. An indication for
flux variability is found at a 2.7${\sigma}$ level ($\chi^{2}$/ndf=22.77/9).
However, at this significance, this variability might well be due to
statistical fluctuations and no firm statements can be made.

###  2.2 CO observations

The molecular complex in this region has been studied in detail in several
lines and scales (Richardson et al., 1985; Chavarría et al., 2008; Bieging et
al., 2009; Zinchenko et al., 2012; Zinchenko et al., 2015; Burns et al., 2016;
Kohno et al., 2022). We investigated the large-scale CO emission on the
surrounding medium of S255 NIRS 3 using deep observations in the region from
the MWISP survey (see details in Su et al. 2019). The 12CO(J=1–0),
13CO(J=1–0), and C18O (J=1–0) lines were simultaneously observed using the
13.7 m telescope. The covered region is towards the northern Galactic plane
with $|b|<5^{\rm o}$. The half-power beam width (HPBW) is $\sim$50 arcsec at
the frequency of $\sim 110-115\,\rm{GHz}$. The rms noise is $\sim 0.5\,$K for
12CO and $\sim$0.3 K for 13CO and C18O, at a velocity resolution of $\sim$0.2
km s-1 with a uniform grid spacing of 30 arcsec. More details about these
observations and a very detailed analysis of the dense cloud around S255 NIRS
3 can be found in Wang et al. (2017). The data cubes were analyzed using the
CARTA image visualization and analysis tool (Comrie et al., 2020).

The velocity profile of the region around S255 NIRS 3 is dominated for a
prominent single peak at v${}_{\rm LST}\sim$7.5 km s-1, which is close to the
H2O maser’s velocity of 6.4 km s-1 for S255 NIRS 3 (see Figures 1-2 in Hirota
et al. 2021). We thus assume that the dense molecular cloud is at a distance
of 1.78 kpc (Burns et al., 2016) based on the physical association between the
filamentary cloud and the masers. The CO intensity is dominated by the stellar
complex S255. A high-density peak can be found at the position of the cluster
at a velocity range between 4 km s-1 and 11 km s-1. In Fig. 3 we show the 13CO
image integrated in this velocity range. Overlaid with green contours we
plotted the integrated intensity of C18O (see also Fig. 8 in Wang et al.
2017). The distribution of significance of the GeV emission is peaked in the
central region of S255 NIRS 3, where the maximum of the CO emission is
observed. Interestingly, the two distinct filamentary structures reported by
Wang et al. 2017; Bieging et al. 2009 are aligned with an angle of $\sim$65o
with the jet structure reported by several authors (Wang et al., 2011;
Zinchenko et al., 2015), but at much larger scale ($\sim$10 pc) (see Fig. 3,
dashed blue line).

Figure 3: Integrated 13CO intensity map of the S255 NIRS 3 region obtained
from the MWISP survey. The green contours are obtained from the C18O intensity
map, tracing the two extense filamentary structures emerging from the core of
the emission. The contours above 12 TS (in steps of 2 TS) obtained from the
LAT Test-Statistic skymap above an energy threshold of 1 GeV are overlaid in
white and the best-fit position in cyan. The direction of the 1 arcmin
molecular bipolar outflow is shown with a blue dashed line. The orange circle
marks the 95 percent CL on the source extension.

The molecular gas mass directly related to the S255 NIRS 3 (the dust lane or
the dense gas filament at VLSR=(4-11) km s-1) is about $10^{3}$M⊙ (in a region
of 17 arcmin length and 4 arcmin width, or 8.8$\times$ 2.1 pc2), while we
estimate a total molecular gas mass of about 104M⊙ in the region where the
gamma-ray source sits (see white contours on Fig. 3). The mean density of the
molecular gas in the filament is about 300-400 H2 cm-3 or even higher. These
rich and dense region provides sufficient target to promote gamma-ray emission
via Bremsstrahlung or proton-proton interaction if relativistic electrons
and/or protons are accelerated somewhere in the region.

##  3 Discussion

The rich multi-wavelength data around the MYSO source S255 NIRS 3 allows us to
investigate the possibility of gamma-ray high energy emission associated with
accretion-ejection phenomena in the massive protostar. This particular MYSO,
with a (steady) luminosity of 5$\times 10^{4}$L⊙, is one of the few objects
that has experienced rare but intense outbursts of radio emission, in which an
increase in the luminosity up to 1.3$\times 10^{5}$L⊙ was measured (Caratti o
Garatti et al., 2017). These outbursts are believed to be recurrent and are
associated with phenomena related to disk accretion and the ejection of a jet.
In 2015 a flare of the Class II methanol maser was reported (Fujisawa et al.,
2015). Follow-up mid-infrared and far-infrared observations in November 2015
showed a sudden increase in the brightness of the central object. Later in
July 2016 (13 months later from the maser outburst), a radio jet burst
triggered by the accretion burst was also detected. This outburst released an
estimated energy of 1.2$\times 10^{46}$erg (calculated from the beginning of
the burst until mid-April 2016, Caratti o Garatti et al. 2017). Using high-
resolution Jansky Very Large Array (JVLA) observations and comparing
observations from 2012 and 2018, Obonyo et al. (2021) reported the presence of
two non-thermal radio lobes associated with the protostar, with energy index
(S${}_{\nu}\propto\nu^{-\alpha}$) of $\alpha\simeq 0.6$, corresponding to an
electron index of p = 2.2. The latest is in good agreement with what is
expected from an uncooled population of relativistic electrons produced by
diffusive shock acceleration (DSA) at a strong non-relativistic shock (e.g.
Protheroe 1998). These lobes are located at a projected distance of 0.02 and
0.04 pc (or $6.2\times 10^{16}$ and $1.2\times 10^{17}$ cm, respectively) from
the compact center. They are also aligned in a direction similar (with an
angle of $\approx$67o) to the orientation of the radio emission of the compact
center, the infrared emission in the field, and the observed CO filaments (see
Fig. 3) pointing to a jet-like emission powered by the central object. Based
on equipartition considerations, Obonyo et al. (2021) estimated a minimum
energy and a minimum magnetic field strength of $4\times 10^{41}$ erg and 1.1
mG, respectively. These radio lobes are probably generated by the strong
terminal shocks off the jets, which have an estimated velocity of $\sim$700 km
s-1 (Obonyo et al., 2021; Cesaroni et al., 2018). This value is within the
range of jet velocity expected from jets in these objects (see Marti et al.
1993, 1995).

Unfortunately, the lightcurve cannot establish an unambiguous link between the
observed flare and the gamma-ray emission reported here. Even if the increase
of the significance and flux towards the time of the flare makes this
possibility attractive, the statistical analysis of the data shows no
significant variability further than 2.7$\sigma$ (without accounting for
trials and systematic errors), preventing such a connection.

However, the presence of a shock in which relativistic electrons are
accelerated and the amplified magnetic field, in a dense region, suggests that
the origin of 4FGL J0613.1+1749c should be related somehow to S255 NIRS 3. The
potential of massive young objects to power gamma-ray sources was investigated
by several authors (see e.g. Bosch-Ramon et al. 2010; Araudo et al. 2007;
Araudo et al. 2021), and applied to several MYSO with similar parameters to
the ones found in S255 NIRS 3. We follow thus the approach described in Bosch-
Ramon et al. (2010) to explore the expected gamma-ray emission in case of
electrons and protons accelerated in jets associated with S255 NIRS 3, and
compare with the observations described above.

Figure 4: Relevant timescales as a function of the particle energy applied to
the massive young stellar object S255 NIRS 3.

S255 NIRS 3 is located at a distance of 1.8 kpc and the total luminosity above
100 MeV amounts L${}_{\gamma}=(4.2\pm 0.5)\times 10^{33}$ erg s-1.
Relativistic particles, hadrons or leptons, can be accelerated in jets, even
for moderated flow speed (H. E. S. S. Collaboration et al., 2022b; Acciari et
al., 2022). For electrons, the relevant mechanisms ruling the observed
spectral energy distribution are synchrotron, inverse Compton, and
relativistic Bremsstrahlung, while for protons, inelastic collision via
proton-proton would originate the observed gamma-ray emission, whereas the
radio emission would correspond to synchrotron emission from secondary
electrons. We used the estimated magnetic field of 1 mG, which is also in
agreement with the typical magnetic field obtained from Zeeman measurement in
dense clouds like the one in S255 NIRS 3 (see Section 2.2). Relativistic
electrons can upscatter the dense FIR photon field in the region to the GeV
regime through inverse Compton. We calculated the IR photon density in the
region of the jet using the one reported during the flare ($1.6\times
10^{5}L_{\sun}$), to account for the maximum-possible contribution. For
radiation mechanisms depending on the gas density (i.e. Bremsstrahlung and
proton-proton collisions), we used a minimum gas density at the location of
the jets of n${}_{\rm H}\simeq 200$ cm-3, which is also in agreement with the
values of the jet density obtained in Bosch-Ramon et al. (2010).

In Fig. 4 we show the cooling times for the mechanisms described above, as
well as the relevant times for the shock evolution in which particles are
believed to be accelerated such as the jet lifetime, diffusive time for
particles in the Bohm limit and acceleration time, for a jet of velocity and
magnetic field of 700 km s-1 and 1 mG. The dominant gamma-ray emission
processes are proton-proton interaction for hadrons and Bremsstrahlung in
electrons (black and orange dashed lines, respectively) with the maximum
energy determined by the acceleration time in the jets (purple solid line).
For inverse Compton to dominate over Bremsstrahlung, the density has to be
lower, requiring a larger emission region of the order of a few parsecs.
However, in such a larger region, the IR photon density is also greatly
suppressed, resulting in a negligible contribution of the inverse Compton
component. This larger region where the large-scale CO filamentary structure
is found is still confined within the angular extension for a point source in
LAT. Alternatively to the emission region located on the jet shock region,
recurrent shocks could have injected particles that diffuse out into a region
of $\sim$10 pc around the protostar. However given the large density in the
sub-parsec region around S255 NIRS 3 only high-energy particles would escape
without losing their energies via Bremsstrahlung/proton-proton radiation,
rendering this scenario less likely.

If the gamma-ray emission is indeed associated with the non-thermal jets
observed in S255 NIRS 3, we can obtain the electron and/or proton population
behind the radio and GeV emission, by fitting the spectral energy distribution
using the Naima library within the framework of the Gammapy software (Zabalza,
2015; Deil et al., 2017). In both cases, the particle spectrum was described
with a power-law function with an exponential cutoff in energy. A single
population of electrons radiating via synchrotron and Bremsstrahlung in the
radio and GeV energy range, respectively, can in principle reproduce the
spectral shape measured. The spectral index of the radio emission in the jets
constrains the local electron power-law index, which is compatible with the
best-fit particle spectral index obtained ($\rm{p}=1.8\pm 0.7$), with an
unconstrained energy cutoff at $\rm{E_{c,e}}\simeq 1.9$ GeV (see dotted line
in Fig. 1). Using the expression in Bosch-Ramon et al. 2010, for primary
electrons:

$B\sim 0.04(L_{35}^{\rm e})^{-2/3}R_{\rm 16}^{-2/3}\nu_{5GHz}^{1/3}d_{\rm
3~{}kpc}^{4/3}~{}F_{\rm\nu~{}mJy}^{2/3}~{}{\rm mG},$ (1)

with $L_{35}^{\rm e}$ the luminosity of the electron population in erg s-1, R
the size of the jets in units of $10^{16}$cm, and $F_{\rm\nu}$ the measured
radio flux at $\nu$=5 GHz ($\sim$0.65 mJy, see Obonyo et al. 2021), results on
a magnetic field strength of less than 0.1 mG, which is below the minimum
magnetic field estimated in the region.

The fit to a proton-proton distribution with particle spectral index fixed to
2 for the uncooled energy spectrum, as expected by DSA in shocks (e.g. Drury
1983), reproduces the GeV data correctly with an energy cutoff of
$\rm{E_{c,p}}=3.1\pm 0.9$ GeV (see dashed lines in Fig. 1). The radio lobes
can in the scenario be attributed to synchrotron emission from the secondary
electrons produced in the GeV hadronic emission. By comparing the radio and
GeV flux (using the expression above, and substituting $L^{e}_{35}$ by its
equivalent for secondary electrons, 0.1(tesc/t${}_{\rm pp})L^{p}_{35}$) a
magnetic field of 1 mG can be estimated, which agrees well with the magnetic
field obtained in Obonyo et al. 2021, and also with the magnetic field
strength considered in theoretical models. The hadronic origin fits well the
data in a one-zone scenario, however, it should be noted that a leptonic
origin cannot be discarded, given that the lower magnetic field required by
the synchrotron emission can be compensated by a moderate increase of the
density on the emission region. In such a scenario though, a departure from
the equipartition hypothesis has to be assumed.

##  4 Conclusions

We investigated the possibility of GeV gamma-ray emission related to the
accretion-ejection phenomena in the massive protostar S255 NIRS 3. We found a
point-like gamma-ray source, dubbed in the 4FGL catalog as 4FGL J0613.1+1749c,
located at a position compatible with the one of the protostar. In general,
associating MYSO objects with gamma-ray sources in terms of positional
coincident is difficult, given the ubiquity of these kinds of objects in
massive clouds, which do correlate with the gamma-ray distribution of sources
in the Galaxy (Pedaletti et al., 2014; Munar-Adrover et al., 2011).
Additionally, several massive OB stars and two non-thermal X-ray sources
associated with background AGNs, are located within the region of the gamma-
ray source (Mucciarelli et al., 2011), adding uncertainty to the
identification of the gamma-ray source with the MYSO. Similar massive stars
are found in high-energy binary systems such e.g. LS 5039 (Motch et al.,
1997), however, no compact object has been found associated with these ones in
the S255 complex. No other types of sources usually associated with Galactic
GeV gamma-ray radiation, such as pulsars or supernova remnants are found
either.

The location of S255 NIRS 3 close to the anticenter, where the stellar
formation regions are not so common, the lack of any other obvious
counterpart, and the point-like nature of the gamma-ray source 4FGL
J0613.1+1749c make the association with S255 NIRS 3 appealing. Additionally,
the observed spectral energy distribution in radio and GeV wavelengths is in
good agreement with the emission expected theoretically, when applied to the
specific conditions of this object. The molecular-rich environment provides
enough targets to explain the gamma-ray emission proton-proton interactions,
with protons accelerated in the observed jet (or in recurrent ones). The
estimated magnetic field in such jets is also in agreement with the one
required to explain the radio emission by synchrotron emission from secondary
electrons.

The CO image also unveiled a large-scale molecular filament, oriented in the
same direction as the radio compact jets, which could be originated from
recurrent ejections powered by the protostar.

A correlation through the lightcurve between the MWL flare and the gamma-ray
emission would leave no doubt about the origin of the gamma-ray source.
However, the hint of variability is too low to claim the flare as powering the
gamma-ray source. Two more of these sources with similar environments, that
is, massive young protostars with jets, have been also proposed as
counterparts of gamma-ray sources, namely HH 211 and HH 80–81 (Yang et al.,
2023; Yan et al., 2022).

Monitoring of this kind of flares should clarify if the observed gamma-ray
radiation can be univocally attributed the MYSO activities, probing the role
of cosmic rays in the evolution of massive young stellar objects.

## Acknowledgments

EdOW acknowledges the support of DESY (Zeuthen), a member of the Helmholtz
Association HGF. R.L.-C. acknowledges the Ramón y Cajal program through grant
RYC-2020-028639-I and the financial support from the grant CEX2021-001131-S
funded by MCIN/AEI/ 10.13039/501100011033. We thank Rolf Bühler for useful
comments on the LAT analysis and Felix Aharonian, Brian Reville and Anabella
Araudo for helpful discussions on the interpretation. We also thank to Alessio
Caratti o Garatti for his insights on the physics of this object. This
research made use of the data from the Milky Way Imaging Scroll Painting
(MWISP) project, which is a multiline survey in 12CO/13CO/C18O along the
northern Galactic plane with the PMO 13.7m telescope. We are grateful to all
the members of the MWISP working group, particularly the staff members at the
PMO 13.7m telescope, for their long-term support. MWISP was sponsored by the
National Key R&D Program of China with grant 2017YFA0402700 and the CAS Key
Research Program of Frontier Sciences with grant QYZDJ-SSW-SLH047.

## Data Availability

The Fermi data underlying this article are available at
https://fermi.gsfc.nasa.gov/ssc/data/access/lat/. The CO data was accessed
from the MWISP survey and can be shared on reasonable request to the
corresponding author. This publication uses data products from the Wide-field
Infrared Survey Explorer, a joint project of the University of California, Los
Angeles, and the Jet Propulsion Laboratory/California Institute of Technology,
funded by the National Aeronautics and Space Administration.

## References

  * Abdo et al. (2010) Abdo A. A., et al., 2010, Science, 329, 817
  * Abdollahi et al. (2020) Abdollahi S., et al., 2020, ApJS, 247, 33
  * Acciari et al. (2022) Acciari V. A., et al., 2022, Nature Astronomy, 6, 689
  * Ackermann et al. (2014) Ackermann M., et al., 2014, Science, 345, 554
  * Aharonian et al. (2019) Aharonian F., Yang R., de Oña Wilhelmi E., 2019, Nature Astronomy, 3, 561
  * Araudo et al. (2007) Araudo A. T., Romero G. E., Bosch-Ramon V., Paredes J. M., 2007, A&A, 476, 1289
  * Araudo et al. (2021) Araudo A. T., Padovani M., Marcowith A., 2021, MNRAS, 504, 2405
  * Bieging et al. (2009) Bieging J. H., Peters W. L., Vila Vilaro B., Schlottman K., Kulesa C., 2009, AJ, 138, 975
  * Bosch-Ramon et al. (2010) Bosch-Ramon V., Romero G. E., Araudo A. T., Paredes J. M., 2010, A&A, 511, A8
  * Brooks et al. (2003) Brooks K. J., Garay G., Mardones D., Bronfman L., 2003, ApJ, 594, L131
  * Burns et al. (2016) Burns R. A., Handa T., Nagayama T., Sunada K., Omodaka T., 2016, MNRAS, 460, 283
  * Caratti o Garatti et al. (2017) Caratti o Garatti A., et al., 2017, Nature Physics, 13, 276
  * Carrasco-González et al. (2010) Carrasco-González C., Rodríguez L. F., Anglada G., Martí J., Torrelles J. M., Osorio M., 2010, Science, 330, 1209
  * Cesaroni et al. (2018) Cesaroni R., et al., 2018, A&A, 612, A103
  * Chavarría et al. (2008) Chavarría L. A., Allen L. E., Hora J. L., Brunt C. M., Fazio G. G., 2008, ApJ, 682, 445
  * Cheung et al. (2022) Cheung C. C., et al., 2022, ApJ, 935, 44
  * Comrie et al. (2020) Comrie A., et al., 2020, CARTA: The Cube Analysis and Rendering Tool for Astronomy, Zenodo, doi:10.5281/zenodo.3746095
  * Deil et al. (2017) Deil C., et al., 2017, in 35th International Cosmic Ray Conference (ICRC2017). p. 766 (arXiv:1709.01751)
  * Drury (1983) Drury L. O., 1983, Reports on Progress in Physics, 46, 973
  * Fermi Science Support Development Team (2019) Fermi Science Support Development Team 2019, Fermitools: Fermi Science Tools, Astrophysics Source Code Library, record ascl:1905.011 (ascl:1905.011)
  * Fujisawa et al. (2015) Fujisawa K., Yonekura Y., Sugiyama K., Horiuchi H., Hayashi T., Hachisuka K., Matsumoto N., Niinuma K., 2015, The Astronomer’s Telegram, 8286, 1
  * Garay et al. (2003) Garay G., Brooks K. J., Mardones D., Norris R. P., 2003, ApJ, 587, 739
  * H. E. S. S. Collaboration et al. (2022a) H. E. S. S. Collaboration et al., 2022a, Science, 376, 77
  * H. E. S. S. Collaboration et al. (2022b) H. E. S. S. Collaboration et al., 2022b, Science, 376, 77
  * Hirota et al. (2021) Hirota T., Cesaroni R., Moscadelli L., Sugiyama K., Burns R. A., Kim J., Sunada K., Yonekura Y., 2021, A&A, 647, A23
  * Kohno et al. (2022) Kohno M., et al., 2022, PASJ, 74, 545
  * Marti et al. (1993) Marti J., Rodriguez L. F., Reipurth B., 1993, ApJ, 416, 208
  * Marti et al. (1995) Marti J., Rodriguez L. F., Reipurth B., 1995, ApJ, 449, 184
  * Motch et al. (1997) Motch C., Haberl F., Dennerl K., Pakull M., Janot-Pacheco E., 1997, A&A, 323, 853
  * Mucciarelli et al. (2011) Mucciarelli P., Preibisch T., Zinnecker H., 2011, A&A, 533, A121
  * Munar-Adrover et al. (2011) Munar-Adrover P., Paredes J. M., Romero G. E., 2011, A&A, 530, A72
  * Obonyo et al. (2021) Obonyo W. O., Lumsden S. L., Hoare M. G., Kurtz S. E., Purser S. J. D., 2021, MNRAS, 501, 5197
  * Padovani et al. (2015) Padovani M., Hennebelle P., Marcowith A., Ferrière K., 2015, A&A, 582, L13
  * Pedaletti et al. (2014) Pedaletti G., de Oña Wilhelmi E., Torres D. F., 2014, A&A, 565, A118
  * Protheroe (1998) Protheroe R. J., 1998, pp 259–297
  * Richardson et al. (1985) Richardson K. J., White G. J., Gee G., Griffin M. J., Cunningham C. T., Ade P. A. R., Avery L. W., 1985, MNRAS, 216, 713
  * Stecklum et al. (2016) Stecklum B., Caratti o Garatti A., Cardenas M. C., Greiner J., Kruehler T., Klose S., Eisloeffel J., 2016, The Astronomer’s Telegram, 8732, 1
  * Su et al. (2019) Su Y., et al., 2019, ApJS, 240, 9
  * Szymczak et al. (2018) Szymczak M., Olech M., Wolak P., Gérard E., Bartkiewicz A., 2018, A&A, 617, A80
  * Wang et al. (2011) Wang Y., et al., 2011, A&A, 527, A32
  * Wang et al. (2017) Wang C., Yang J., Xu Y., Li F., Su Y., Zhang S., 2017, ApJS, 230, 5
  * Wood et al. (2017) Wood M., Caputo R., Charles E., Di Mauro M., Magill J., Perkins J. S., Fermi-LAT Collaboration 2017, in 35th International Cosmic Ray Conference (ICRC2017). p. 824 (arXiv:1707.09551)
  * Yan et al. (2022) Yan D.-H., Zhou J.-N., Zhang P.-F., 2022, Research in Astronomy and Astrophysics, 22, 025016
  * Yang et al. (2023) Yang R.-z., Li G.-X., Wilhelmi E. d. O., Cui Y.-D., Liu B., Aharonian F., 2023, Nature Astronomy,
  * Zabalza (2015) Zabalza V., 2015, in 34th International Cosmic Ray Conference (ICRC2015). p. 922 (arXiv:1509.03319), doi:10.22323/1.236.0922
  * Zinchenko et al. (2012) Zinchenko I., Liu S. Y., Su Y. N., Kurtz S., Ojha D. K., Samal M. R., Ghosh S. K., 2012, ApJ, 755, 177
  * Zinchenko et al. (2015) Zinchenko I., et al., 2015, ApJ, 810, 10