Deciphering the gamma-ray emission in the Cygnus region

TL;DR

This paper models gamma-ray emission in the Cygnus region, suggesting a supernova remnant explains high-energy observations, while inverse Compton processes fit lower energies, highlighting complex astrophysical sources.

Contribution

It introduces a comprehensive 3D modeling approach combining hydrodynamical simulations and gas distribution to interpret gamma-ray data in Cygnus.

Findings

A ~50 kyr-old supernova remnant explains gamma-ray morphology and spectrum from 10 TeV to PeV.

Inverse Compton emission from stellar-wind shocks fits gamma-ray data below 10 TeV.

Hadronic models are incompatible with the superbubble structure around Cygnus OB2.

Abstract

The Cygnus region is a vast star-forming complex harbouring a population of powerful objects, including massive star clusters and associations, Wolf-Rayet stars, pulsars, and supernova remnants. The multi-wavelength picture is far from understood, in particular the recent LHAASO detection of multi-degree scale diffuse gamma-ray emission up to PeV energies. We aim to model the broadband gamma-ray data, discriminating plausible scenarios amongst all candidate accelerators. We consider in particular relic hadronic emission from a supernova remnant expanding in a low-density environment and inverse Compton emission from stellar-wind termination shocks in the Cygnus OB2 stellar association. We first estimate the maximum particle energy from a 3D hydrodynamical simulation of the supernova remnant scenario. The transport equation is then solved numerically to determine the radial distribution of non-thermal protons and electrons. In order to compute synthetic gamma-ray spectra and emission maps, we develop a 3D model of the gas distribution. This includes, firstly, a HI component with a low-density superbubble around Cygnus OB2 and, secondly, molecular clouds lying at the edge of the superbubble and in the foreground. We find that a powerful, ~50 kyr-old supernova remnant can account for both the morphology and spectrum from 10 TeV-PeV. At PeV energies, the microquasar Cygnus X-3 and diffuse Galactic cosmic rays might also contribute to the flux. Below about 10 TeV, hadronic models are incompatible with the expected existence of a superbubble centred on Cygnus OB2. Instead, the spectrum is well fitted with inverse Compton emission from electrons accelerated at stellar-wind termination shocks in Cygnus OB2 in line with existing multi-wavelength limits.

Tables

Table 1. Table 1: Molecular clouds in the Cygnus region.

	Name	l [∘]	Distance [kpc]	Distance ref
A	Foreground Rift	74–86	0.7–1.1	Zhang et al. (2024)
B	Cygnus North-West	79.5–81	$\sim 1$	Zhang et al. (2024)
C	Cygnus 1.3 North	80–83.5	1.3–1.4	Zhang et al. (2024)
D	Cygnus 1.3 South	77–80	1.2–1.3	Zhang et al. (2024)
E	Cygnus North Filament	80.5–82	1.4–1.5	Rygl et al. (2012)
F	Cygnus Central Clumps	79–81	1.6–1.8	Orellana et al. (2021); Xu et al. (2013b)
G	Cygnus West	78–80	$>1.4$
H	Cygnus SFR South	77–80	$>1.4$
J	Cygnus South Filament	77–79	$>1.4$
K	Cygnus Far South	76–78	1.6–1.8	Schneider et al. (2007)

Table 2. Table 2: Parameters of the neutral gas model

Parameter	Rift	SFR
centre $d$ [kpc]	1.0	1.7
centre $l$ [∘]	78	79
centre $b$ [∘]	0.5	2
total mass [${10}^6{M}_{\odot }$]	4.2	3.4
${\sigma }_l$ [pc]	200	200
${\sigma }_b$ [pc]	50	80
${\sigma }_d$ [pc]	200	200
rotation [radian]	0.05	6

Equations

$$E_{\max }\approx \left(\frac{B}{5\text{leavevmode}\text{textmuG}}\right)\left(\frac{V_{\mathrm{SN}}}{10^4\text{leavevmode}\text{kms−1}}\right){\left(\frac{M_{\mathrm{ej}}}{1M_{\odot }}\right)}^{\frac{1}{3}}{\left(\frac{n_{\mathrm{SB}}}{0.01\text{leavevmode}{\mathrm{cm}}^{-3}}\right)}^{-\frac{1}{3}}\text{leavevmodenobreak PeV}.$$

$${\xi }_{\mathrm{cr}}{u}^3>v_{\mathrm{A}}^{2}c,$$

$$\delta B^{2}_{\rm sat}\approx\left\{\begin{array}[]{cc}2\xi_{\mathrm{cr}}u\sqrt{4\pi\rho}\,B_{\mathrm{eq}}&\enspace\mbox{RSI}\\ 4\pi\xi_{\mathrm{cr}}\,\rho u^{3}/c&\enspace\mbox{NRSI}\end{array}\right.\,,$$

$$\frac{\mathrm{d}p}{\mathrm{d}t}\approx \frac{u^{2}p}{20\kappa },\kappa =pc/(3qB),$$

$$R_{\text{b}}=183{\left(\frac{{\xi }_{\mathrm{b}}{L}_{\text{w}}}{10^{38}\mathrm{erg}{\mathrm{s}}^{-1}}\right)}^{1/5}{\left(\frac{n_{\mathrm{ISM}}}{1{\mathrm{cm}}^{-3}}\right)}^{-1/5}{\left(\frac{t}{5\mathrm{Myr}}\right)}^{3/5}\mathrm{pc},$$

$${\partial }_{t}f$$

Full text

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) $\gamma$-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, making 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” for “extended”). On top of the diffuse emission, Downes & Rinehart (1966) identified 20 individual HII regions, later assigned DR4–DR23, which appear as prominent radio hotspots. The only source of non-thermal (synchrotron) radio emission is the supernova remnant $\gamma$-Cygni (Xu et al., 2013a). The diffuse thermal component comes from material ionised by massive stars, in particular 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^{\circ}\times 13^{\circ}$) 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 (Uyanıker et al., 2001). Another hint of large-scale feedback in Cygnus is the presence of a ring of H$\alpha$ 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\mbox{--}3\times 10^{38}\,\mathrm{erg}\,\mathrm{s}^{-1}$, 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, as well as 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^{\circ}$) and the southern ($l<80^{\circ}$) 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, in particular 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, concluding 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 $\gamma$-ray sky. Extended GeV sources in Cygnus were first detected by EGRET (1995) and later associated with the $\gamma$-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 (2002) and MILAGRO (2007) and later associated with the pulsar wind nebula powered by PSR J2032+4127 (Camilo et al., 2009).

On top of these three identified $\gamma$-ray sources, the Fermi-LAT collaboration discovered an extended diffuse $\gamma$-ray source coincident with the radio emission in Cygnus X, filling a region delimited by bright infrared photodissociation regions (Fermi-LAT, 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 (2014), VERITAS (2018), Tibet {\mathrm{AS}}_{\gamma} (2021), HAWC (2021) and LHAASO (2024a). The centroid of the diffuse emission lies at the heart of Cygnus X, approximately coincident with the position of Cygnus OB2 below 1 TeV, and loosely correlated with the northern and southern molecular clouds at ultra-high energy. The spectrum is hard in the GeV band ($\mathrm{d}N/\mathrm{d}E\sim E^{-2}$) but turns over around 1 TeV ($\mathrm{d}N/\mathrm{d}E\sim E^{-3}$). The steep component measured by LHAASO (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 related to PSR J2032+4127 or $\gamma$-Cygni at their current stages of evolution, as the spectra of these sources cut off at energies ${\lesssim}100$ TeV and their integrated flux is subdominant with respect to the extended diffuse component (Fleischhack, 2019; HAWC, 2021; MAGIC, 2023; LHAASO, 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 model a cavity of radius 55 pc filled with a turbulent plasma of average velocity $1500\,\mathrm{km}\,\mathrm{s}^{-1}$, 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 is however unlikely to produce multi-PeV protons (see Sect. 2), and it has further 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 $\gamma$-ray emission in star-forming environments generally suffer from the fundamental caveat that they simultaneously need i) a low-density environment which allows 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 $\gamma$-ray emission in Cygnus was excluded by HAWC (2021) arguing that in a magnetic field of 20 µG and a density of 30 cm*-3* 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*-3* within a few tens of parsecs from Cygnus OB2. The values quoted by the HAWC collaboration are taken from Fermi-LAT (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 promotes consideration for a leptonic origin of the diffuse $\gamma$-ray emission, in particular 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) cannot however 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 $\gamma$-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, disentangling the foreground components. We then proceed to build a lepto-hadronic model to show how the diffuse $\gamma$-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, 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\mbox{--}3\times 10^{38}\,\mathrm{erg}\,\mathrm{s}^{-1}$ (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.

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^{38}\,\mathrm{erg}\,\mathrm{s}^{-1}$ (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 ${\sim}3\times 10^{36}\,\mathrm{erg}\,\mathrm{s}^{-1}$), which powers the Dragonfly pulsar wind nebula, PSR J2021+4026 (spin-down power ${\sim}10^{35}\,\mathrm{erg}\,\mathrm{s}^{-1}$), and finally PSR J2032+4127 (spin-down power ${\sim}10^{35}\,\mathrm{erg}\,\mathrm{s}^{-1}$), 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 $\gamma$-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 $\gamma$-ray emission, although it might contaminate the significance map around Cygnus OB2.

Of the three supernova remnants observed within 2 kpc toward the Cygnus region, W63 and the Cygnus Loop are older remnants that, while perhaps efficient accelerators in the distant past, are both too far away from the centre of emission to account for the integrated flux. $\gamma$-Cygni is a younger remnant (about 5 kyr, Lozinskaya et al., 2000) which is bright in radio, X-rays, and $\gamma$-rays. It expands in a low-density environment (about 0.1 cm*-3*) to a present size of about 15 pc with a shock velocity of nearly $1000\,\mathrm{km}\,\mathrm{s}^{-1}$ (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 $\gamma$-ray spectrum associated with $\gamma$-Cygni extends to tens of TeV (Fleischhack, 2019) and is thought to originate from proton-proton interactions with a nearby cloud (MAGIC, 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 $\gamma$-ray emitting particles were accelerated in the early expansion phase of the remnant and already escaped the $\gamma$-Cygni region to fill a network of cavities in the Cygnus star-forming complex, eventually reaching 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 $\gamma$-ray emission. Similar microquasars have been recently recognised as efficient particle accelerators up to 100s of TeV (LHAASO, 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. 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}_{\mathrm{acc}}\approx V_{\rm w}^{2}p/(4D)$, where $V_{\rm 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 $\dot{p}_{\mathrm{loss}}\geq 2V_{\mathrm{w}}p/(3R_{\rm WTS})$, where $R_{\rm WTS}$ is the radius of the wind termination shock. The particle energy is therefore bounded as $D(E)<3R_{\rm WTS}V_{\rm 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/(\xi qB)$ with $q$ the charge of the particle, $c$ the speed of light and $\xi\leq 3$. The maximum energy is then limited as $E_{\rm max}<qBR_{\rm WTS}V_{\rm 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_{\rm w}/M_{\rm A}$, where $\rho$ is the gas density and efficient particle acceleration requires $M_{\rm A}$ substantially exceeds unity. Conservation of mass in a spherical wind satisfies $\rho=\dot{M}/(4\pi R_{\rm WTS}^{2}V_{\mathrm{w}})$, which leads to $E_{\rm max}\ll q\sqrt{V_{\rm w}L_{\rm w}}/c$, where $L_{\rm w}=0.5\dot{M}V_{\mathrm{w}}^{2}$ is the mechanical power and $\dot{M}$ 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$\,\mathrm{km}\,\mathrm{s}^{-1}$, mass-loss rates around $5\times 10^{-6}\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}$, stellar radii around 15–20 $R_{\odot}$ and projected rotational velocities around $100\mbox{--}150\,\mathrm{km}\,\mathrm{s}^{-1}$ (Berlanas et al., 2020). Assuming a strong surface field ${\sim}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 ${\sim}200\,\mathrm{km}\,\mathrm{s}^{-1}$, surface magnetic field ${\sim}1$ kG and stellar radius ${\sim}5\,R_{\odot}$, 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\times 10^{37}\,\mathrm{erg}\,\mathrm{s}^{-1}$, but whose radius is only about $1\,R_{\odot}$ (Sander et al., 2019), the maximum energy falls short of the PeV regime even for extremely large surface fields or rotational velocities.

On top of 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_{\mathrm{sat}}^{2}\approx 4\pi\xi_{\mathrm{cr}}\,\rho V_{\mathrm{w}}^{3}/c$ which, for an acceleration efficiency $\xi_{\mathrm{cr}}=0.1$, a density $\rho=0.01m_{\mathrm{p}}\,\mathrm{cm}^{-3}$ and a wind velocity $V_{\rm w}=3000\,\mathrm{km}\,\mathrm{s}^{-1}$, gives $\delta B_{\mathrm{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^{38}\,\mathrm{erg}\,\mathrm{s}^{-1}$, 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\,\mathrm{MeV}<E<100$ GeV) to very high-energy ($100\,\mathrm{GeV}<E<100$ TeV) diffuse $\gamma$-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). 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 assume Solar metallicity and a single stellar population with a total initial mass of $1.65\times 10^{4}\,\mathrm{M}_{\odot}$ (Wright et al., 2015). We assume an initial mass function index $\alpha=-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 $\sim 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^{51}$ 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 kyr 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). Besides, the shock is never expected to enter the radiative phase but rather slow down to subsonic velocities after a few 100 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):

[TABLE]

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 $\gamma$-ray spectrum.

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 $\dot{M}=10^{-5}\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}$ and $V_{\mathrm{w}}=2500\,\mathrm{km}\,\mathrm{s}^{-1}$. 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 $\rho\propto r^{-9}$, as described in Whalen & Norman (2008) and Badmaev et al. (2024). The explosion energy is set to $5\times 10^{51}$ erg and the core velocity to 16 000$\,\mathrm{km}\,\mathrm{s}^{-1}$, 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*-3*), 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 15 pc and the velocity has decreased to about $3000\,\mathrm{km}\,\mathrm{s}^{-1}$, 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_{\mathrm{eq}}=u\sqrt{4\pi\eta_{\mathrm{B}}\rho}$ with $\eta_{\mathrm{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)

[TABLE]

where $u$ is the shock velocity, $v_{\mathrm{A}}$ the upstream background Alfvén speed, and $\xi_{\mathrm{cr}}$ is the acceleration efficiency that we set to 10%. The background field, $B_{\mathrm{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_{\mathrm{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)

[TABLE]

where RSI denotes the resonant and NRSI the non-resonant streaming instability. In all cases we assume that diffusive shock acceleration operates at maximum efficiency (Bohm diffusion), in which case the acceleration rate is written as (Lagage & Cesarsky, 1983)

[TABLE]

where $B$ is the turbulent magnetic field. We repeat the calculation assuming first that streaming instabilities do not operate ($B=B_{\mathrm{eq}}$), then that RSI operates ($B=\delta B_{\mathrm{sat,RSI}}$) and finally that NRSI operates ($B=\delta B_{\mathrm{sat,NRSI}}$ if the left hand-side of Eq. 2 is ten times the right hand-side, $B=B_{\mathrm{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. Note that we naively expect the maximum energy to scale as the square root of the explosion energy such that a nominal supernova ($10^{51}$ 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 magnetized stellar winds within superbubbles (see for example Härer et al., 2025). These findings provide confidence that the Hillas limit could be achieved without resorting to the Bohm scattering limit.

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 $\gamma$-ray emission (see Sect. 5 and Sect. 6).

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 which 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 S106 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 12CO and 13CO emission lines. Because the Cygnus region is seen tangentially to the galactic rotation, kinematic velocities do not translate into reliable distance measurements. Restricting velocities to the range $(-20,30)\,\mathrm{km}\,\mathrm{s}^{-1}$ 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^{\circ}$ in the velocity range $(-7,1)\,\mathrm{km}\,\mathrm{s}^{-1}$. 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)\,\mathrm{km}\,\mathrm{s}^{-1}$, and was placed in the foreground by early absorption studies (Dickel & Wendker, 1978). Parallax measurements of masers (Rygl et al., 2012) place DR20 and DR21 (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\mbox{--}0.8$ kpc. Schneider et al. (2007) completed the census in the southernmost part up to the S106 HII region and showed that several substructures in this area are shaped by stellar feedback in a star-forming region at $1.6\mbox{--}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.

Cygnus-X North

Unfortunately, the extinction becomes too high beyond ${\sim}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 ${\sim}1.6\mbox{--}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 12CO 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)\,\mathrm{km}\,\mathrm{s}^{-1}$ and $(7,15)\,\mathrm{km}\,\mathrm{s}^{-1}$, respectively. These are shown in red and orange in Fig. 5. 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”. 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)\,\mathrm{km}\,\mathrm{s}^{-1}$, $(2.2,4.5)\,\mathrm{km}\,\mathrm{s}^{-1}$ and $(4,8.5)\,\mathrm{km}\,\mathrm{s}^{-1}$ around $78\mbox{--}79^{\circ}$ in longitude. The result defines the “Cygnus SFR South” region shown in orange at $l<80^{\circ}$ and $-1^{\circ}\lesssim b\lesssim 1^{\circ}$ in Fig. 5. In the range $(12,20)\,\mathrm{km}\,\mathrm{s}^{-1}$, a filamentary cloud appears, which we call the “South Filament” (shown in pink in Figure 5). Finally, another region at high latitude, $b\approx 4^{\circ}$ in the velocity range $(-5,3)\,\mathrm{km}\,\mathrm{s}^{-1}$ 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 DR18 and DR20 in the velocity range $(-11,1)\,\mathrm{km}\,\mathrm{s}^{-1}$ around $l=81^{\circ}$. 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_{\mathrm{CO}}$ factor of $4\times 10^{20}$ H/cm2 (K km s*-1*)-1. For later convenience, we introduce shorthands A–K to refer to the clouds.

5 Modelling particle propagation and $\gamma$-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 $\gamma$-ray emission. We construct both hadronic and leptonic models and apply them to supernova and stellar-wind scenarios in Sect. 6 and 7.

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 (c.f. 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\,\mathrm{eV}\,\mathrm{cm}^{-3}$. 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\pi R^{2}c)$. At a distance of 20 pc, the energy density of the cluster photon field is $25\,\mathrm{eV}\,\mathrm{cm}^{-3}$, 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\,\mathrm{cm}^{-3}$, 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.

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 Tab. 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 Tab. 2. Model parameters and total masses are 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

[TABLE]

where $L_{\mathrm{w}}$ is the wind power, $n_{\mathrm{ISM}}$ the ambient density, $t$ the age of the system, and $\xi_{\mathrm{b}}$ 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^{\circ}$, $b=0.8^{\circ}$). Inside the superbubble, we set the density to a constant value of $0.05\,\mathrm{cm}^{-3}$. This value lies in the range typical for superbubbles (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 evaporation at the shell is hence 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_{\rm SB}$, is expected to be smaller than its average interstellar value. We assume $D_{\rm SB}=10^{25}p_{\rm GeV}^{1/2}$ cm2 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_{\rm ISM}=10^{28}p_{\rm GeV}^{0.4}$ cm2 s*-1*.

We model large-scale advection in a mean $r^{-2}$ wind profile with a speed of $100\,\mathrm{km}\,\mathrm{s}^{-1}$ 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 $f(r,p,t)$ reads,

[TABLE]

where we adopt spherical symmetry. Adiabatic losses vanish for our adopted $r^{-2}$ radial velocity profile. $Q(r,p,t)$ is a source term which 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 $f(r,p,t)$ which can be post-analysed to determine the non-thermal emission.

5.4 Gamma-ray emission

We compute non-thermal emission using the GAMERA code (Hahn et al., 2022)333http://libgamera.github.io/GAMERA/docs/main_page.html. In Sect. 4, we provided the details of our model for the 3D distribution of target material, which includes both the molecular clouds and neutral medium. Convolving this with $f(r,p,t)$ yields both the full spectra integrated over the region and synthetic $\gamma$-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 $\gamma$-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\,\mathrm{cm}^{-3}$. 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 Setup 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 $\gamma$-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\propto p^{-s}\exp(-p/p_{\mathrm{max}})$, where $p_{\mathrm{max}}=2\,\mathrm{PeV}c^{-1}$, the injection momentum is set to 10 $\mathrm{MeV}c^{-1}$ 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 $\sigma=0.3^{\circ}$ 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\times 10^{50}\,$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 $\gamma$-ray morphology between the 25–100 TeV and ${>}100\,$TeV bands (see Fig. 15 in the appendix). We note that the parameters $t_{\mathrm{SN}}$ and $D_{\mathrm{SB}}$ are degenerate as the proton distribution scales with $\sqrt{D_{\mathrm{SB}}t_{\mathrm{SN}}}$. A different assumption on $D_{\mathrm{SB}}$ will therefore affect the age estimate.

The required energy corresponds to a powerful supernova with $3\times 10^{51}\,$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\,\mathrm{cm}^{-3}$, 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 superbubble to maintain the relative normalisation between the neutral gas and molecular cloud component, which is constrained by the $\gamma$-ray morphology (see Fig. 14 in the appendix). For typical superbubble densities, $n_{\mathrm{SB}}\lesssim 0.1\,\mathrm{cm}^{-3}$, the model predicts the clouds to be within ${\sim}10\mbox{--}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 $\gamma$-ray maps. We highlight however 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 $\gamma$-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.

6.3 Origin of the ultra-high energy $\gamma$-rays

Substantial interest has been taken in the origin of the ultra-high energy $\gamma$-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_{\mathrm{max}}=2\,$PeV and $t_{\mathrm{SN}}=50$ kyr gives a model within the $1\sigma$ 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_{\mathrm{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^{51}$ erg assuming $E_{\mathrm{max}}$ is roughly proportional to the square root of the explosion energy). However, as argued in Sec. 3.1, the young age of Cygnus OB2 favours Type Ic supernovae whose energy likely exceeds $10^{51}$ 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\gtrsim 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 (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\sigma$.

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\sigma$ uncertainties of the data above 1 PeV, if the injection spectrum is a $p^{-4}$ power-law with exponential cut-off.

6.4 Sensitivity to other parameters

We set $D_{\rm ISM}=10^{28}p_{\rm GeV}^{0.4}\mathrm{cm}^{2}\,\mathrm{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_{\mathrm{ISM}}\gg D_{\mathrm{SB}}$, which remains true even if $D_{\mathrm{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\mbox{--}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 $\gamma$-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 $\gamma$-ray spectrum as described in Sect. 5, taking into account the radial dependence of the cluster photon field for computing losses and $\gamma$-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_{\mathrm{max}}=50\,\mathrm{TeV}c^{-1}$ and a spectral index $s=4.1$. The injection momentum is set to 10 $\mathrm{MeV}c^{-1}$, as in the supernova model. The choice for the spectral index is motivated by the observation of comparatively steep $\gamma$-ray spectra from star-forming environments (e.g. H.E.S.S., 2021, 2024; LHAASO, 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 $\gamma$-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 ${\sim}50\,$pc in our model. This is consistent with the morphology observed by LHAASO (2024a) in the 2–20 TeV band.

To reproduce the flux level observed in $\gamma$-rays, the required power in non-thermal electrons is $3.5\times 10^{36}\,\mathrm{erg}\,\mathrm{s}^{-1}$ 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\times 10^{38}\,\mathrm{erg}\,\mathrm{s}^{-1}$, 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 (HAWC, 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 ${\sim}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 ${\lesssim}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\gtrsim 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 $\gamma$-ray flux at 100 MeV–1 GeV. Alternatively, the superbubble flow velocity could have a profile shallower than $r^{-2}$, 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*-1* 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.

7.2 Hadronic component

The grey lines in Fig. 12 show $\gamma$-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 $\gamma$-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_{\mathrm{SB}}>0.5\,\mathrm{cm}^{-3}$. For instance, previous claims that Cygnus OB2 produces hadronic $\gamma$-rays matching the observed flux assume densities of $20\mbox{--}30$ cm*-3* (HAWC, 2021; Bykov & Kalyashova, 2022; Banik & Ghosh, 2022; Menchiari et al., 2024). Within the framework of Weaver et al. (1977) $n_{\mathrm{SB}}>1\,\mathrm{cm}^{-3}$ is only obtained assuming Cygnus OB2 was born in a parent molecular cloud of density ${>}10^{3}$ cm*-3*. This is two orders of magnitude larger than the mean value estimated from the HI column density.

8 Conclusion

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

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. 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 $\gamma$-ray emission. We argue that an undetected supernova with above-average power has likely occurred within the last few 100 kyr. Extending the 3D hydrodynamical simulations presented in Vieu et al. (2024b), we determine the maximum particle energy at the forward shock of such a supernova remnant to be 1–2 PeV. 3. 3.

We consider a two-zone diffusion model, where the first zone represents the highly turbulent superbubble and the second zone the interstellar medium. A low-diffusion zone is required to confine the non-thermal particles within the $\gamma$-ray emission region. We obtain the particle distribution function by numerically solving the transport equation in spherical symmetry. Reviewing the literature, we identify six molecular clouds within 100–200 pc from Cygnus OB2, which serve as target material for $\gamma$-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 remain at the present day, which, reaching a size of several 10s 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 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 Sec. 2, it is also not excluded that the very-high energy particles were accelerated during the early expansion phase of the $\gamma$-Cygni supernova remnant which, despite its offset, might be embedded in the Cygnus OB2 superbubble. 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 $\gamma$-Cygni can contribute to the non-thermal content of the Cygnus OB2 superbubble at very-high energy requires a detailed, asymmetric, propagation model, which is beyond the scope of the present work.

While a powerful supernova can account for the $\gamma$-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 $\gamma$-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 $\gamma$-rays below ${\sim}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 $\gamma$-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 $\gamma$-ray regime and the role of massive-star forming regions in the galactic cosmic-ray ecosystem remains an open problem. In this work, we 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/10.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.51670220, 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 $\gamma$-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 $\gamma$-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 13 compiles distance estimates found in the literature. Very Long Baseline Interferometry (VLBI) measurements of masers in Cygnus North suggest a distance of $1.3\mbox{--}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\pm 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\pm 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\pm 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\pm 0.15$ kpc, in good agreement with the Gaia EDR3 parallax measurement ($1765\pm 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\mbox{--}1.75$ kpc. Three stellar foreground groups are identified by Orellana et al. (2021), one at around 1.46 kpc, which is associated with the open cluster UBC 585, and two at around 1.28 kpc.

Appendix B Supplementary $\gamma$-ray maps

Figures 14 and 15 show the synthetic $\gamma$-ray emission for the supernova remnant model (see Sect. 6). Figure 14 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 15 shows the model at 80 kyr. While the morphology in the ${>}100\,$TeV band is unchanged compared to the model with $t_{\mathrm{SN}}=50\,$kyr, the 25–100 TeV band shows decreased brightness in the northern clouds at $l>80^{\circ}$. 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).