Probing the origin of cosmic rays in Cygnus Cocoon using ultrahigh-energy gamma-ray and neutrino observations

Recent ultrahigh energy gamma-ray observations by the HAWC up to 100 TeV and LHAASO observatories up to 1.4 PeV energies from the direction of Fermi-LAT 4FGL source 4FGL J2028.6+4110e (Cygnus Cocoon), are indicative of a hadronic origin over a leptonic process for their creation. The IceCube Neutrino Observatory has reported IceCube-201120A, a neutrino event coming from the same direction, suggesting that the Cygnus Cocoon may correspond to one of the most plausible sources of high-energy cosmic rays. The apparent relationship of the neutrino event with the observed ultra high energy gamma-rays from Cygnus Cocoon is investigated in this work to study if it can be explained consistently in hadronic interactions of accelerated cosmic rays with ambient matter. Our findings reveal that leptonic mechanisms, together with pure hadronic mechanisms, make a considerable contribution to the understanding of the total electromagnetic spectrum as well as the observed neutrino event. The estimate of expected muon neutrino events from the Cygnus cocoon agrees with the one muon neutrino event detected so far in IceCube multi-year observations. Thus, our results are indicative of the potential of the Cygnus Cocoon to be a galactic cosmic ray source capable of accelerating at least up to PeV energies.


INTRODUCTION
The Milky Way is known to accelerate cosmic rays with energies up to a few PeV (PeVatrons), however the origin of Galactic cosmic rays has yet to be proven (Hillas 1984;Berezinskii et al. 1990). Supernova remnants (SNRs) are commonly regarded as the most likely origins of Galactic cosmic rays since they are powerful and abundant enough to sustain the intensity of observed cosmic rays (Baade & Zwicky 1934;Blasi 2013). However, normal conditions make it difficult for SNR to accelerate particles to PeV energies (Bell et al. 2013;Aharonian 2013). Moreover, there is just no observational evidence to support SNRs as sources of hadrons with energy more than a few tens of TeV (Helder et al. 2012;Aharonian et al. 2019). Supernova explosions tend to cluster in space (within a few parsecs) and time because large OB stars (the progenitors of corecollapse supernovae) are formed in clusters and live short lives (within a few 10 5 yr) (Higdon & Lingenfelter 2005). As a result, Galactic cosmic rays with energy up to PeV are anticipated to be accelerated by overlapping shocks from SNRs and massive stellar winds (referred to as superbubbles (Tenorio-Tagle & Bodenheimer 1988)) generated around OB associations (Bykov & Fleishman 1992;Parizot, E. et al. 2004). The energy spectra and radial distribution of the calculated cosmic ray flux give evidence for particles accelerated to near-PeV energies in large star clusters .
The accelerated cosmic rays in the superbubbles may interact with ambient matter and radiation within the source, producing ultrahigh energy gamma-rays and neutrinos with energies up to PeV. Observations of neutrinos produced in association with ultrahigh energy gamma-rays would unambiguously identify superbubbles as a Galactic cosmic ray PeVatrons. The IceCube collaboration recently announced a candidate track-like neutrino event with an estimated energy of 154 TeV through the standard BRONZE alert procedure (Blaufuss et al. 2019) on November 20, 2020 (IceCube Collaboration 2020; Dzhappuev et al. 2021). The neutrino event is likely linked to an extended gamma-ray source Cygnus Cocoon which is a superbubble surrounding a region of OB2 massive star formation (IceCube Collaboration 2020). The Carpet−2 experiment reported a 3.1σ (post-trial) excess of atmospheric air showers from the same direction, consistent with a few months flare in photons over 300 TeV, in temporal correlation with the neutrino event (Dzhappuev et al. 2021). Implication of this observation has been discussed in detail in results section. This is the first evidence for a neutrino event being correlated with a Galactic object, despite being determined with considerable uncertainty (Dzhappuev et al. 2021).
The Fermi Gamma-Ray Space Telescope's Large Region Telescope (LAT) first identified Cygnus Cocoon emitting hard, multi-GeV gamma-rays in the nearby star-forming area known as Cygnus X (Ackermann et al. 2011). The ARGO experiment first identified it at TeV energies (Bartoli et al. 2014). The High Altitude Water Cherenkov (HAWC) observatory has reported observations of 1 − 100 TeV gamma-rays originating from the Cygnus Cocoon, which may be represented by a power-law below 10 TeV and exhibits spectrum softening around 10 TeV (Abeysekara et al. 2020;Abeysekara et al. 2021). The LHAASO collaboration recently reported the discovery of ultrahigh energy photons with energies up to 1.4 PeV from this location, indicating that the spectrum can extend up to ∼ 1 PeV (Cao et al. 2021;Li 2021). Amenomori et al. (2021) revealed that gamma-ray sources in the Cygnus area contribute significantly to the Galactic-plane diffuse gamma radiation above 400 TeV. These findings imply that the Cygnus Cocoon may be acting as galactic hadronic Pe-Vatron, and could provide significant evidence for understanding the knee of the observed cosmic ray energy spectrum.
The star-forming regions such as Cygnus constellation have been proposed as potential sites for cosmic-ray acceleration, as well as gamma-ray and neutrino production . It has been suggested that the high-energy neutrino flux from Cygnus Cocoon will be close to the IceCube sensitivity (Yoast-Hull et al. 2017). The detected gamma-ray flux might be both leptonic and hadronic in nature. The HAWC collaboration reported that the observed 1-100 TeV gamma-rays from Cygnus Cocoon were unlikely to be explained by a single electron population emitting gamma-rays from GeV to the highest energy via inverse-Compton emission without its synchrotron radiation exceeding the flux limits set by radio and X-ray studies (Abeysekara et al. 2021). Therefore, a significant contribution in the gamma-ray spectrum above a few TeV likely to have hadronic origin for their production. However, the appearance of a cut-off or a break in the measured gamma-ray spectrum at a few TeV is thought to be due to either cosmic ray leakage from the Cocoon or a cut-off in the cosmic ray spectrum injected from the source (Abeysekara et al. 2021). Again, the most recent LHASSO observations of ultrahigh-energy gamma-rays up to 1.4 PeV from the Cygnus Cocoon significantly disfavored a leptonic origin for their formation and strongly imply acceleration of cosmic rays at energies greater than PeV (Cao et al. 2021). Therefore, the apparent association of the observed neutrino event IceCube-201120A with the measured ultra high energy gamma-rays up to highest energies demands new explanation for their productions.
Under such circumstances, we would like to examine the apparent link of the neutrino event IceCube-201120A with the observed ultra high energy gammarays from Cygnus Cocoon to study whether it can be explained consistently in the framework of hadronic interactions of accelerated cosmic rays with ambient matter. We will also investigate the possibility of a leptonic origin contribution to the total gamma-ray spectra, in addition to a hadronic origin, taking into account Ice-Cube's non-detection of multiple neutrino events. We would also like to inspect the maximum energy that a cosmic ray particle can achieve in the Cygnus Cocoon.
The following is the article's structure: The next section describes the method for evaluating gammaray and neutrino fluxes from Cygnus Cocoon. Section 3 shows numerical estimates of the fluxes of multiwavelength electromagnetic (EM) spectral energy distribution (SED) and high-energy neutrinos produced by the Cygnus Cocoon. The results are also discussed in the same section. Finally, we will conclude in the section 4.

METHODOLOGY
It can be generally assumed that the electrons and protons are co-accelerated via the diffusive shock acceleration mechanism in the interacting winds created by the collective activity of massive stars in the Cygnus Cocoon. We can assume that a fraction of the stellar wind energy, η e , may be used to accelerate electrons and a fraction, η p , can be used to accelerate hadrons (Bednarek 2007).
In this case, we consider a broken power-law energy distribution of shock-accelerated electrons with spectral indices α 1 and α 2 before and after the spectral break at Lorentz factor γ b , as illustrated below (Katarzyński et al. 2001) (1) where K e represents the normalization constant which is related with the available stellar wind power (L w ) as given bellow where γ e = E e /m e c 2 represents the Lorentz factor of electrons of energy E e , r is the radius of the Cygnus Cocoon, and v w is a typical stellar wind velocity. For a radiative cooling break in a uniform magnetic field B, the electron distribution breaks in its index by one power (i.e, ∆α = α 2 − α 1 ≈ 1) above the spectral break Lorentz factor γ ′ b (Longair 1994). The low-energy component of the EM SED extending from radio to x-ray energies generated in the Cygnus Cocoon is represented by synchrotron radiation of primary accelerated electrons, which is estimated here using the methodology given by Böttcher et al. (2013). The inverse Compton (IC) scattering of primary accelerated electrons with target photons contributes significantly to the observed EM spectrum in the MeV to GeV ranges and may be estimated using the formulas presented in Blumenthal & Gould (1970). Here, we consider four target radiation fields following Ackermann et al. (2011) for gamma-ray generation via IC scattering, including synchrotron photons, strong stellar light fields around Cyg OB2 and NGC 6910 (a star cluster in the neighborhood of OB2), and a more diffuse dust radiation field over the whole Cocoon. The Bremsstrahlung scattering of primary accelerated electrons with ambient matter of density n H is found to explain the observed gamma-ray spectrum in GeV to few TeV energy band and may be estimated by following Blumenthal & Gould (1970).
The cosmic ray (protons) production spectrum is also expected to follow a power law (Malkov & Drury 2001) where E p denotes the energy of the cosmic ray proton, α p is the spectral index, and the proportionality constant is K p , which may be derived by using the fraction of stellar wind energy carried by cosmic ray protons as follows where E p,min and E p,max denotes the minimum and maximum energies of accelerated cosmic ray protons, respectively. The interaction of shock-accelerated cosmic ray protons with ambient matter (protons) of density n H can explain the high energy gamma-ray emission from the Cygnus Cocoon. Such hadronic (pp) interaction produces neutral and charged pions, which decay to create high-energy gamma-rays and neutrinos, respectively. To estimate the high energy gamma rays and neutrino emissivities (Q pp γ & Q pp ν respectively) produced in hadronic interactions (pp) in the Cygnus Cocoon, we follow Refs. (Kelner et al. 2006;Banik & Bhadra 2017;Banik & Bhadra 2019). The associated differential flux of high-energy gamma-rays and muon neutrinos reaching Earth from the Cygnus Cocoon may be represented as where ζ is a constant equal to 1 for gamma rays and equal to 1/3 for muon neutrinos due to neutrino oscillation, V = 4 3 πr 3 represents the volume of the emission region, and d is the distance between the Cygnus Cocoon and the Earth. The number of expected muon neutrino events in the IceCube detector in time t may be computed using the following formula where A ef f is the effective area of the IceCube detector (IC86) at the declination of the source (IceCube Collaboration 2021). We can choose E ν,min ≈ 30 TeV which is in good agreement with the effective energy threshold of the IceCube detector for astrophysical neutrinos (Taboada, Ignacio 2016).

RESULTS AND DISCUSSION
The Cygnus Cocoon is situated at RA 307.17 deg and Dec 41.17 deg (J2000) (IceCube Collaboration 2020). It has an angular size of around 2.1 • , corresponding to a radius of r = 55 pc at a distance of d = 1.4 kpc from Earth (Abeysekara et al. 2021). The Cygnus Cocoon comprises of two star clusters, Cyg OB2 and NGC 6910, with total wind power estimates of (2 − 3) ×10 38 erg/s and (1−1.5) ×10 36 erg/s, respectively (Ackermann et al. 2011). Here, we consider a typical stellar wind with velocity v w = 10 3 km/s (Ackermann et al. 2011) and wind power L w = 3 × 10 38 erg/s to estimate stellar wind energy density.
The IceCube neutrino observatory has reported a neutrino event IceCube-201120A, which is likely to be associated to Cygnus Cocoon. The Carpet−2 experiment observed an excess of gamma-ray events consistent with a few months flare in photons above 300 TeV from the direction of the Cygnus region, in temporal and spatial coincidence with the IceCube neutrino alert (Dzhappuev et al. 2021). As the gamma-ray emission from the cocoon (the radius ∼ 55 pc) cannot vary on a time scale of months, this neutrino event is more likely to be due to any compact source within the cocoon rather than the diffused emission. Therefore, if this event is truly related with the few months gamma-ray flare, it cannot be used to limit hadronic gamma-ray emission from the cocoon. However, the neutrino event may not be linked to such a flare because no excess GeV-TeV gamma-ray flux was measured from the Cygnus region during the neutrino arrival period by Fermi-Lat and HAWC (Garrappa et al. 2020;Ayala & HAWC Collaboration 2020).
Recently, the LHAASO observatory detected gammaray flux of 0.54(0.10) CU at 100 TeV (CU, the Crab Nebula flux at 100 TeV; 1 CU = 6.1 × 10 −17 photons TeV −1 cm −2 s −1 ) from the direction of the source LHAASO J2032+4102 (Cygnus Cocoon) when only half of its KM2A detectors were operational (Cao et al. 2021). The LHAASO observatory found 45 on-source events (with 6.7 number of background events) with energies above 100 TeV up to 1.4 PeV during an exposure period of 2648.2 hr from the direction of Cygnus Cocoon (Cao et al. 2021). We estimated the corresponding gamma-ray flux seen by LHAASO from Cygnus Cocoon by assuming a power-law gamma-ray spectrum of f γ = N 0 E −Γ with a photon index of Γ ≈ 2.7 in the energy range of 100 TeV to 1.4 PeV. The normalization constant N 0 of the observed gamma-ray flux from Cygnus Cocoon by LHAASO may be computed as (Aharonian et al. 2020) where S γ represents the number of gamma-ray signal event in LHAASO detector, A γ ef f denotes the effective area when only half of its KM2A detectors were operational , and T ex represents corresponding exposure time (Cao et al. 2021) for Cygnus Cocoon. Here, ǫ = 0.68 is the fraction of observed event counts within the angular resolution of the instrument (He et al. 2019). The combined gamma-ray spectra from GeV to highest energies (1.4 PeV) from Cygnus Cocoon as observed by Fermi-LAT (Ackermann et al. 2011;Astiasarain et al. 2021), ARGO (Bartoli et al. 2014), HAWC (Abeysekara et al. 2021), and LHAASO (Cao et al. 2021;Li 2021) observatories indicate hadronic origin for their generation (Abeysekara et al. 2021;Cao et al. 2021).
The likelihood of a leptonic origin contribution to the overall gamma-ray spectra, in addition to a hadronic origin, is examined below, taking into consideration Ice-Cube's non-detection of multiple neutrino events from Cygnus Cocoon.

Pure hadronic origin
We have estimated the gamma-ray flux produced in the hadronic interaction (pp) of a single population of accelerated cosmic rays (protons) with the ambient protons within the astronomical object Cygnus Cocoon.
Because the Cygnus area contains a massive molecular cloud complex with a total mass of 8 × 10 6 M ⊙ (Ackermann, M. et al. 2012), the interstellar gas density should be more than 10 cm −3 (Bartoli et al. 2014). To explain the observed EM SED, we choose an ambient matter density of n H = 30 cm −3 in the area, as suggested by HI and HII observations (Abeysekara et al. 2021). Here, we adopt a magnetic field of B = 20 µG as inferred from pressure balance with the gas throughout the Cygnus Cocoon region (Ackermann et al. 2011). The red dashed and magenta dash-single-dotted lines represent the acceleration timescales with ξp = 9 × 10 −5 and 5 × 10 −7 respectively. The blue dash-double-dotted, and green dash-triple-dotted lines denote the energy loss timescale for protons in pp−interaction and diffusion timescale respectively. The black continuous line indicates the lifetime of the Cocoon. The points, as indicated by the cyan and brown arrows, represent the maximum achievable energy by a cosmic ray proton with ξp = 9 × 10 −5 and 5 × 10 −7 respectively.
We compared the acceleration time-scale of a proton with its energy loss time-scale in pp−interaction, diffusion timescale, and the age of the Cocoon (t age ∼ 2 Myr) to understand the maximum achievable energy by a cosmic ray particle within the Cocoon. The timescale of acceleration of cosmic ray protons can be represented as (t acc = Ep ξpeBc ), where ξ p (≤ 1) is the proton acceleration coefficient. The energy loss timescale for protons in pp−interaction can be written as t pp = 1 kppσppnH c , where k pp = 0.45 and σ pp represent the in-elasticity (Gaisser 1990) and the interaction cross section (Kelner et al. 2006), respectively. The diffusion timescale can be represented by t dif f = r 2 D dif f (Bednarek 2007), where D dif f = D 0 ( Ep 10 GeV ) δ (Giuliani, A. et al. 2010;Berezinskii et al. 1990) denotes the diffusion coefficient of accelerated protons. We may choose D 0 = 1.2 × 10 27 cm 2 /s as the diffusion coefficient at 10 GeV energy since the dense gaseous medium has a slower diffusion than the galactic medium (≈ 10 28 cm 2 /s in our Galactic medium) (Berezinskii et al. 1990;Aharonian & Atoyan 1996). We have taken into account δ = 0.33, as recently found by the Alpha Magnetic Spectrometer (AMS-02) when measuring the boron to carbon flux ratio in cosmic rays (Aguilar et al. 2016). The aforementioned timescales of relativistic protons as functions of proton energy are displayed in Fig. 1 (also see Bednarek (2007)). By comparing the acceleration time-scale of a proton with its diffusion timescale (e.g. Bednarek (2007)), we found that hadrons can be accelerated up to 5 × 10 15 eV (or 10 14 eV) within the cocoon with an acceleration coefficient of ξ p = 9 × 10 −5 (or 5 × 10 −7 ). To match the observed gamma-ray spectrum, we first consider a primary cosmic-ray production spectrum obtained by assuming an acceleration efficiency of η p = 100% of cosmic ray protons with a spectral index of α p = −2.6, and the maximum energy of E p,max = 5×10 15 eV. According to our findings, an estimated gamma-ray flux based on a single power-law distribution of accelerated cosmic rays cannot properly explain the observed overall gamma-ray spectrum up to energies of 1.4 PeV. This is mostly due to a spectral break in the gamma-ray spectrum around 10 TeV energy, which has been seen in both ARGO and HAWC detector studies.
When we consider η p = 1.6%, i.e. the portion of the stellar wind energy carried by cosmic ray protons with a spectral index of α p = −2.35 and a maximum achievable energy of E p,max = 10 14 eV, the observed gamma-ray spectrum can be reproduced well from GeV to few tens of TeV energy range. However, it was unable to explain the reported PeV gamma-ray flux by the LHAASO detector. The estimated differential gammaray flux reaching the Earth from the Cygnus Cocoon for two stated scenarios are displayed in the Fig. 2 along with the observations.
In the next section, we have investigated whether viable leptonic processes, in conjunction with pure hadronic mechanisms, can explain the observed spectrum of gamma-rays, as well as the reported neutrino event from Cygnus Cocoon.

Lepto-hadronic origin
The source is thought to accelerate both electrons and protons at the same time. The acceleration timescale of electrons is expressed similarly to that of protons, but with ξ e ≤ 1 as the acceleration coefficient. To estimate the diffusion timescale, we use D dif f = 1.2 × 10 27 ( Ee 10 GeV ) 0.33 cm 2 /s as the diffusion coefficient for electrons, which is the same as that for protons at 10 GeV energy. The advection timescale can be estimated as t adv = r/v w (Bednarek & Sitarek 2007). The maximum electron Lorentz factor (γ e,max ) was calculated by matching the acceleration timescale with the synchrotron energy loss timescale (t cool = 3mc 4uB σT γe ) and was found to be 8 × 10 7 with an acceleration coefficient of ξ e = 10 −5 . Here, u B = B 2 8π represents the magnetic field energy density and σ T is the Thomson scattering cross section. The synchrotron cooling timescale begins to take precedence over the diffusive timescale, i.e. the average time spent by electrons with energy E e inside the cocoon region, at Lorentz factor 1.7 × 10 5 as shown in Fig. 3, which can be regarded as the spectral break Lorentz factor γ b . The red dashed and blue dash-single-dotted lines represent the acceleration timescale of electron with ξe = 10 −5 and synchrotron cooling timescale respectively. The green dashdouble-dotted, and magenta dash-triple-dotted lines denote the diffusion, and advection timescales respectively. The points, as indicated by the cyan and brown arrows, represent the maximum achievable Lorentz factor, and spectral break Lorentz factor of a electron, respectively.
We consider a broken power-law distribution of accelerated primary electrons with spectral indices of α 1 = 2.1 and α 2 = 3.1 before and after the spectral break at the Lorentz factor γ b = 1.7 × 10 5 , derived by assuming a fraction of wind energy η e = 9% carried by the electrons. The synchrotron emission of the accelerated electrons has been computed and compared to multi-wavelength observations of the Cygnus Cocoon from radio to X-ray energies. The IC scattering of primary relativistic electrons with synchrotron photons co-moving within the source is found to have no substantial contribution to the observed EM spectrum from the source at energies ranging from MeV to TeV. A significant contribution to the observed EM spectrum in the MeV to few tens of TeV energy range is found to be produced by IC scattering of primary relativistic electrons with strong star light fields around Cyg OB2 and NGC 6910 (a star cluster in the neighborhood of OB2), and a more diffuse dust radiation field over the entire Cocoon. Furthermore, we have found that the gamma-ray flux produced by Bremsstrahlung scattering of primary relativistic electrons in ambient matter with a density of n H = 30 cm −3 may contribute significantly to the measured gamma-ray energies ranging from MeV to TeV from the source.
The interactions of relativistic primary cosmic rays with ambient protons in the source can contribute significantly to EM SED above 100 TeV energies, as detected by the LHAASO observatory, as well as generate high energy neutrinos. The required fraction of star wind energy carried by the accelerated primary protons is determined to be η p = 8% with the best fitting spectral slope of α p = −2.4 and the maximum attainable energy E p,max ≃ 5 × 10 15 eV. The left panel of Fig. 4 shows the estimated differential gamma-ray spectra escaping from the Cygnus Cocoon along with the different satellite and ground-based observational data. Using Eq. 7 and the model-estimated gamma-ray flux from Cygnus Cocoon, we again computed the predicted gamma-ray signal events in the LHAASO detector from the source and found that they were compatible with the observations.
The estimated neutrino flux from Cygnus Cocoon reaching at Earth is displayed in the right panel of the Fig. 4. The expected muon neutrino event in the Ice-Cube detector from Cygnus Cocoon is estimated to be roughly N νµ = 0.65 event above 30 TeV energy in 10 years using Eq. (6). Because a possible contribution to the event rates due to interactions of tau-neutrinos that create muons with a branching ratio of 17.7% was not addressed, the estimated number of neutrino events is conservative (Ansoldi et al. 2018;Banik et al. 2020). As a result, the total muon like neutrino events may be computed as N like µ = N µ + 17.7% × N µ , resulting in a 0.77 event. Our estimate of expected muon neutrino events is consistent with the only one muon neutrino event reported so far from the cocoon direction in IceCube multi-year observations. The model fitted parameters are displayed in Table 1.

CONCLUSION
The observed 1-100 TeV gamma-rays from Cygnus Cocoon is unlikely to be explained by a single electron population emitting gamma-rays from GeV to the highest energy via inverse-Compton, and Bremsstrahlung emission without exceeding the flux limits established by radio and X-ray studies. Our findings show that the combined gamma-ray spectra from GeV to maximum energies (1.4 PeV) from Cygnus Cocoon as reported by Fermi-LAT, ARGO, HAWC, and LHAASO observatories cannot be explained by pure hadronic (pp) interactions of relativistic cosmic rays with ambient matter. Our results suggest that leptonic processes, in combination with pure hadronic mechanisms, are necessary to consistently represent the complete electromagnetic  Figure 4. Left: The estimated differential EM SED reaching earth from the Cygnus Cocoon in the Lepto-hadronic origin scenario. The red dotted and green small-dashed lines represent the EM spectrum produced by synchrotron emission and Bremsstrahlung scattering of primary relativistic electrons in ambient matter, respectively. The EM spectrum produced by IC scattering of relativistic electrons with star light fields around NGC 6910 and Cyg OB2, as well as a dust radiation field, is denoted by the gray long-dash-single-dotted, magenta long-dash-double-dotted, and violet long-dash-triple-dotted lines, respectively. The gamma-ray flux created in pp interactions between relativistic protons and ambient protons is indicated by the brown smalldash-single-dotted line. The black continuous line shows the estimated overall differential multiwavelength EM SED coming from the Cygnus Cocoon. The cyan very-long-dash-single-dotted, and orange very-long-dashed lines indicate the detection sensitivity of the LHAASO, and e-ASTROGAM detectors for one years of observation, respectively. The blue long dashed line represents our estimated gamma-ray flux limit detected by LHAASO. Right: The estimated corresponding all flavor neutrino flux reaching the Earth from the Cygnus Cocoon. spectrum. Particularly, the detected gamma-ray flux in sub-PeV energies by LHAASO is found to be best explained by hadronic interaction of cosmic rays, which originated in the Cygnus Cocoon, with ambient matter. The single muon neutrino event detected so far from the cocoon direction in IceCube multi-year data agrees with our estimate of expected muon neutrino events. Thus, the Cygnus Cocoon might be one of the long-suspected galactic PeVatrons, capable of accelerating cosmic rays with energies at-least upto few PeV, providing strong evidence for the origin of knee in the observed cosmic ray energy spectrum. Future gamma ray telescopes with better sensitivity than current generation gamma ray telescopes, such as e-ASTROGAM (de Angelis et al.