© The Authors 2025

1 Introduction

The Cygnus region is the most massive star-forming complex in the Local Arm. It displays rich multi-wavelength emission features from radio to ultra-high-energy (> 100 TeV) γ-rays, revealing a highly structured environment shaped by the feedback of numerous stellar objects. Unfortunately, the region lies in a high extinction patch of the sky, hidden behind a dark lane of gas and dust called the Cygnus Rift (see Reipurth & Schneider 2008, for a review). At a galactic longitude of about 80°, it is also seen tangentially to the direction of galactic rotation, which complicates distance estimations and makes it difficult to disentangle superimposed features along the line of sight. Despite decades of extensive studies, the observations remain poorly understood, in particular the nature and origin of the structures visible on large scales.

Cygnus was first detected as a diffuse radio source by Piddington & Minnett (1952), who named it Cygnus-X (‘X’ stands for ‘extended’). In addition to the diffuse emission, Downes & Rinehart (1966) identified 20 individual HII regions, later assigned as DR4–DR23, which appear as prominent radio hotspots. The only source of non-thermal (synchrotron) radio emission is the supernova remnant γ-Cygni (Xu et al. 2013a). The diffuse thermal component comes from material ionised by massive stars, especially the stellar association Cygnus OB2, which is thought to be responsible for the diffuse emission in a Strömgren sphere of 2^° in radius. In addition to the diffuse emission, the compact radio hotspots coincident with dense clumps originate from massive stars and protostars within the surrounding clusters and star-forming clouds (Knödlseder et al. 2002; Knödlseder 2004). Enclosing the diffuse radio emission is a large (18° × 13°) elliptic ring of soft X-rays (0.5–2 keV, Cash et al. 1980). This X-ray shell, historically referred to as the ‘Cygnus Superbubble’, is thought to have been created by the collective feedback of massive stars and past supernovae in Cygnus X over the past tens of millions of years (Kimura et al. 2013; Bluem et al. 2020), though the possibility of a chance coincidence of X-ray emission along the line of sight is not excluded (Uyaniker et al. 2001). Another hint of large-scale feedback in Cygnus is the presence of a ring of Hα filaments (Dickel et al. 1969) coincident with the X-ray shell. These observations are consistent with the existence of an ionised low-density cavity of radius 100–150 pc roughly centred on Cygnus OB2, which could act as the source of ionising radiation and mechanical power. This notion is consistent with the current stellar-wind power in Cygnus OB2 (2−3 × 10³⁸ erg s⁻¹, Vieu et al. 2024b) integrated over several Myr of stellar activity. Detection of diffuse X-ray emission in the close vicinity of Cygnus OB2 at 1–3 keV suggests| that stellar winds are indeed interacting close to the centre of the association (Albacete-Colombo et al. 2023). Despite its relevance, the relation between the small-scale stellar feedback and the large-scale features is still not settled and has received scarce attention in the recent literature.

Stellar feedback in such a complex environment is not expected to produce spherically symmetric bubbles but rather intricate, asymmetric structures. Critical factors include the wind properties of the most massive stars, the inhomogeneous interstellar medium, and the past activity of supernovae and older star clusters in the region (e.g. NGC6913 and NGC6910; see Almeida et al. 2023). Mid-infrared observations of Cygnus X reveal a rich morphology in the warm dust emission with bright photodissociation regions, molecular ridges, extended pillars, filaments, and parsec-scale globules (Schneider et al. 2006). Most likely, these structures were heated and shaped by the Cygnus OB2 association. Far-infrared imaging (Schneider et al. 2016) shows a complex web of colder molecular clouds interacting with multiple HII regions, including nascent star clusters close to Cygnus OB2, which hints at triggered star-formation events. These infrared observations correlate well with the CO intensity, which traces molecular hydrogen. The molecular complex is usually divided between the northern (l > 80°) and the southern (l < 80°) clouds. In Cygnus North lies the well-studied dense star-forming molecular ridge DR21 (e.g. Hennemann et al. 2012; Bonne et al. 2023). In Cygnus South, the CO emission is more diffuse yet tends to coincide with bright mid-infrared emission, particularly at the location of the molecular cloud L889 and radio hotspots (e.g. DR9, DR12, DR13, and DR15; see Schneider et al. 2006). It is possible that the molecular material was shaped by interacting superbubble cavities blown by Cygnus OB2, NGC 6913 (or the more extended Cygnus OB9 association), and NGC 6910 (e.g. Schneider et al. 2007), although the difficulty in measuring accurate distances to both cloud material and star clusters complicates any such interpretation. Recent progress has been made by Zhang et al. (2024), who applied an extinction jump method to high resolution CO data and concluded that the molecular gas is organised in multiple layers between 0.8 and 2 kpc from the Sun. In particular, the method successfully disentangles the foreground rift component, which appears to extend out to approximately 1.1 kpc and seems to contain more gas than previously anticipated in earlier work by Schneider et al. (2006). An intermediate layer at 1.3 kpc, with counterparts both in Cygnus North (W75N and DR17, see also Gottschalk et al. 2012) and Cygnus South, also appears in the foreground of the active star-forming complex, the latter being located around 1.6–1.8 kpc. Unfortunately the extinction is too high in Cygnus to track extinction jumps out to 2 kpc. It is not presently possible to precisely disentangle the foreground clouds from the molecular web interacting with the star-forming complex at larger distances.

Finally, identifying coherent stellar groups is not straightforward. The stellar population is scattered and diverse, with a large age-spread and possible past supernova events (e.g. Reipurth & Schneider 2008; Comerón & Pasquali 2012; Comerón et al. 2020; Berlanas et al. 2020). Modern clustering algorithms do not identify Cygnus OB1, OB8, and OB9 as genuine associations (Quintana & Wright 2021), and Gaia measurements indicate that even Cygnus OB2, the largest and most powerful group of stars, is not a single coherent group but a superposition of several groups along the line of sight (Berlanas et al. 2019; Orellana et al. 2021), the most massive being located at about 1.6–1.7 kpc. The latter is what we refer to as ‘Cygnus OB2’ in the remainder of the present work (see Section 2.1).

It might be unsurprising that such a complex region with a long history of star-forming activity is a bright source in the γ-ray sky. Extended GeV sources in Cygnus were first detected by Thompson et al. (1995) and later associated with the γ-Cygni supernova remnant (Brazier et al. 1996) and the Dragonfly pulsar wind nebula (around PSR J2021+3651, Roberts et al. 2002). An extended TeV source coincident with Cygnus OB2 was discovered by HEGRA Collaboration (2002) and Abdo et al. (2007) and later associated with the pulsar wind nebula powered by PSR J2032+4127 (Camilo et al. 2009).

In addition to these three identified γ-ray sources, the Fermi-LAT collaboration discovered an extended diffuse γ-ray source coincident with the radio emission in Cygnus X, filling a region delimited by bright infrared photodissociation regions (Ackermann et al. 2011). They demonstrated that in-situ particle acceleration was necessary to account for the observation. This breakthrough revived the high-energy community’s interest in star-forming regions as cosmic-ray accelerators and was followed up with detections by ARGO-YBJ Collaboration (2014), Abeysekara et al. (2018), Tibet AS_γ Collaboration (2021), Abeysekara et al. (2021), and LHAASO Collaboration (2024a). The centroid of the diffuse emission lies at the heart of Cygnus X, is approximately coincident with the position of Cygnus OB2 below 1 TeV, and is loosely correlated with the northern and southern molecular clouds at ultra-high-energy. The spectrum is hard in the GeV band (dN/dE ~ E⁻²) but turns over around 1 TeV (dN/dE ~ E⁻³). The steep component measured by LHAASO Collaboration (2024a) extends up to 2 PeV with a flux substantially above the estimate for the galactic diffuse background, which suggests the presence of an extreme particle accelerator in the region. This accelerator is unlikely to be related to PSR J2032+4127 or γ-Cygni at their current stages of evolution, as the spectra of these sources cut off at energies ≲ 100 TeV, and their integrated flux is subdominant with respect to the extended diffuse component (Fleischhack 2019; Abeysekara et al. 2021; MAGIC Collaboration 2023; LHAASO Collaboration 2024a).

The nature of candidate PeV accelerators in Cygnus has been debated in recent years. Bykov & Kalyashova (2022) proposed a scenario of particle acceleration in supersonic turbulence excited by the wind feedback in Cygnus OB2. They modeled a cavity of radius 55 pc filled with a turbulent plasma of average velocity 1500 km s⁻¹, though this requires an energy budget that cannot be supplied by the massive stars in Cygnus OB2 alone for any reasonable plasma density. Such a scenario also fails to account for the ultra-high-energy detection by LHAASO because of fundamental limitations of second-order acceleration mechanisms (Vieu et al. 2022b, 2024a). Menchiari et al. (2024) postulated the existence of a large-scale wind termination shock around Cygnus OB2 and applied a particle acceleration and transport model in spherical symmetry (Morlino et al. 2021), showing good agreement with Fermi-LAT and HAWC data. This mechanism, however, is unlikely to produce multi-PeV protons (see Sect. 2), and it has been shown that Cygnus OB2, as a loose association, is not compact enough to produce a large-scale cluster wind (Vieu et al. 2024b). In addition, hadronic models for γ-ray emission in star-forming environments generally suffer from the fundamental caveat that they simultaneously need i) a low-density environment, which allows for the persistence of strong shocks and inside which particles are confined, and ii) a high-density environment for efficient proton-proton interactions. This is why leptonic scenarios were previously favoured in star-forming regions (see Härer et al. 2023, for the case of Westerlund 1). A leptonic origin for the γ-ray emission in Cygnus was excluded by Abeysekara et al. (2021), who argue that in a magnetic field of 20 µG and a density of 30 cm⁻³, one would expect a detection of synchrotron and bremsstrahlung radiation beyond the limits set by radio and X-ray observations (see also Guevel et al. 2023). These values for the magnetic field and gas density are however far in excess of what one would expect in a region dominated by strong stellar feedback. Magnetohydrodynamic simulations of star clusters show that the magnetic field drops rapidly in the cluster wind and is no more than a few µG in the wind-blown superbubble (Härer et al. 2025). Hydrodynamic simulations, along with analytic theory (Weaver et al. 1977) predict a density of at most a few 0.1 cm⁻³ within a few tens of parsecs from Cygnus OB2. The values quoted by the HAWC collaboration are taken from Ackermann et al. (2011), where they were used only to estimate the properties of the Cygnus molecular clouds, not those of the ionised gas around Cygnus OB2. Unfortunately, the properties of this gas remain unknown – though there is a strong theoretical argument in favour of a low-density cavity excavated by winds from massive stars in Cygnus OB2 (e.g. Bluem et al. 2020; Vieu et al. 2024b), in addition to the observational hints outlined above. This lack of gas promotes consideration for a leptonic origin of the diffuse γ-ray emission, especially given the strong photon field in the vicinity of Cygnus OB2 (see Schneider et al. 2016, and Sect. 5.1). Such models (e.g. Banik & Ghosh 2022; Astiasarain et al. 2023), however, cannot account for the ultra-high-energy component of the spectrum since the cooling time of PeV electrons is much shorter than the acceleration time at a non-relativistic shock under usual circumstances.

The aim of the present work is to clarify the origin of the broadband γ-ray emission in Cygnus in light of the current multi-wavelength picture. In Sect. 2, we start by reviewing the stellar objects observed in the region and conclude that all are disfavoured as the source of the non-thermal PeV radiation at their current stages of evolution. In Sect. 3, we make the case for a past supernova event in Cygnus OB2 and explore this scenario by performing a detailed hydrodynamic simulation of a supernova remnant in a Cygnus-like stellar association. Section 4 is devoted to a detailed modelling of the gas and disentangling the foreground components. We then proceed with the building of a lepto-hadronic model to show how the diffuse γ-ray emission can be produced by the joint action of stellar winds and a past supernova. Section 5 introduces the methods. The supernova remnant model is then discussed in Sect. 6 and the wind model in Sect. 7. We present our conclusions in Sect. 8.

2 Candidate particle accelerators in Cygnus

The Cygnus region hosts a diverse population of powerful objects, including massive star clusters and associations (e.g. Comerón et al. 2020; Quintana & Wright 2021), Wolf-Rayet stars, pulsars (e.g. Lyne et al. 2015; Jin et al. 2023), supernova remnants (e.g. Leahy et al. 2013) and microquasars (e.g. MAGIC Collaboration 2010). In the following section we review the extreme objects in Cygnus and discuss possible acceleration mechanisms.

2.1 Cygnus OB2

The most notable stellar group in Cygnus is the Cygnus OB2 association, which is by far the most massive and lies roughly at the centre of the region of diffuse radio emission (‘Cygnus X’). Other OB associations (Quintana & Wright 2021) or young star clusters (e.g. Berkeley 86, NGC6913, NGC6910) are substantially less massive and most likely older according to estimates by Cantat-Gaudin et al. (2020). Hosting a rich population of more than 70 O stars and three Wolf-Rayet stars, Cygnus OB2 has a mechanical wind luminosity of 2−3 × 10³⁸ erg s⁻¹ (see Vieu et al. 2024b for a review of the main properties of the association). It is worth noting that Cygnus OB2 is not a compact star cluster but a rather loose OB association, with a half-mass radius of about 15 pc (Wright et al. 2015). As such, the stellar population is heterogeneous and likely not coeval, with star ages ranging between 3 and 7 Myr and hints of multiple star-forming episodes (e.g. Berlanas et al. 2020). Part of the heterogeneity is due to the superposition of several stellar groups along the line-of-sight at about 1.3, 1.5 and 1.7 kpc (Orellana et al. 2021). This complication can also explain some earlier confusion in distance estimates (see Appendix A). In light of Gaia DR2, it seems that the main group is located at about 1.6–1.7 kpc, the uncertainty being mostly due to the Gaia zero-point offset (e.g. Lindegren et al. 2018). We stress that the distance of 1.4 kpc inferred by Rygl et al. (2012) employing maser measurements, which is frequently quoted in the high-energy community, is no longer believed to be the distance to Cygnus OB2, but rather the distance to the molecular complex in Cygnus North, which lies in the foreground of the association.

Fig. 1 Powerful objects in the Cygnus region. Stars represent Wolf-Rayet stars from Rosslowe & Crowther (2015, http://pacrowt.her.staff.shef.ac.uk/WRcat/index.php). Diamonds indicate pulsars. The distance to PSR J2032+4127 has been revised to 1.76 kpc by association with its Be star companion MT91 213. Rings show the location of supernova remnants from Green (2024, 2025, https://www.mrao.cam.ac.uk/surveys/snrs). Different methods give distances to γ-Cygni between 0.98 to 2.3 kpc. The circles labelled ‘10 pc’ and ‘100 pc’ are centred on Cygnus OB2, assuming a distance of 1.65 kpc to the association.

2.2 Other objects

A map showing the location of the powerful objects in the Cygnus region is provided in Fig. 1. Tens of Wolf-Rayet stars are scattered throughout the region, either clustered or isolated (Rosslowe & Crowther 2015). The most powerful are WR144 (most likely a member of Cygnus OB2), WR142 (likely a member of Berkeley 87) and WR147, each one having a mechanical power of about 10³⁸ erg s⁻¹ (Sander et al. 2019; Hamann et al. 2019). In addition, over 20 red supergiants are observed (Comerón et al. 2020). Three pulsars have been detected in the region: PSR J2021+3651 (spin-down power ~3 × 10³⁶ erg s⁻¹), which powers the Dragonfly pulsar wind nebula, PSR J2021+4026 (spin-down power ~10³⁵ erg s⁻¹), and finally PSR J2032+4127 (spin-down power ~10³⁵ erg s⁻¹), which likely lies within the Cygnus OB2 association and suggests past supernova activity at the heart of the star-forming complex. However, only two supernova remnants are observed in the 1–3 kpc range in the direction of Cygnus: the bright source γ-Cygni and the supernova remnant W63.

In the following we address whether any of these objects could potentially accelerate particles beyond 1 PeV, as required to explain the recent LHAASO detection of ultra-high-energy (> 100 TeV) photons from the region. A contribution from PSR J2021+3651 is excluded as the Dragonfly nebula is masked from the LHAASO analysis. A source coincident with PSR J2032+4127 in Cygnus OB2 is detected as a subdominant component with a low-energy (30 TeV) cut-off and is therefore unlikely to significantly contribute to the integrated diffuse γ-ray emission, although it might contaminate the significance map around Cygnus OB2.

Of the three supernova remnants observed within 2 kpc towards the Cygnus region, W63 and the Cygnus Loop are older remnants that, while perhaps efficient accelerators in the distant past, are too far away from the centre of emission to account for the integrated flux. γ-Cygni is a younger remnant (about 5 kyr, Lozinskaya et al. 2000) which is bright in radio, X-rays, and γ-rays. It expands in a low-density environment (about 0.1 cm⁻³) to a present size of about 15 pc with a shock velocity of nearly 1000 km s⁻¹ (Lozinskaya et al. 2000). This low-density environment could be part of the superbubble blown by Cygnus OB2 or the cavity from NGC 6910 (depending on its position along the line-of-sight). The γ-ray spectrum associated with γ-Cygni extends to tens of TeV (Fleischhack 2019) and is thought to originate from proton-proton interactions with a nearby cloud (MAGIC Collaboration 2023). This component is already subtracted in the data analysis of Fermi-LAT, HAWC, and LHAASO. It is, however, not excluded that ultra-high-energy γ-ray emitting particles were accelerated in the early expansion phase of the remnant, have already escaped the γ-Cygni region to fill a network of cavities in the Cygnus star-forming complex, and have eventually reached the molecular clouds on its edge.

Finally, the microquasar Cygnus X-3 is seen in coincidence with Cygnus OB2 and therefore lies in the centre of the diffuse γ-ray emission. Similar microquasars have been recently recognised as efficient particle accelerators up to 100s of TeV (LHAASO Collaboration 2024b). At an estimated distance of 9 kpc, Cygnus X-3 is unlikely to contribute to the bulk of the emission. Particles accelerated in the system would have to be confined in a physical radius of over 300 pc to reproduce the 2° extent of the diffuse emission. However, a subdominant contribution at ultra-high energies appears plausible, as we discuss in Sect. 6. Given these details, the best accelerator candidates in the Cygnus region are most likely linked to the energy input from massive stars, including the possibility of a supernova remnant in the low-density cavity around the association.

2.3 Particle acceleration

Winds from massive stars or star clusters are thought only able to accelerate PeV particles if they are exceptionally powerful and strongly magnetised. The physical requirements are detailed in Vieu et al. (2022b) and Härer et al. (2025). A key limitation in wind acceleration models comes from the necessity to overcome adiabatic losses when accelerated particles propagate in the radially expanding region upstream of the wind termination shock. We note that the argument applies both in the case of an isolated massive star or a compact star cluster powering a collective outflow. Disregarding the downstream residence time, which might be shortened by field amplification, shock induced turbulence or wind-wind collisions, the acceleration rate at a strong shock reads ${\dot{p}}_{\text{acc}}\approx V_{\text{W}}^{2}p/\left(4D\right)$, where V_W is the wind speed and D is the upstream diffusion coefficient (Lagage & Cesarsky 1983; Drury 1983). This must always be larger than the adiabatic loss rate in the upstream, which reads ṗ_loss ≥ 2V_Wp/(3R_WTS), where R_WTS is the radius of the wind termination shock. The particle energy is therefore bounded as D(E) < 3R_WTS V_w/8. The most efficient acceleration is reached when the diffusion is in the Bohm limit in a magnetic field of strength B: D(E) = Ec/(ξqB) with q the charge of the particle, c the speed of light and ξ ≤ 3. The maximum energy is then limited as E_max < qBR_WTSV_w/c, which is recognised as the Hillas limit in the upstream region (Hillas 1984).

Independently of the assumption made for the generation or amplification of the magnetic field, the magnetic field strength is limited by the requirement that the wind must be super-Alfvénic for acceleration to proceed. This sets a loose limit on B as B = $\sqrt{4\pi \rho }{V}_{\text{w}}/M_{\text{A}}$, where ρ is the gas density and efficient particle acceleration requires M_A substantially exceeds unity. Conservation of mass in a spherical wind satisfies $\rho =\dot{M}/\left(4\pi R_{\text{WTS}}^2{V}_{\text{W}}\right)$, which leads to $E_{\max }\ll q\sqrt{V_{\text{W}}{L}_{\text{W}}}/c$, where $L_{\text{W}}=0.5\dot{M}{V}_{\text{W}}^2$ is the mechanical power and Ṁ the mass-loss rate of the central source. Remarkably, this upper limit only depends on the properties of the star and is in particular independent of the geometry of the shock. The limit barely reaches 1 PeV for the most powerful stars in the heart of Cygnus OB2 (WR144, WR146, Schulte 7). In practice, this limit is most likely never reached as it requires unreasonable assumptions for the amplification of magnetic fields.

The background magnetic field can be more accurately estimated from a standard model of stellar magnetic fields, that is, a dipolar structure below the Alfvén point and a primarily toroidal structure at large distances (e.g. Usov & Melrose 1992). The most powerful O stars in Cygnus OB2 (e.g. Schulte 7, Schulte 11) have wind terminal velocities around 2500–3000 km s⁻¹, mass-loss rates around 5 × 10⁻⁶ M_⊙ yr⁻¹, stellar radii around 15–20 R_⊙ and projected rotational velocities around 100–150 kms⁻¹ (Berlanas et al. 2020). Assuming a strong surface field ~ 1 kG, the limit on the maximum energy is at most 300 TeV. To reach energies above 1 PeV requires a surface field higher than 10 kG which is arguably not realistic, in particular since high magnetic fields are correlated with low rotational velocities (e.g. Grunhut et al. 2017). For Wolf-Rayet stars, the mass-loss rate is typically enhanced by a factor ten but the stellar radius is only a few solar radii (Sander et al. 2019). Single Wolf-Rayet stars seem also to be slow rotators (Meynet & Maeder 2003). Assuming a stellar rotation velocity ~200 km s⁻¹, surface magnetic field ~1 kG and stellar radius ~5 R_⊙, the limit on the maximum energy is 200 TeV. In the case of the very powerful star WR144, which dominates the power output in Cygnus OB2 with a mechanical luminosity over 9 × 10³⁷ erg s⁻¹, but whose radius is only about 1 R_⊙ (Sander et al. 2019), the maximum energy falls short of the PeV regime even for extremely large surface fields or rotational velocities.

In addition to the background field, cosmic-ray-induced instabilities are believed to amplify magnetic field fluctuations upstream of strong shocks (e.g. Bell & Lucek 2001; Bell 2004). However, the efficiency of these instabilities is limited in low-density environments. For instance, the non-resonant streaming instability at saturation is argued to generate a turbulent field $\delta B_{\text{sat}}^2\approx 4\pi {\xi }_{\text{cr}}\rho V_{\text{W}}^2/c$ which, for an acceleration efficiency ξ_cr = 0.1, a density ρ = 0.01m_p cm⁻³ and a wind velocity V_w = 3000 km s⁻¹, gives δB_sat = 4 µG which falls short of what would be necessary to accelerate PeV particles. These arguments dismiss individual stellar winds as well as compact clusters (NGC 6910, NGC 6913) as sources of the ultra-high-energy photons observed from Cygnus. We stress that these fundamental arguments apply independently of the boundary conditions downstream of the stellar wind and therefore cannot be circumvented by invoking downstream confinement effects (e.g. via wind-wind collisions, Bykov et al. 2013). Based on the argument in the previous paragraph, only Cygnus OB2, with a mechanical power well above 10³⁸ erg s⁻¹, could conceivably accelerate PeV particles, though this requires the association as a whole is able to generate a collective large-scale wind termination shock. However, Cygnus OB2 is not a compact cluster: its stellar distribution extends over a diameter of 30 pc with no clear sign of mass segregation. The mechanical power is therefore diluted in a large region and direct wind-wind interactions are not efficient enough to generate a large-scale cluster wind termination shock, as demonstrated by Vieu et al. (2024b). As such, Cygnus OB2 should be seen as a loose collection of individual stellar winds embedded in a low-density superbubble powered by thermalisation of these stellar winds. The only possible sites of particle acceleration are the individual wind termination shocks around the most powerful O and Wolf–Rayet stars, and these cannot accelerate PeV particles, at least not within our current understanding of microphysical processes near non-relativistic shocks.

In summary, the most likely sources of high-energy (100 MeV < E < 100 GeV) to very-high-energy (100 GeV < E < 100 TeV) diffuse γ-ray emission in the Cygnus region are the stellar winds of the most powerful O and Wolf-Rayet stars, in particular stars lying at the heart of Cygnus OB2, which can provide a stationary supply of non-thermal particles over millions of years. On the other hand, at their current stages of evolution, none of the known sources can reasonably account for the ultra-high-energy (E > 100 TeV) emission observed by LHAASO. The ultra-high-energy accelerator must therefore be either hidden or transient.

3 Particle acceleration in a past supernova remnant

3.1 The case for a recent supernova in Cygnus OB2

Cygnus OB2 is likely not coeval. Age estimates for its sub-populations range from 3 to 7 Myr (see e.g. Wright et al. 2015, and references therein). The presence of three Wolf-Rayet stars, 24 supergiants and one pulsar is evidence that late phases of stellar evolution have already been reached and it is likely that supernovae have exploded in the past. To make this argument more quantitative, we show an example calculation of the expected supernova rate using the Hoki package (Stevance et al. 2020), which is an interface for models computed with the Binary Populations and Spectral Synthesis (BPASS) code (Stanway & Eldridge 2018). The BPASS models calculate stellar evolution tracks for a wide range of initial parameters, include binary evolution effects, and allow for the calculation of transient rates such as supernovae. For the calculation of supernova event rates, we assumed Solar metallicity and a single stellar population with a total initial mass of 1.65 × 10⁴ M_⊙ (Wright et al. 2015). We assumed an initial mass function index α = −2.0. This is steeper than that suggested for the massive stellar population by Wright et al. (2015), but given that the known OB population has increased by ~25% since that estimate (Berlanas et al. 2018) this value is reasonable. The calculated rates are shown in Fig. 2. For an age of 3–5 Myr, type Ic supernovae are by far the most likely type of supernova, and the rate is consistent with several supernovae per Myr. Type Ic supernovae are likely to exceed the canonical supernova power of 10⁵¹ erg by a factor of a few (Fryer et al. 2018). The assumption that a powerful supernova has occurred in the last few tens of kiloyears in Cygnus OB2 is thus well motivated.

Despite indirect hints of past explosions (e.g. the pulsar PSR J2032+4127 located in Cygnus OB2, Lyne et al. 2015, or runaway stars, Comerón & Pasquali 2007) no supernova remnant signature has been confirmed in the close vicinity of Cygnus OB2, as discussed in Sect. 2. This is perhaps unsurprising. A supernova remnant expanding in the low-density environment carved by the stellar winds around the core of Cygnus OB2 would quickly reach a radius of several tens of parsec and fade below the detection capabilities of radio, optical, and X-ray observatories after a few ten-thousand years (Chu 1997). Moreover, the shock is never expected to enter the radiative phase but rather slow down to subsonic velocities after a few hundred thousand years (Parizot et al. 2004). For these reasons, old supernova remnants expanding in hot bubbles are only expected to leave extended, diffuse, and faint signatures. These would be difficult to detect, in particular in a complex, inhomogeneous environment with, as in the case of Cygnus, a superposition of extended multi-wavelength signals along the line-of-sight.

Assuming a supernova exploded a few tens of thousands of years ago in the Cygnus OB2 association, the maximum energy of particles accelerated at the forward shock is estimated as (Vieu et al. 2022a, b) $$E_{\max }\approx \left(\frac{B}{5\mu \text{G}}\right)\left(\frac{V_{\text{SN}}}{{10}^4\text{km}{\text{s}}^{-1}}\right){\left(\frac{M_{\text{ej}}}{1M_{\odot }}\right)}^{\frac{1}{3}}{\left(\frac{n_{\text{SB}}}{0.01{\text{cm}}^{-3}}\right)}^{-\frac{1}{3}}\text{PeV}.$$(1)

This estimate assumes that the remnant expands in a low-density cavity powered by a loose OB association (in contrast to the case of a compact star cluster). For loose associations, the maximum energy reaches at most a few PeV for powerful supernova remnants (see also Vieu & Reville 2023). Yet, we show in the following that this can be sufficient to account for the steep component of the ultra-high-energy γ-ray spectrum.

Fig. 2 Event rates for different types of core-collapse supernovae for a cluster with a total initial mass of 1.65 × 104 M⊙, as appropriate for Cygnus OB2, computed using Hoki.

3.2 Hydrodynamic simulation of a supernova remnant in Cygnus OB2

To estimate quantitatively the maximum energy of particles accelerated by a supernova remnant shock in Cygnus OB2, we require a prescription for the dynamical evolution of the remnant, accounting for the inhomogeneous superbubble properties ahead of the shock. To this end, we explode a supernova in the previously published hydrodynamic simulation of a Cygnus OB2-like association (Vieu et al. 2024b). The simulation is resumed at 1.3 Myr, at which point the high-mass supernova progenitor is turned into a Wolf-Rayet star with Ṁ = 10⁻⁵ M_⊙ yr⁻¹ and V_w = 2500 km s⁻¹. The supernova is then initialised at 1.6 Myr with a linear radial velocity profile and a flat density profile with a power-law cut-off ρ ∝ r⁻⁹, as described in Whalen & Norman (2008) and Badmaev et al. (2024). The explosion energy is set to 5 × 10⁵¹ erg and the core velocity to 16 000 km s⁻¹, values consistent with a powerful type Ic supernova (Fryer et al. 2018; Barbarino et al. 2021). We neglect thermal conduction and cooling, which are not expected to affect the early phases of expansion.

The forward shock of the simulated supernova remnant expands beyond the inner core of Cygnus OB2, as shown in Fig. 3, and is progressively decelerated. Interaction with the pre-supernova medium, in particular with powerful Wolf–Rayet winds, slows down the expansion further along some directions, producing asymmetries, though the shock expansion is not blocked by the stellar winds. The remnant expands well-beyond the inner core after several thousand years. As the density of the background medium is low (about 0.01 cm⁻³), the free-expansion phase proceeds for a long time, up to about 2 kyr, providing favourable conditions for particle acceleration.

3.3 Maximum energy of accelerated particles

The simulation was analysed to track the properties of the forward shock as a function of time. The top row of Fig. 4 illustrates how the properties are distributed across the shock surface at 2 kyr. At this point, the shock has reached a radius of about 15pc and the velocity has decreased to about 3000 km s⁻¹, such that diffusive shock acceleration is no longer effective at the required energies. In order to infer the maximum energy of particles accelerated up to this time, we integrate the acceleration rate in time. Since the upstream magnetic field is not calculated in the hydrodynamic simulation, we assume that 10% of the kinetic energy in the superbubble goes into turbulent magnetic energy, defining $B_{\text{eq}}=u\sqrt{4\pi {\eta }_{\text{B}}\rho }$ with η_B = 0.1 and u the pre-supernova velocity field in the superbubble. Cosmic-ray induced instabilities are further expected to enhance the turbulent component of the background magnetic field and generate a flat turbulence spectrum, although these instabilities grow slowly in a low-density medium. In particular, the non-resonant streaming instability operates only if the cosmic-ray current is strong enough to overcome the magnetic field tension. This condition reads (Bell 2004) $${\xi }_{\text{cr}}{u}^3>{\upsilon }_{\text{A}}^{2}c,$$(2)

where u is the shock velocity, υ_A is the upstream background Alfvén speed, and ξ_cr is the acceleration efficiency that we set to 10%. The background field, B_eq, inferred from the simulation is typically about 1–10 µG. These values are consistent with those we inferred from magnetohydrodynamic simulations in Härer et al. (2025). For such B_eq, the condition 2 is satisfied most of the time upstream of the supernova remnant shock. Provided the streaming instabilities operate, the saturated magnetic field can be estimated as (e.g. Bell & Lucek 2001; Bell 2004; Matthews et al. 2017) $$\delta B_{\text{Sat}}^2\approx \{\begin{array}{cc}2{\xi }_{\text{cr}}u\sqrt{4\pi \rho }{B}_{\text{eq}}& \text{RSI}\\ 4\pi {\xi }_{\text{cr}}\rho u^3/c& \text{NRSI}\end{array},$$(3)

where RSI denotes the resonant and NRSI the non-resonant streaming instability. In all cases, we assumed that diffusive shock acceleration operates at maximum efficiency (Bohm diffusion), in which case the acceleration rate is written as (Lagage & Cesarsky 1983) $$\frac{\text{d}p}{\text{d}t}\approx \frac{u^{2}p}{20\kappa },\kappa =pc/\left(3qB\right),$$(4)

where B is the turbulent magnetic field. We repeat the calculation assuming first that streaming instabilities do not operate (B = B_eq), then that RSI operates (B = δB_sat,RSI) and finally that NRSI operates (B = δB_sat,NRSI if the left hand-side of Eq. (2) is ten times the right hand-side, B = B_eq otherwise). We integrate the acceleration rate as function of time along a large number of independent lineouts, reading the local shock velocity, density and inferred magnetic field at every timestep. An implicit assumption of this procedure is that particles do not propagate along azimuthal directions from one lineout to another. This is a reasonable assumption if the magnetic field is parallel to the shock normal. At every time-step we check that the Hillas confinement limit, E < qBRu/c, where R is the distance travelled by the forward shock, is not violated, otherwise integration is not performed. Efficient particle acceleration stops after about 2 kyr when the supernova remnant enters the Sedov-Taylor phase and the highest energy particles are no longer effectively confined.

The second row of Fig. 4 shows the maximum energy obtained for protons at the end of the phase of efficient particle acceleration. Without streaming instabilities, the supernova remnant shock is overall able to accelerate efficiently up to about 500 TeV. With streaming instabilities, a maximum proton energy of a few PeV is reached. We note that we naively expect the maximum energy to scale as the square root of the explosion energy such that a nominal supernova (10⁵¹ erg) is expected to accelerate protons up to 0.5–1 PeV with streaming instabilities.

In the approach adopted above, we chose optimal parameters for confinement and acceleration of the highest energy particles, as required to reach the Hillas limit. This is usually motivated by the fact that cosmic rays self-excite the turbulence. Independently of this mechanism, studies by Bell et al. (2025) (see also Lazarian & Xu 2021) find that particle mirroring in large-scale fields - those larger than the gyroscale of the accelerating particles – can exceed quasi-linear theory estimates of particle confinement near shocks, even with moderate (sub-Bohm) scattering rates. Large-scale field structures naturally result from the turbulent interactions of magnetised stellar winds within superbubbles (see e.g. Härer et al. 2025). These findings provide confidence that the Hillas limit could be achieved without resorting to the Bohm scattering limit.

Fig. 3 Density slices from the simulation of a supernova remnant expanding in Cygnus OB2. (For details, see Sect. 3.2.)

Fig. 4 Top: supernova remnant properties across the forward shock surface at 2 kyr, extracted from the hydrodynamical simulation shown in Fig. 3. Bottom: integrated maximum energy at 2 kyr for the resonant and non-resonant streaming instabilities (RSI and NRSI, respectively) and without streaming instabilities. Beq was computed assuming that 10% of the kinetic energy in the superbubble goes into a turbulent magnetic field.

4 Molecular gas

In this section, we reconstruct the 3D molecular gas distribution in the Cygnus region from previous studies. This is subsequently used to compute the γ-ray emission (see Sects. 5 and 6).

The CO intensity maps reveal a huge molecular complex from about 76° to about 85° in longitude and from about −3° to about +4° in latitude. Within this complex, five regions are especially prominent: the cloud in the far-north that separates the North America and the Pelican nebulae, the northern complex between 80–83° (Cygnus-X North), the southern complex between 77–80° (Cygnus-X South), the north-west clouds at high latitudes, and the lane extending in the far-south, including the S 106 region. It has been established that this molecular complex is not a coherent structure but a superposition of clouds along the line-of-sight, including the well-known Cygnus Rift in the foreground (e.g. Reipurth & Schneider 2008).

Schneider et al. (2006) established the first comprehensive census of the molecular clouds in Cygnus-X, using ¹²CO and ¹³CO emission lines. Because the Cygnus region is seen tan-gentially to the galactic rotation, kinematic velocities do not translate into reliable distance measurements. Restricting velocities to the range (−20, 30) km s⁻¹ allows to exclude the Perseus Arm and beyond; however, the 3D distribution of the gas in the Local Arm (within 3 kpc) cannot be directly modelled. Nevertheless, scanning the intensity map along the velocity axis can give some information about superimposed components, at least at the qualitative level. The underlying assumption is that coherent structures correlate in position-position-velocity space.

After carefully scanning all prominent structures in the CO map, Schneider et al. (2006) extracted four seemingly coherent groups, in addition to the foreground Rift. Group I is a filamentary cloud north of Cygnus which hosts several dense clumps forming protoclusters seen as radio sources, for example DR21, DR20, DR22, and DR23 around l = 81° in the velocity range (−7, 1) km s⁻¹. In the following we refer to this group as the ‘North Filament’. Group II partly overlaps with the filament but is significantly decorrelated in velocity space with a range (6, 20) km s⁻¹, and was placed in the foreground by early absorption studies (Dickel & Wendker 1978). Parallax measurements of masers (Rygl et al. 2012) place DR2O and DR2l (North Filament) at about 1.5 kpc, and group II, which includes in particular W75N and DR17, in the foreground at about 1.3 kpc (see also Gottschalk et al. 2012). Group III gathers compact clumps and globules around Cygnus OB2, notably DR18 and the southern arm of DR20. The star-forming region AFGL 2591 was originally included in this group; however it has since been relegated to the background (about 3.3 kpc) by a maser measurement (Rygl et al. 2012). Finally, most of the southern clouds were gathered in a single Group IV, and the remaining unassociated features were attributed to the foreground Cygnus Rift (Group V) between 0.6–0.8 kpc. Schneider et al. (2007) completed the census in the southernmost part up to the S 106 HII region and showed that several substructures in this area are shaped by stellar feedback in a star-forming region at 1.6–1.7 kpc, which hosts notably the compact star clusters NGC 6913, Berkeley 86 and Berkeley 87.

The groups identified by Schneider et al. (2006) can serve as a basis to identify substructures in Cygnus; however, the relative masses of the groups are highly uncertain, and accurate distances cannot be directly inferred. In addition, the assumption that coherent structures in three-dimensional space correlate in position-position-velocity space is doubtful. Turbulence can dramatically spread the velocity component across a given molecular cloud, especially near active star-forming-regions, which leads to over-decomposition of coherent clouds. On the other hand, overlapping clouds can have similar velocities by coincidence, which leads to confusion between the various components along the line-of-sight. The advent of the Gaia era, enabling precise parallax measurements, has shed new light onto the structure of the molecular clouds in Cygnus. Zhang et al. (2024) employed a clustering algorithm to identify CO clouds and measured their distance via extinction jumps along the line-of-sight. This successfully disentangles foreground components, showing that the Rift actually extends up to 1.1 kpc and is much more massive than thought by Schneider et al. (2006). In particular, it becomes clear that the Group IV in Schneider et al. (2006) is largely confused with the foreground Rift. In addition, Zhang et al. (2024) measure a massive cloud in the south with a distance estimate of 1.3 kpc, including DR13 and its prominent filament (red triangular cloud in Fig. 5). This brings a southern counterpart to Schneider Group II. We propose to name these two groups in the 1.3 kpc layer ‘Cygnus 1.3 North’ and ‘Cygnus 1.3 South’ for clarity.

Fig. 5 Main molecular clouds in the Cygnus region. Red: Cygnus 1.3 North at l > 80° (Schneider Group II, 1.3 kpc) and Cygnus 1.3 South at l < 80° (part of Schneider Group IV, 1.3 kpc). Orange: Cygnus North Filament at l > 80° (Schneider Group I, 1.5 kpc) and Cygnus SFR South at l < 80° (part of Schneider Group IV, unknown distance >1.4 kpc). Pink: Cygnus South Filament (part of Schneider Group IV, unknown distance > 1.4 kpc). Green: Cygnus Central Clumps (Schneider Group III, 1.6–1.8 kpc). Dark blue at high latitude: Cygnus West, unknown distance (>1.4 kpc). Turquoise: S106 regions (described in Schneider et al. 2007, 1.6–1.8 kpc). The grey background shows foreground clouds: the North-America and Pelican nebulae on the left, the North-West clouds, and the Rift lane extending to the right in the direction to the Galactic Centre.

Cygnus-X North

Unfortunately, the extinction becomes too high beyond ~1.4 kpc, preventing Zhang et al. (2024) to probe more distant molecular clouds, and in particular the clouds in the vicinity of the star-forming region at ~ 1.6–1.8 kpc (including Cygnus OB2). In Cygnus North, Zhang et al. (2024) were able to associate a distance to about 50% of the total mass, which represents about 115 thousand solar masses in the foreground of the North Filament. On the other hand, Schneider et al. (2006) associated a mass of 167 thousand and 55 thousand, respectively for the North Filament (Group I) and Cygnus 1.3 North (Group II) which, corrected for the updated distance measurements, adds up to about 160 thousand solar masses. Within mass inference uncertainties, it is therefore reasonable to assume that the census is complete in the north, with Cygnus 1.3 North and the North Filament lying behind the Rift. Using the ¹²CO data from Zhang et al. (2024), we extract the shape of the North Filament and Cygnus 1.3 North in the velocity ranges (−7, 1) km s⁻¹ and (7, 15) km s⁻¹, respectively. These are shown in red and orange in Fig. 5. We note that the most recent distance estimates place Cygnus OB2 at about 1.65 kpc, which is behind the North Filament (1.5 kpc). The North Filament might actually lie at the edge of the putative Cygnus OB2 superbubble, which could explain its peculiar shape (Li et al. 2023).

Cygnus-X South

The situation is less clear in the south. Zhang et al. (2024) were able to associate a distance to approximately 150 thousand solar masses, corresponding to only 36% of the total mass. The missing mass most likely lies beyond 1.4 kpc, possibly embedded in the interfaces of photodissociation regions at the distance where several active star-forming regions are located, beyond the reach of the extinction jump method (>1.3 kpc). We suggest to call this hidden group ‘Cygnus SFR South’ (SFR standing for ‘star-forming region’). In order to extract this group, we identify regions of the position-position-velocity space which do not correspond to the clouds catalogued by Zhang et al. (2024). These regions appear in the velocity ranges (−10, −5) km s⁻¹, (2.2, 4.5) km s⁻¹ and (4, 8.5)km s⁻¹ around 78–79° in longitude. The result defines the ‘Cygnus SFR South’ region shown in orange at l < 80° and −1° ≲ b ≲ 1° in Fig. 5. In the range (12, 20) km s⁻¹, a filamentary cloud appears, which we call the ‘South Filament’ (shown in pink in Figure 5). Finally, another region at high latitude, b ≈ 4° in the velocity range (−5, 3) km s⁻¹ is extracted and labelled ‘Cygnus West’. We assume in the following that all these southern clouds beyond the reach of the extinction jump method are located within the active star-forming region at about 1.7 kpc, loosely coincident with the photodissociation regions. This aligns with a plausible scenario where these clouds are at the intersections of multiple cavities that were carved in the past by the joint feedback of Cygnus OB2, NGC 6913, NGC 6910, and possibly other clusters and associations (see Sect. 1).

Central clumps

In addition to the massive molecular clouds in the north and south, a number of clumps, sometimes called ‘globules’ or ‘condensations’, are seen around Cygnus OB2, with hints of direct interaction with the radiative feedback from the association (e.g. Schneider et al. 2016). A number of these clumps are listed in Table 3 of Zhang et al. (2024), to which we further add the prominent regions DR 18 and DR2O in the velocity range (− 11, 1) km s⁻¹ around l = 81°. We group these clumps under the common label ‘Cygnus Central Clumps’ and assume a distance similar to that of Cygnus OB2 (1.6–1.8 kpc).

Other groups

The remaining foreground groups (North-West, 1.3 South) can be directly identified from the distance measurements given in Zhang et al. (2024). Finally, Cygnus Far-South can be extracted manually according to its well-defined shape described in Schneider et al. (2007). Table 1 summarises our census of the molecular cloud complexes in the Cygnus region. The listed masses are inferred by assuming a X_CO factor of 4 × 10²⁰ H/cm² (K km s⁻¹)⁻¹. For later convenience, we introduce shorthands A–K to refer to the clouds.

5 Modelling particle propagation and γ-ray emission

In the previous sections, we motivated that a powerful past supernova in Cygnus is a plausible source of very-high-energy particles. We estimated the maximum energy of particles accelerated at the forward shock and discussed the 3D distribution of target material in molecular clouds. Here, we outline our approach to model particle transport and γ-ray emission. We construct both hadronic and leptonic models and apply them to supernova and stellar-wind scenarios in Sects. 6 and 7.

Fig. 6 Loss timescales for electrons in the putative superbubble surrounding Cygnus OB2. The photon field from stars in Cygnus OB2 is evaluated at 20 pc. Diffuse dust and starlight are taken from the model by Popescu et al. (2017) and enhanced by a factor of two. Coulomb losses are shown for a hot, low-density, ionised medium (T = 106 K, nSB = 0.05 cm−3) following Schlickeiser 2002 (Sect. 4.5). All other timescales are computed with the GAMERA code (Hahn et al. 2022).

5.1 Losses, photon fields, and magnetic field

Knowledge of ambient photon fields and the magnetic field is required to compute leptonic losses and inverse Compton emission. We model galactic diffuse starlight and dust-scattered starlight following the axisymmetric model of Popescu et al. (2017), with an enhancement factor of two for both photon fields, to reflect increased activity within the spiral arm (cf. Breuhaus et al. 2021). The cosmic microwave background is included as a blackbody with a temperature of 2.7 K and energy density of 0.26 eV cm⁻³. We model the photon fields of individual massive stars in Cygnus OB2 with blackbody spectra, taking luminosities and temperatures from Wright et al. (2015). The total photon field is obtained by summing all spectra. The energy density at a distance R is obtained from the luminosity, L, as L/(4πR²c). At a distance of 20 pc, the energy density of the cluster photon field is 25 eV cm⁻³, which is comparable to the model value obtained by Schneider et al. (2016), who take into account the spatial distribution of stars. Figure 6 shows cooling time-scales for electrons, for a magnetic field of 5 µG and a density of 0.05 cm⁻³, which is representative of the superbubble interior. We assume that the magnetic field is constant over the superbubble interior. In the GeV range, the dominating photon field is the cluster photon field. At higher energies, losses are dominated by synchrotron cooling. Losses for nuclei are neglected as their timescales are large in the low-density superbubble environment.

Fig. 7 Neutral gas model. The model consists of two 3D-Gaussian components, the Cygnus Rift (left) and the Cygnus SFR (middle), with a spherical (radius 150 pc), low-density (0.05 cm−3) cavity around Cygnus OB2. The total mass is set to match the mass obtained from integrating the neutral hydrogen column density obtained by Astiasarain et al. (2023). For model parameters, see Table 2.

5.2 Target material

We consider proton-proton interactions in molecular clouds and neutral gas. The distances of clouds A–B and C–D in Table 1 are fixed to 1 pc and 1.3 kpc respectively. We note that, assuming a distance of 1.6 kpc for Cygnus OB2, the density of particles at clouds A–D is too low to materially contribute to the flux. These clouds however contribute to diffuse emission from galactic cosmic rays. Beyond 1.4 kpc, extinction introduces considerable uncertainty in cloud positions (see Sect. 4). Clouds E–K are initially assumed to be equidistant to Cygnus OB2 and close to the superbubble edge. This is a plausible scenario: density increases in the superbubble shell, which is likely to be thin and fragmented due to cooling (see Weaver et al. 1977). Such a scenario is also plausible considering the CO data (see Sect. 4).

We construct a neutral gas model based on measurements of the column density of atomic hydrogen (HI) by Astiasarain et al. 2023 (see their Fig. 2). We associate the neutral gas with known molecular cloud complexes: the Cygnus Rift at a distance of about 1 kpc and the Cygnus SFR at about 1.6–1.8 kpc. Each component is modelled as an ellipsoid with parameters listed in Table 2. Column density maps are shown in Fig. 7. Model parameters and total masses were chosen to account for the total column density. Thus the model assumes that all mass is contained in the two components mentioned above.

As discussed in Sect. 1, there are a number of observations hinting at the existence of a low-density ionised cavity around Cygnus OB2, which shapes the nearby molecular clouds (e.g. Li et al. 2023). The shape and full extent of this cavity cannot be directly inferred from observations. According to Weaver et al. (1977), assuming a continuous injection of energy and a homogeneous ambient medium, the radius of such a wind blown bubble is expected to be $$R_{\text{b}}=183{\left(\frac{{\xi }_{\text{b}}{L}_{\text{W}}}{{10}^{38}\text{erg}{\text{s}}^{-1}}\right)}^{1/5}{\left(\frac{n_{\text{ISM}}}{1{\text{cm}}^{-3}}\right)}^{-1/5}{\left(\frac{t}{5\text{Myr}}\right)}^{3/5}\text{pc,}$$(5)

where L_w is the wind power, n_ISM the ambient density, t the age of the system, and an empirical correction factor (see Vieu et al. 2022a). Modelling a spherical superbubble is a necessary simplification given the difficulty to observationally determine its real boundaries. Based on the current power, past stellar activity (see Sect. 2.1) and observational hints, we set the superbubble radius to 150 pc from the centre of Cygnus OB2 (l = 80.2°, b = 0.8°). Inside the superbubble, we set the density to a constant value of 0.05 cm⁻³. This value lies in the range typical for super-bubbles (Weaver et al. 1977, see also Sect. 1). We note that the density obtained in the hydrodynamical simulation (see Sect. 3) is lower because thermal conduction is not included, and hence evaporation at the shell is omitted.

5.3 Particle injection and propagation

We assume that at the end of the efficient particle acceleration phase, particles are released into the low-density superbubble carved by the winds from massive stars in Cygnus OB2. In the resulting turbulent environment, the diffusion coefficient, D_SB, is expected to be smaller than its average interstellar value. We assume $D_{\text{SB}}={10}^{25}{p}_{\text{GeV}}^{1/2}{\text{cm}}^2{\text{s}}^{-1}$ which, in the framework of quasi-linear theory, is expected for a magnetic field around 10 µG and a coherence length of a couple of pc. Beyond the low-density superbubble, the diffusion coefficient is increased to its fiducial value in the interstellar medium, $D_{\text{ISM}}={10}^{28}{p}_{\text{GeV}}^{0.4}{\text{cm}}^2{\text{s}}^{-1}$.

We model large-scale advection in a mean r⁻² wind profile with a speed of 100 km s⁻¹ at 10 pc. We note that advection is only relevant at low energies, especially for our chosen initial velocity, which is appropriate for loose associations such as Cygnus OB2 (see Vieu et al. 2024b). For electrons, we take inverse Compton, synchrotron, and bremsstrahlung losses into account, as described in Sect. 5.1. The resulting inhomogeneous transport equation for the particle distribution function ƒ(r, p, t) reads $${\partial }_{t}f+u{\partial }_{r}f-\frac{1}{r^2}{\partial }_r\left(r^{2}D{\partial }_{r}f\right)-\frac{1}{p^2}{\partial }_p\left(p^2\frac{\text{d}p}{\text{d}t}f\right)=Q\left(r,p,t\right),$$(6)

where we adopt spherical symmetry. Adiabatic losses vanish for our adopted r⁻² radial velocity profile. Q(r, p, t) is a source term that can represent, for instance, impulsive injection from a supernova remnant or stationary injection from wind termination shocks. The transport equation is an example of a advection-diffusion-reaction equation. We solve it with a finite element discontinuous Galerkin method, adapting the method presented in Di Pietro & Ern (2012, Sect. 4.6) for time-dependent problems, using a Crank-Nicolson integrator and an upwind-flux for the advection term. For efficient computations, we use logarithmic p and r variables. We use zero Dirichlet boundary conditions everywhere, suppressing the diffusion term below 10 pc to avoid propagation to the centre. We made sure that the solution was not impacted by the choice of boundary conditions, especially checking that the outer r boundary was far enough away for the solution to sufficiently decay.

The solver was implemented in C++ using the open-source finite-element library deal.II (Africa et al. 2024; Arndt et al. 2021). The code was validated against a number of tests for which analytic solutions are known, including homogeneous synchrotron losses in a box, advection-diffusion without losses, and an impulsive Gaussian profile with homogeneous diffusion coefficient. The code provides the solution for ƒ(r, p, t) which can be post-analysed to determine the non-thermal emission.

5.4 Gamma-ray emission

We computed non-thermal emission using the GAMERA code (Hahn et al. 2022)¹. In Sect. 4, we provide the details of our model for the 3D distribution of target material, which includes both the molecular clouds and neutral medium. Convolving this with ƒ(r, p, t) yields both the full spectra integrated over the region and synthetic γ-ray emission maps. For the leptonic model, we evaluate the UV photon field from Cygnus OB2 as a function of radius and keep all other photon fields and the magnetic field strength constant.

We fitted the γ-ray morphology by adjusting the positions of clouds E–J. Parameters of the neutral gas model were kept fixed, including the superbubble interior density of 0 05 cm⁻³. Emission from proton-proton interactions in molecular clouds was evaluated at the barycentre of each individual cloud. We evaluated the neutral gas distribution on a 3D grid with cell size 16.8′ by 16.8′ in l and b and 24 pc in distance.

The diffuse galactic cosmic-ray background is modelled following Schwefer et al. (2023). We evaluate their model for atomic and molecular hydrogen and increase the resulting flux by a factor of 1.7 to account for contributions from heavier nuclei (Kafexhiu et al. 2014). We assume a distance of 1.6 kpc to Cygnus OB2 to convert luminosity to flux.

6 Hadronic emission from a supernova remnant

6.1 Set-up and results

In Sect. 3, we concluded that, despite the absence of a detectable remnant, a supernova likely occurred in Cygnus OB2 in the past 10–100 kyr and inferred the maximum energy of accelerated particles from a hydrodynamic simulation. Based on this estimate (see Fig. 4), we assume a maximum energy of 2 PeV for the following calculation. To model the γ-ray spectrum of protons accelerated at the forward shock of the remnant, we follow the approach outlined in Sect. 5, assuming an exponential cut-off power-law injection, Q ∝ p^−s exp(−p/p_max), where p_max = 2 PeVc⁻¹, the injection momentum is set to 10 MeVc⁻¹ and s = 4. Particles are impulsively injected over a sphere of radius 20 pc around the centre of Cygnus OB2.

Figure 8 shows the region-integrated spectra for a supernova remnant age of 50 kyr and Fig. 9 the particle flux map, smoothed using a Gaussian kernel of width σ = 0.3° to compare to the data. To reproduce the morphology, we place clouds F–K at a distance of 137 pc behind Cygnus OB2 and the North Filament (cloud E) 140 pc in front of Cygnus OB2. For this distribution of target material, the model requires a total energy of 3 × 10⁵⁰ erg in non-thermal protons above 0.01 GeV. At 100 TeV–PeV, the spectrum is steep with slight curvature. In the TeV range, the bulk of the emission comes from interactions with gas inside the superbubble. At higher energies, contributions from both molecular clouds and neutral gas become comparable. The molecular cloud component shows a low-energy cut-off, which stems from the fact that no clouds are located close to the particle injection site and that the diffusion length of protons at these energies is significantly smaller than the distance to the nearest cloud.

Overall, the observed morphology with a bright spot in the south and an additional weaker bright spot in the north is well reproduced (see Fig. 9). In addition, the model reproduces the energy dependent morphology between the 25–100 TeV and >100 TeV LHAASO bands. The morphology of the emission in the upper and middle LHAASO bands is dominated by emission from molecular clouds. The bulk of the emission comes from the North Filament and Central Clumps in the north. The south is dominated by the SFR South, with some contribution from the Cygnus Far South and Cygnus West clouds. The neutral gas component is required to reproduce the extent and steepness of the radial profile of the emission and the diffuse halo surrounding the clouds (see Appendix B). The energy dependence in the morphology stems the proton distribution, which, as shown in Fig. 10, is more peaked at small radii for 100 TeV particles than it is for 1 PeV particles. Consequentially, the brightness of the neutral gas emission varies with energy. The brightness change is less in molecular cloud emission since the clouds are located close to the superbubble edge.

6.2 Supernova remnant age, energy, and target material

The energy dependence of the spatial profile of the proton distribution function is greater for younger supernova remnants. Ages >80 kyr produce minimal energy dependence in the γ-ray morphology between the 25–100 TeV and > 100 TeV bands (see Fig. B.2). We note that the parameters t_SN and D_SB are degenerate as the proton distribution scales with $\sqrt{D_{\text{SB}}{t}_{\text{SN}}}$.

A different assumption on D_SB will therefore affect the age estimate.

The required energy corresponds to a powerful supernova with 3 × 10⁵¹ erg, if a proton acceleration efficiency of 10% is assumed. A powerful supernova is consistent with the assumption made in Sect. 3 to motivate a cut-off energy of a few PeV. The energy requirement can be reduced by assuming a superbubble density >0.05 cm⁻³, a younger supernova remnant, or, in principle, a higher cut-off energy. For younger (more recent) supernova remnants, particle escape at high energies is reduced, which leads to higher normalisations, harder spectra, and a cut-off at higher energy. A lower bound for the supernova age is probably 10–20 kyr, given that no remnant is detected in Cygnus OB2 (see Sect. 3.1).

In the scenario of higher superbubble density, molecular clouds would have to be slightly further embedded in the super-bubble to maintain the relative normalisation between the neutral gas and molecular cloud component, which is constrained by the γ-ray morphology (see Fig. B.1). For typical superbubble densities, n_SB ≲ 0.1 cm⁻³, the model predicts the clouds to be within ~ 10–20 pc of the edge. We note that the normalisation of the proton distribution changes rapidly across the superbubble edge (see Fig. 10). This is a consequence of the inhomogeneous two-zone diffusion model. A change in cloud position within a few parsecs can therefore produce measurable differences in the γ-ray maps. However, we highlight that in fitting the cloud positions, we only deviated from the assumption of equiradial clouds in the case of the North Filament, which is located at a radius of 140 pc instead of 137 pc as is the case for the remaining clouds in the vicinity of Cygnus OB2. In other words, except for the North Filament, the γ-ray brightness of clouds in our model is well correlated with the cloud mass. We remark in addition that cloud masses are subject to uncertainty (see Sect. 4). As highlighted in Sect. 5, approximately equiradial clouds close to the superbubble edge are expected. The fact that this assumption holds well within our model supports the validity of the supernova scenario.

Fig. 8 Preferred model for the region-integrated γ-ray spectrum consisting of hadronic emission from a 50 kyr-old supernova remnant along with inverse Compton emission from electrons accelerated at stellar-wind termination shocks in Cygnus OB2. The top panel shows both components in red and yellow as well as emission from diffuse galactic cosmic rays in blue. Diffuse galactic cosmic rays are only shown above 300 GeV since they are included in the Fermi-LAT background model. The middle panel shows hadronic model components. The lower panel shows contributions from inverse Compton scattering on individual photon fields. The data are from ARGO-YBJ Collaboration (2014), Abeysekara et al. (2021), Abdollahi et al. (2020), and LHAASO Collaboration (2024a).

6.3 Origin of the ultra-high-energy γ-rays

Substantial interest has been taken in the origin of the ultra-high-energy γ-rays from the Cygnus region. In the following, we explore how well the supernova remnant model describes the data up to the highest energies. Figure 11 highlights relevant uncertainties in the modelling: age of the supernova remnant and cut-off energy, the normalisation of the galactic cosmic ray background, and the contribution of Cygnus X-3 to the flux above 1 PeV. Our default assumption of E_max = 2 PeV and t_SN = 50 kyr gives a model within the σ uncertainty of the two highest-energy data points. Models with a lower cut-off or older supernova remnants miss these points (see Fig. 11, left panel). Our choice of E_max = 2 PeV is motivated for a powerful supernova (cf. Fig. 4). For a nominal supernova, a lower cut-off energy might be more appropriate (about 0.5–1 PeV for an explosion energy of 10⁵¹ erg assuming E_max is roughly proportional to the square root of the explosion energy). However, as argued in Sect. 3.1, the young age of Cygnus OB2 favours Type Ic super-novae whose energy likely exceeds 10⁵¹ erg (Fryer et al. 2018).

The background from diffuse galactic cosmic rays and the possibility of PeV particle acceleration in Cygnus X-3 introduce further uncertainty into the model. The uncertainty in the galactic cosmic-ray component is about a factor of two (Schwefer et al. 2023). With an increase in the background within this range, the constraints on supernova remnant age and cut-off energy can be relaxed slightly. However, models with high galactic cosmic-ray background tend to over-predicted the flux at TeV energies if the Fermi-LAT data is fit with inverse Compton emission from stellar winds in Cygnus OB2 (see Sect. 7) and require enhanced synchrotron cooling, B ≳ 10 µG.

Microquasars have recently received increased attention as potential accelerators of PeV particles. While at a larger distance, Cygnus X-3 is in spatial coincidence with Cygnus OB2 and two out of the eight photons reported in LHAASO Collaboration (2024a) above 1 PeV are spatially correlated with its location. In the right panel of Fig. 11, we naively scale the flux points above 1 PeV by a factor 2/8, leaving the error bars unchanged. This is meant only to illustrate how attributing some of the detected photons to Cygnus X-3 might affect the flux expected from the supernova remnant and should not be seen as proper modelling. A cut-off of 2 PeV provides the best description of the scaled data points (unscaled: 3 PeV). We note, however, that the uncertainties are large and that the flux from diffuse galactic cosmic rays is almost within 1σ.

In conclusion, considering that an association of at least some of the ultra-high-energy photons with Cygnus X-3 is likely and that the significance at high energies is low, the supernova remnant scenario is conceivable even for a lower cut-off energy. With standard assumptions (no contribution from Cygnus X-3 and standard galactic cosmic-ray background) a cut-off energy of 2 PeV is required to remain within the 1σ uncertainties of the data above 1 PeV, if the injection spectrum is a p⁻⁴ power-law with exponential cut-off.

Fig. 9 Synthetic γ-ray emission maps (top row) and significance maps from LHAASO Collaboration (2024a) (bottom row). Both the model and data have been smoothed using a Gaussian kernel of σ = 0.3°. The top of the colour scale is set to the maximum and the bottom to 10% of the maximum. The white cross marks the position of the stellar association Cygnus OB2. The figure demonstrates that the model is in good agreement with the γ-ray morphology and its energy dependence.

Fig. 10 Proton distribution function as a function of radius at 100 TeV (blue-dashed line) and 1 PeV (red-solid line). The dash-dotted vertical line marks the superbubble radius, which separates the low-diffusion and interstellar-medium diffusion zones. The solid vertical lines indicate the positions of molecular clouds in the model. The relative normalisation between the two displayed proton distributions changes strongly at small radii, which results in energy dependent morphology (see Fig. 9).

Fig. 11 Dependence of the model spectrum on cut-off energy and age of the supernova remnant (left panel) and uncertainty in the flux of galactic cosmic rays (middle panel). The right-hand panel shows the spectrum for the case that two of the eight photons with energy > 1 PeV can be attributed to the microquasar Cygnus X-3. In addition to a change in spectral shape, the energy required to match the normalisation of the data increases for lower cut-off energies and older supernova remnants. The displayed models require 4.5, 6, and 2.1 × 1050 erg in non-thermal protons for an age of 80 kyr and a cut-off energy of 1 PeV and 3 PeV, respectively.

6.4 Sensitivity to other parameters

We set $D_{\text{ISM}}={10}^{28}{p}_{\text{GeV}}^{0.4}{\text{cm}}^2{s}^{-1}$ and the superbubble radius to 150 pc. A change in these parameters within realistic bounds only affects the model marginally. Key to the behaviour of the model is that D_ISM ⪢ D_SB, which remains true even if D_ISM is reduced, for example, by a factor ten. Considering the discussion in Sect. 5.2, the superbubble radius might vary by as much as 20–30%, depending on properties of the environment and star cluster. A change of this magnitude does not affect the shape of the particle distribution significantly.

Finally, in our model we fix the spectral index at injection to the nominal value s = 4 obtained from test-particle diffusive shock acceleration theory (e.g. Drury 1983). This parameter is degenerate with other model parameters, in particular the supernova age, the positions of molecular clouds relative to the superbubble edge, and the diffusion coefficient in the superbubble. If other model parameters are kept fixed, steepening, be it from injection or propagation effects, reduces the fit quality at ultra-high energies. As discussed in Sect. 6.3, this can be compensated for by contributions from Cygnus X-3 or diffuse galactic cosmic rays.

6.5 Leptonic component

Electrons accelerated by the supernova remnant shock will produce inverse-Compton emission. This results in a component with a comparatively sharp peak at about 1 TeV, which cuts off at higher energies due to synchrotron losses, as shown in pink in Fig. 12 for an electron acceleration efficiency of 1% and assuming the same photon fields as described in Sect. 5.1. Leptonic emission from the supernova is subdominant compared to the much broader underlying inverse Compton component from stellar winds, such that a discernable observational signature in γ-rays is not expected.

7 Leptonic emission from Cygnus OB2

7.1 Results and discussion

Cygnus OB2 harbours several early-type and Wolf-Rayet stars, whose stellar-wind termination shocks are expected to accelerate particles up to energies of the order of 50 TeV (see Sect. 2). We obtain the particle distribution and γ-ray spectrum as described in Sect. 5, taking into account the radial dependence of the cluster photon field for computing losses and γ-ray emission. Particles are injected continuously over 3 Myr within a sphere of radius 10 pc centred on Cygnus OB2. The injection spectrum is a power-law with an exponential cut-off at p_max = 50 TeVc⁻¹ and a spectral index s = 4.1. The injection momentum is set to 10 MeVc⁻¹, as in the supernova model. The choice for the spectral index is motivated by the observation of comparatively steep γ-ray spectra from star-forming environments (e.g. H.E.S.S. Collaboration 2021, 2024; LHAASO Collaboration 2025). Figure 8 shows the stellar-wind inverse Compton component in yellow. The model reproduces the hardening of the spectrum in the Fermi-LAT band and the softening above 100 GeV. The γ-ray morphology is expected to peak at Cygnus OB2 and extend to a radius corresponding to the distance travelled by electrons within their cooling time (see Fig. 6), which is ~50 pc in our model. This is consistent with the morphology observed by LHAASO Collaboration (2024a) in the 2–20 TeV band.

To reproduce the flux level observed in γ-rays, the required power in non-thermal electrons is 3.5 × 10³⁶ erg s⁻¹ across all energies, which represents about 1% of the currently inferred stellar wind power in Cygnus OB2. While this value is higher than that typically inferred from multi-wavelength models of supernova remnants (e.g. Berezhko et al. 2009), the shock conditions at termination shocks of isolated stellar winds are quite different to those at the external shocks of young supernova remnants. For example, they have a rather slow upstream flow speed and high upstream temperature compared to supernova remnants. Stellar wind termination shocks are also expected to show high magnetic field obliquity. Furthermore, since electrons are injected continuously over the lifetime of the association, the expected flux depends on past stellar wind activity. It is conceivable that the wind power of Cygnus OB2 was higher in the past than the currently inferred 3 × 10³⁸ erg s⁻¹, considering the complex star-formation history of the association (see Sect. 2.1).

The bottom panel of Fig. 8 shows the contributions from scattering on different photon fields. In the Fermi-LAT band, the scattering on the cluster photon field dominates. The TeV range is dominated by scattering on dust emission. Figure 12 shows the multi-wavelength spectrum, including radio upper limits for the Cygnus Cocoon (Abeysekara et al. 2021) and X-ray upper limits from Suzaku. For our assumption of a homogeneous magnetic field of 5 µG, the synchrotron model remains within the upper limits, while bremsstrahlung is always subdominant in the low-density superbubble.

Our choice of parameters provides a good fit to the LHAASO and Fermi-LAT data in the ~1 GeV–10 TeV range while remaining within the multi-wavelength limits. Considerable uncertainty stems from the fact that the magnetic field is not well constrained and most certainly varies within the region. Due to the non-trivial nature of stellar wind interaction (Vieu et al. 2024b; Härer et al. 2025), it is challenging to make an informed choice for a prescription for the spatial dependence of the magnetic field. It is conceivable that the magnetic field in the inverse Compton emission region (radii ≲50 pc) is larger than 5 µG but decreases beyond to result in the given radio upper-limits. Increased synchrotron losses soften the spectrum above 1 TeV. Best-fit models for B ≳ 10 µG require either harder injection spectra (s = 4) or increased dust emission. As mentioned we already assume that dust emission is enhanced by a factor of two over the model by Popescu et al. (2017) as is appropriate in spiral arms. It is however conceivable that the enhancement factor is higher in the vicinity of Cygnus OB2 due to local dust production in the star-forming environment (see Green et al. 2019). Models with B < 5 µG provide a good fit for softer particle injection spectra, but tend to over-predict the γ-ray flux at 100 MeV–1 GeV Alternatively, the superbubble flow velocity could have a profile shallower than r⁻², in which case adiabatic losses would reduce the flux back within range. This is one of several factors that limits the predictive power of the model at low energies. The model also depends on the choice of lowest particle momentum (10 MeV c⁻¹ here) and the particle injection history over several millions of years. Cygnus OB2 harbours several distinct stellar populations, with ages in the range of 3–7 Myr (see Vieu et al. 2024b). Particle injection was almost certainly not continuous. Jointly, these factors can lead to a more complex spectral shape than predicted by our model. The fact that the bump around 1 GeV is not captured by the model is therefore not surprising. We note that assumptions on the low-energy electron population also impact synchrotron emission in the radio.

Fig. 12 Spectral energy distribution model showing additional components: synchrotron emission (purple) and bremsstrahlung (cyan) from electrons accelerated at stellar-wind termination shocks in Cygnus OB2, hadronic emission from the same source (grey), and inverse Compton emission from the supernova remnant shock (pink). The latter is normalised assuming a supernova energy of 3 × 1051 erg and an acceleration efficiency of electrons of 1%. We show radio upper limits from Abeysekara et al. (2021), which were derived for the entirety of the Cygnus Cocoon. We assume uniform B = 5 µG and n = 0.05 cm−3 in the superbubble.

7.2 Hadronic component

The grey lines in Fig. 12 show γ-rays produced in molecular clouds and in the neutral medium by protons accelerated at stellar-wind termination shocks in Cygnus OB2. We assume a hadronic acceleration efficiency of 10%. The model is dominated by the flux from hadronic interactions within the superbubble interior, because particles diffuse away quickly once they leave the superbubble and enter the interstellar medium. Therefore, they only produce low levels of γ-ray flux within the higher-density regions beyond the cavity. Simultaneously, the continuous supply of fresh particles from Cygnus OB2 ensures that the density of particles inside the superbubble remains high. The model is below the measured flux by about one order of magnitude. A fit to the data would require an unusually high superbubble interior density, n_SB>0.5 cm⁻³. For instance, previous claims that Cygnus OB2 produces hadronic γ-rays matching the observed flux assume densities of 20–30 cm⁻³ (Abeysekara et al. 2021; Bykov & Kalyashova 2022; Banik & Ghosh 2022; Menchiari et al. 2024). Within the framework of Weaver et al. (1977) n_SB > 1 cm⁻³ is only obtained assuming Cygnus OB2 was born in a parent molecular cloud of density >10³ cm⁻³. This is two orders of magnitude larger than the mean value estimated from the HI column density.

8 Conclusion

We have demonstrated that a powerful supernova that exploded about 50 kyr ago in the stellar association Cygnus OB2 is the best candidate for the origin of the diffuse γ-ray emission at ~10 TeV–PeV observed in the Cygnus region. Our model reproduces the observed energy dependence of the γ-ray morphology. With an additional inverse Compton component from electrons accelerated at stellar-wind termination shocks in Cygnus OB2, we modelled the full γ-ray spectrum from 100 MeV–PeV. This joint lepto-hadronic scenario is based on the following considerations:

1. Cygnus is a vast, close-by star-forming complex. Multi-wavelength observations reveal large-scale structures shaped by feedback from generations of massive stars. The main source of mechanical power is the loose stellar association Cygnus OB2. Even though the connection between large-scale structures observed at different wavelengths is far from being understood, they support the notion of a low-density superbubble about 150 pc in radius around Cygnus OB2;

2. Cygnus harbours powerful objects, such as Wolf-Rayet stars, pulsars, and supernova remnants. None of the detected objects stand out as a clear candidate source of ultra-high-energy γ-ray emission. We argue that an undetected supernova with above-average power has likely occurred within the last few 100 kyr; By extending the 3D hydrodynamical simulations presented in Vieu et al. (2024b), we determined the maximum particle energy at the forward shock of such a supernova remnant to be 1–2 PeV;

3. We considered a two-zone diffusion model, where the first zone represents the highly turbulent superbubble and the second zone represents the interstellar medium; A low-diffusion zone is required to confine the non-thermal particles within the γ-ray emission region. We obtained the particle distribution function by numerically solving the transport equation in spherical symmetry. Through reviewing the literature, we identified six molecular clouds within 100–200 pc from Cygnus OB2, which serve as target material for γ-ray production by non-thermal protons. In our preferred model, these clouds are located about equiradially at the edge of the Cygnus OB2 superbubble.

In the proposed scenario, an old supernova remnant would still be present at the current time, which due to having reached a size of several 10 s of parsecs in a hot low-density superbubble may have escaped detection. Nevertheless, its past feedback would not only impact the production of very-high-energy cosmic rays but would also have implications for our understanding of the dynamics of star-forming complexes. In the case of Cygnus, this might be key to understanding how the complex was carved up to large scales. As mentioned in Sect. 2, it is also not excluded that the very-high-energy particles were accelerated during the early expansion phase of the γ-Cygni supernova remnant, which despite its offset might be embedded in the Cygnus OB2 super-bubble. Very-high-energy particles could have leaked out of the remnant, filled the superbubble, and reached the complex of molecular clouds in Cygnus X. Determining under which conditions γ-Cygni can contribute to the non-thermal content of the Cygnus OB2 superbubble at very-high energy requires a detailed, asymmetric propagation model, but constructing such a model is beyond the scope of the present work.

While a powerful supernova can account for the γ-ray photons observed at ultra-high energies, an association of all flux above 1 PeV with such an event is not the only plausible scenario. Some ultra-high-energy photons are spatially correlated with the microquasar Cygnus X-3, and the contribution from diffuse galactic cosmic rays at ultra-high energies is uncertain. Taking these factors into account allows for a cut-off in the supernova remnant γ-ray spectrum below 2 PeV; in other words, the existence of an extreme accelerator in Cygnus X is not essential to explain the observations. At hundreds of TeV, the supernova remnant scenario is preferred over the wind termination shock scenario since Cygnus OB2 is too extended to form a collective cluster wind and the maximum particle energy at individual stellar-wind termination shocks is at most 100 TeV. Particles accelerated at wind shocks are however the preferred source of γ-rays below ~ 10 TeV. Because of the continuous supply of freshly accelerated particles, strong photon fields, and the low-density superbubble environment, inverse Compton dominates over hadronic emission and allows to explain the data with an electron acceleration efficiency of 1%.

Given the current statistical uncertainties, future observations with LHAASO and CTAO are necessary to constrain the origin of ultra-high-energy γ-ray photons observed in Cygnus. In addition, a dedicated multi-wavelength effort is needed to decipher the connection between large-scale structures across the electromagnetic spectrum in the Cygnus region. The relationship between measurements in the ultra-high-energy γ-ray regime and the role of massive star forming regions in the galactic cosmic-ray ecosystem remains an open problem. In this work, we have established a connection between an energetic supernova in Cygnus OB2 and multi-PeV cosmic rays.

Data availability

The column density maps for the ten molecular clouds listed in Table 1 are available at https://doi.org/00.17617/3.G9LUS1. They were inferred using the CO data and the catalogue of cloud distances from Zhang et al. (2024) and are available at https://doi.org/10.57760/sciencedb.16716 and https://www.doi.org/10.26093/cds/vizier.5167S22S, respectively. The initial conditions and output of the hydrodynamic simulation (Sect. 3.2) may be shared on reasonable request to TV. Details on the transport equation solver (Sect. 5.3) may be shared on reasonable request to FS (florian.schulze@mpi-hd.mpg.de) or TV.

Acknowledgements

This work made use of the MHD code PLUTO. Computations were performed on the HPC system Raven at the Max Planck Computing and Data Facility. We thank L. Olivera Nieto and Q. Remy for numerous helpful discussions on the γ-ray data and the gas distribution. We acknowledge R. Yang for providing LHAASO significance maps and details on the LHAASO analysis. CJKL gratefully acknowledges support from the International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg in the form of an IMPRS PhD fellowship.

Appendix A Distance to Cygnus OB2

Cygnus OB2 is by far the most massive and powerful association of massive stars in the Cygnus region. The determination of its distance is critical to understanding feedback processes on nearby molecular clouds and surrounding multi-wavelength emission, from radio to ultra-high-energy γ-rays. Spectrophotometric studies usually show a large spread in the distance modulus of the association and do not allow positioning of the association accurately within the Cygnus-X region (e.g. Hanson 2003).

Figure A.1 compiles distance estimates found in the literature. Very long baseline interferometry (VLBI) measurements of masers in Cygnus North suggest a distance of 1.3–1.5 kpc for these molecular clouds (Rygl et al. 2012); however, this might not necessarily be the distance to Cygnus OB2. VLBI measurements of circumstellar masers around the red hypergiant NML Cyg provide a distance of 1.61 ± 0.12 kpc, but NML Cyg, with a latitude of −1.9° , does not exactly correlate with Cygnus OB2; its projected separation is about 90 pc.

Kiminki et al. (2015) measure the distance to four spectroscopic binaries at 1.33 ± 0.6 kpc. However, Gaia EDR3 provides compatible parallaxes for three of these binaries, which are consequently no longer considered to be in the main group of Cygnus OB2. Only for one binary, MT91 372, there is tension in the distance determination (1765 ± 44 pc according to Gaia EDR3, still placing it within Cygnus OB2). A dedicated study of the eclipsing blue-supergiant binary V^★ V1827 Cyg finds a distance of 1.65 ± 0.15 kpc, in good agreement with the Gaia EDR3 parallax measurement (1765 ± 70 pc). This star is not close to the centre of Cygnus OB2, but at a projected distance of 0.45° (13 pc at a distance of 1.7 kpc). Still, it is considered to be part of the association in modern censuses.

Despite the systematic bias, Gaia parallaxes seem to provide the most direct, robust and accurate distance measurements. The latest analyses (Quintana & Wright 2021; Orellana et al. 2021) are compatible with distances of 1.6–1.75 kpc. Three stellar fore-ssociated with the open cluster UBC 585, and two at around 1.28 kpc.

Fig. A.1 Distance estimates for Cygnus OB2 quoted in the literature. Distances to foreground substructures are shown in fading colours. A systematic uncertainty of 0.03 mas is added to quoted Gaia parallaxes.

Appendix B Supplementary γ-ray maps

Figures B.1 and B.2 show the synthetic γ-ray emission for the supernova remnant model (see Sect. 6). Figure B.1 shows components of the model presented in Fig. 9. The figure demonstrates that the diffuse neutral gas component (HI) is required to reproduce the morphology of the emission, since emission from the clouds alone is too localised. In addition, the figure highlights that the HI component is brighter at lower energies. This is discussed in Sect. 6 and also shown in Fig. 10. Figure B.2 shows the model at 80 kyr. While the morphology in the >100 TeV band is unchanged compared to the model with t_SN = 50 kyr, the 25–100 TeV band shows decreased brightness in the northern clouds at l > 80°. As a consequence, the northern clouds are only marginally brighter in the lower-energy band than in the higher-energy band and match the LHAASO observations less well (see Fig. 9).

Fig. B.1 Components of the synthetic γ-ray emission maps shown in Fig. 9. The colour-scale is individually normalised to the maximum in each energy band. Maps are smoothed using a Gaussian function with a half-width of 0.3°.

Fig. B.2 Model at 80 kyr, which shows that older supernova remnants produce less prominent energy dependence in the northern clouds at l > 80° (cf. Fig. 9). The colour-scale is individually normalised to the maximum in each energy band. Maps are smoothed using a Gaussian function with a half-width of 0.3°.

References

All Tables

All Figures

Fig. 1 Powerful objects in the Cygnus region. Stars represent Wolf-Rayet stars from Rosslowe & Crowther (2015, http://pacrowt.her.staff.shef.ac.uk/WRcat/index.php). Diamonds indicate pulsars. The distance to PSR J2032+4127 has been revised to 1.76 kpc by association with its Be star companion MT91 213. Rings show the location of supernova remnants from Green (2024, 2025, https://www.mrao.cam.ac.uk/surveys/snrs). Different methods give distances to γ-Cygni between 0.98 to 2.3 kpc. The circles labelled ‘10 pc’ and ‘100 pc’ are centred on Cygnus OB2, assuming a distance of 1.65 kpc to the association.

In the text

	Fig. 2 Event rates for different types of core-collapse supernovae for a cluster with a total initial mass of 1.65 × 104 M⊙, as appropriate for Cygnus OB2, computed using Hoki.
In the text

	Fig. 3 Density slices from the simulation of a supernova remnant expanding in Cygnus OB2. (For details, see Sect. 3.2.)
In the text

Fig. 4 Top: supernova remnant properties across the forward shock surface at 2 kyr, extracted from the hydrodynamical simulation shown in Fig. 3. Bottom: integrated maximum energy at 2 kyr for the resonant and non-resonant streaming instabilities (RSI and NRSI, respectively) and without streaming instabilities. Beq was computed assuming that 10% of the kinetic energy in the superbubble goes into a turbulent magnetic field.

In the text

Fig. 5 Main molecular clouds in the Cygnus region. Red: Cygnus 1.3 North at l > 80° (Schneider Group II, 1.3 kpc) and Cygnus 1.3 South at l < 80° (part of Schneider Group IV, 1.3 kpc). Orange: Cygnus North Filament at l > 80° (Schneider Group I, 1.5 kpc) and Cygnus SFR South at l < 80° (part of Schneider Group IV, unknown distance >1.4 kpc). Pink: Cygnus South Filament (part of Schneider Group IV, unknown distance > 1.4 kpc). Green: Cygnus Central Clumps (Schneider Group III, 1.6–1.8 kpc). Dark blue at high latitude: Cygnus West, unknown distance (>1.4 kpc). Turquoise: S106 regions (described in Schneider et al. 2007, 1.6–1.8 kpc). The grey background shows foreground clouds: the North-America and Pelican nebulae on the left, the North-West clouds, and the Rift lane extending to the right in the direction to the Galactic Centre.

In the text

Fig. 6 Loss timescales for electrons in the putative superbubble surrounding Cygnus OB2. The photon field from stars in Cygnus OB2 is evaluated at 20 pc. Diffuse dust and starlight are taken from the model by Popescu et al. (2017) and enhanced by a factor of two. Coulomb losses are shown for a hot, low-density, ionised medium (T = 106 K, nSB = 0.05 cm−3) following Schlickeiser 2002 (Sect. 4.5). All other timescales are computed with the GAMERA code (Hahn et al. 2022).

In the text

	Fig. 7 Neutral gas model. The model consists of two 3D-Gaussian components, the Cygnus Rift (left) and the Cygnus SFR (middle), with a spherical (radius 150 pc), low-density (0.05 cm−3) cavity around Cygnus OB2. The total mass is set to match the mass obtained from integrating the neutral hydrogen column density obtained by Astiasarain et al. (2023). For model parameters, see Table 2.
In the text

Fig. 8 Preferred model for the region-integrated γ-ray spectrum consisting of hadronic emission from a 50 kyr-old supernova remnant along with inverse Compton emission from electrons accelerated at stellar-wind termination shocks in Cygnus OB2. The top panel shows both components in red and yellow as well as emission from diffuse galactic cosmic rays in blue. Diffuse galactic cosmic rays are only shown above 300 GeV since they are included in the Fermi-LAT background model. The middle panel shows hadronic model components. The lower panel shows contributions from inverse Compton scattering on individual photon fields. The data are from ARGO-YBJ Collaboration (2014), Abeysekara et al. (2021), Abdollahi et al. (2020), and LHAASO Collaboration (2024a).

In the text

Fig. 9 Synthetic γ-ray emission maps (top row) and significance maps from LHAASO Collaboration (2024a) (bottom row). Both the model and data have been smoothed using a Gaussian kernel of σ = 0.3°. The top of the colour scale is set to the maximum and the bottom to 10% of the maximum. The white cross marks the position of the stellar association Cygnus OB2. The figure demonstrates that the model is in good agreement with the γ-ray morphology and its energy dependence.

In the text

Fig. 10 Proton distribution function as a function of radius at 100 TeV (blue-dashed line) and 1 PeV (red-solid line). The dash-dotted vertical line marks the superbubble radius, which separates the low-diffusion and interstellar-medium diffusion zones. The solid vertical lines indicate the positions of molecular clouds in the model. The relative normalisation between the two displayed proton distributions changes strongly at small radii, which results in energy dependent morphology (see Fig. 9).

In the text

Fig. 11 Dependence of the model spectrum on cut-off energy and age of the supernova remnant (left panel) and uncertainty in the flux of galactic cosmic rays (middle panel). The right-hand panel shows the spectrum for the case that two of the eight photons with energy > 1 PeV can be attributed to the microquasar Cygnus X-3. In addition to a change in spectral shape, the energy required to match the normalisation of the data increases for lower cut-off energies and older supernova remnants. The displayed models require 4.5, 6, and 2.1 × 1050 erg in non-thermal protons for an age of 80 kyr and a cut-off energy of 1 PeV and 3 PeV, respectively.

In the text

Fig. 12 Spectral energy distribution model showing additional components: synchrotron emission (purple) and bremsstrahlung (cyan) from electrons accelerated at stellar-wind termination shocks in Cygnus OB2, hadronic emission from the same source (grey), and inverse Compton emission from the supernova remnant shock (pink). The latter is normalised assuming a supernova energy of 3 × 1051 erg and an acceleration efficiency of electrons of 1%. We show radio upper limits from Abeysekara et al. (2021), which were derived for the entirety of the Cygnus Cocoon. We assume uniform B = 5 µG and n = 0.05 cm−3 in the superbubble.

In the text

	Fig. A.1 Distance estimates for Cygnus OB2 quoted in the literature. Distances to foreground substructures are shown in fading colours. A systematic uncertainty of 0.03 mas is added to quoted Gaia parallaxes.
In the text

	Fig. B.1 Components of the synthetic γ-ray emission maps shown in Fig. 9. The colour-scale is individually normalised to the maximum in each energy band. Maps are smoothed using a Gaussian function with a half-width of 0.3°.
In the text

	Fig. B.2 Model at 80 kyr, which shows that older supernova remnants produce less prominent energy dependence in the northern clouds at l > 80° (cf. Fig. 9). The colour-scale is individually normalised to the maximum in each energy band. Maps are smoothed using a Gaussian function with a half-width of 0.3°.
In the text