Hydrodynamic simulation of Cygnus OB2: the absence of a cluster wind termination shock

T. Vieu,1 C. J. K. Larkin,1,2,3 L. Härer,1 B. Reville,1 A.A.C. Sander,2 V. Ramachandran2
1Max-Planck-Institut für Kernphysik, Saupfercheckweg 1, D-69117 Heidelberg, Germany
2Zentrum für Astronomie der Universität Heidelberg, Astronomisches Rechen-Institut, Mönchhofstr. 12-14, 69120 Heidelberg, Germany
3Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany
E-mail: thibault@mpi-hd.mpg.de
(Accepted XXX. Received YYY; in original form ZZZ)
Abstract

We perform a large-scale hydrodynamic simulation of a massive star cluster whose stellar population mimics that of the Cygnus OB2 association. The main-sequence stars are first simulated during 1.6 Myr, until a quasi-stationary state is reached. At this time the three Wolf-Rayet stars observed in Cygnus OB2 are added to the simulation, which continues to 2 Myr. Using a high-resolution grid in the centre of the domain, we can resolve the most massive stars individually, which allows us to probe the kinetic structures at small (parsec) scales. We find that, although the cluster excavates a spherical “superbubble” cavity, the stellar population is too loosely distributed to blow a large-scale cluster wind termination shock, and that collective effects from wind-wind interactions are much less efficient than usually assumed. This challenges our understanding of the ultra-high energy emission observed from the region.

keywords:
open clusters and associations: individual: Cygnus OB2 – shock waves – stars: winds, outflows – hydrodynamics – cosmic rays
pubyear: 2024pagerange: Hydrodynamic simulation of Cygnus OB2: the absence of a cluster wind termination shockHydrodynamic simulation of Cygnus OB2: the absence of a cluster wind termination shock

1 Introduction

The Cygnus X complex is a nearby region of ongoing star-formation at a distance of about 1.6 kpc. It hosts a number of OB associations, including Cygnus OB2, a young massive star cluster containing thousands of OB stars (Massey & Thompson, 1991; Knödlseder, 2000) among which tens are known to be powerful early O-stars Wright et al. (2015); Berlanas et al. (2020). The strong mechanical feedback imparted by the cluster is expected to shape the surrounding environment and impact stellar formation. This in turn is anticipated to influence non-thermal particle acceleration and emission from the cluster. However, despite decades of observational and theoretical investigations, there is still a limited understanding of the physical processes that account for observations of the region. One of the chief difficulties is the presence of several foreground structures which hinder the interpretation of multi-wavelength observations (Uyanıker, B. et al., 2001; Zhang et al., 2024) and the reconstruction of the gas distribution. Radio and X-ray observations hint at the presence of a large-scale bubble (e.g., Albacete-Colombo et al., 2023), most likely heated by the stellar feedback, although the boundaries do not seem to coincide with the dense “shell” structure that is theoretically expected (Weaver et al., 1977).

The Cygnus X region has become a source of great interest in the field of high-energy astrophysics since the discovery of coincident diffuse gamma-ray emission first measured by the Fermi satellite (Ackermann et al., 2011), and subsequently detected by the HAWC and LHAASO observatories (Abeysekara et al., 2021; Cao et al., 2021). In particular, the LHAASO collaboration reported the detection of gamma-rays above 1 PeV, which suggests that the Cygnus region hosts a hadronic “PeVatron” (LHAASO Collaboration, 2024). Identifying and understanding acceleration sites within the Cygnus X region could unveil a new class of cosmic-ray sources which, given the growing number of observational constraints challenging the classical supernova paradigm (e.g., Gabici et al., 2019), is key to understanding the origin of galactic cosmic rays at very-high energies.

Star-forming regions are intricate environments with non-trivial feedback amongst individual components. The massive stars blow powerful supersonic winds which create, in the early-stage of the cluster evolution, a collection of expanding cavities called wind-blown bubbles (Weaver et al., 1977). After a few tens of kyrs, the wind-blown bubbles percolate to form a large-scale superbubble filled with a hot, high pressure and rarefied plasma. Inside the superbubble, each massive star continues to blow a supersonic wind, which terminates at a distance of approximately 1 pc away from the star, at its stellar wind termination shock. If the cluster is compact enough (typically if its half-mass radius is of the order of 1 pc), the stellar winds interact so efficiently that they are expected to merge into a single, large-scale (1020102010-2010 - 20 pc) supersonic outflow, which terminates at a so-called cluster wind termination shock. On the other hand, if the stars are too loosely distributed, the interactions between the supersonic winds might not be enough to drive the formation of such a large-scale cluster wind termination shock (Gupta et al., 2020). In the end, the environment in and around a stellar cluster might contain several strong shocks from the interactions between supersonic outflows, which is also expected to generate a high level of turbulence (Vieu et al., 2024). This makes these regions favourable sites for particle acceleration.

Several scenarios of efficient particle acceleration in star-forming regions have been proposed in the past, including re-acceleration processes in the core (e.g. Bykov & Fleishman, 1992; Klepach et al., 2000; Bykov et al., 2013), re-acceleration processes in the low-density superbubble (e.g. Ferrand & Marcowith, 2010; Vieu et al., 2022a), acceleration at the cluster wind termination shock (e.g. Morlino et al., 2021), acceleration by a supernova remnant expanding in the close vicinity of the cluster core (Vieu & Reville, 2023). However, most of these scenarios fail to accelerate PeV protons (Vieu et al., 2022b), which are necessary to explain the LHAASO data. Shortcomings of the re-acceleration processes in the core have been especially addressed in recent work (Vieu et al., 2024), where it was shown that re-acceleration processes become inefficient in expected cluster-core conditions. It seems therefore that a scenario such as the one investigated by Bykov & Kalyashova (2022), where particles are re-accelerated in an ensemble of stochastic shocks interacting in a very compact cluster core, is excluded for Cygnus OB2. The presence of turbulence in the low density cavity beyond the cluster core might enhance the acceleration processes, although this is not expected to produce PeV protons (Vieu et al., 2022b), especially in the case of Cygnus OB2 which, with an estimated mechanical power of a few 1038superscript103810^{38}10 start_POSTSUPERSCRIPT 38 end_POSTSUPERSCRIPT erg s-1, is not the most powerful stellar cluster in our Galaxy.

Another possible scenario to account for the gamma-ray observations is that of efficient particle acceleration at a large-scale cluster wind termination shock (WTS). Menchiari et al. (2024) recently explored this possibility and argued that it could provide a satisfactory explanation to the data, although the authors disregarded the latest LHAASO measurements. This model seems promising, however it entirely relies on the hypothesis that Cygnus OB2 is able to blow a large-scale cluster WTS, which is as yet unclear.

With an estimated power of 23×103823superscript10382-3\times 10^{38}2 - 3 × 10 start_POSTSUPERSCRIPT 38 end_POSTSUPERSCRIPT erg s-1 and a age of a few Myr (Wright et al., 2015), the one-dimensional hydrodynamic theory of Weaver et al. (1977) predicts a WTS size of about 13 pc, which could in principle be enough to accelerate PeV protons if Bohm diffusion is assumed (Morlino et al., 2021). However, Cygnus OB2 is a rather peculiar stellar cluster, closer to a loose association of massive stars than a young compact cluster such as Westerlund 1. Indeed, the stellar power is not concentrated in a compact region, but distributed within about 15 pc (in particular, the three Wolf-Rayet stars, which dominate the mechanical output, are strongly off-centred). In this case, the putative cluster WTS is a minima weak. Using the calculations by Cantó et al. (2000), one finds a sonic Mach number of 2.6, which, in a textbook diffusive shock acceleration scenario, would produce a spectrum much steeper (f(p)p4.7similar-to𝑓𝑝superscript𝑝4.7f(p)\sim p^{-4.7}italic_f ( italic_p ) ∼ italic_p start_POSTSUPERSCRIPT - 4.7 end_POSTSUPERSCRIPT) than what is required to explain the gamma-ray observations (f(p)p4.2similar-to𝑓𝑝superscript𝑝4.2f(p)\sim p^{-4.2}italic_f ( italic_p ) ∼ italic_p start_POSTSUPERSCRIPT - 4.2 end_POSTSUPERSCRIPT according to Menchiari et al. 2024). On the other hand, it is not clear if there is at all such a large-scale collective termination shock around Cygnus OB2. It is usually believed that a collective wind can form if and only if the cluster WTS predicted by Weaver et al. (1977) extends beyond the extension of the cluster core. While this is surely a necessary condition, it may not be sufficient. In particular, if the distance between the O stars is larger than the individual stellar WTSs blown over a few Myr, one should not expect efficient wind-wind interactions, and the notion of a collective wind becomes ill-defined.

Beyond phenomenological arguments, detailed simulations are key to properly understand star-forming regions, in particular to probe the properties of the plasma, understand the kinetic structures and identify the sites of efficient particle acceleration. There have been many efforts in the recent years towards global hydrodynamic or magnetohydrodynamic simulations of stellar clusters forming superbubbles, often aiming at probing the feedback of the superbubble onto the interstellar medium at very large (100s pc) scales. For this purpose, the star clusters are often modelled as point-like injections of energy, and the impact of supernovae is studied in detail (e.g. Agertz et al., 2013; Krause et al., 2013; Kim et al., 2017; El-Badry et al., 2019; Lancaster et al., 2021). These studies provide information on the expansion of the superbubble shell, the impact of thermal conduction, radiative cooling, photoionisation, density fluctuations, the occurrence of shell instabilities, or the impact of the presence of clumps and turbulence in the vicinity of the cluster. This is relevant not only to understand how the stellar feedback carves the parent molecular cloud, but also to investigate the dynamics of the Galaxy at large scales, e.g., the distribution of hot gas and the regulation of star-formation (e.g. Gatto et al., 2017). However, in order to investigate structures at intermediate scales (10similar-toabsent10\sim 10∼ 10 pc), refined simulations must be set up, which is particularly relevant when the cluster wind and supernovae are blown inside an inhomogeneous, dense, and turbulent molecular cloud (e.g. Rogers & Pittard, 2013; Geen et al., 2021). In these simulations, the feedback is found to be highly asymmetric and deviates strongly from the theoretical picture, although the cluster is still modelled as a point-like source.

In order to understand the structure of the supersonic flow in the close vicinity of the cluster, and therefore the shape of the cluster wind termination shock, one needs to go beyond a point-like injection of energy. A possible improvement is to inject each individual star as a source of thermal energy (Gupta et al., 2018b), which effectively creates an extended, continuous, injection region. Although this might be a good approximation to model the feedback of a very compact and powerful star cluster – such as Westerlund 1 where several hundreds of OB stars are gathered within a core of 1 pc (Clark et al., 2005) – it is less relevant in the case of a more extended cluster like Cygnus OB2, where approximately 80 O stars are distributed over about 15 pc. Besides, as will be discussed in detail in Section 2, the wind power in Cygnus OB2 is dominated by a few very powerful O and Wolf-Rayet stars. In such a case, it is necessary to simulate each star individually using a resolution that is high enough to resolve the regions of wind-wind interactions at sub-parsec scales (Badmaev et al., 2022; Vieu et al., 2024). It is computationally challenging to resolve such small scales while keeping track of the medium-size structures (1020102010-2010 - 20 pc) inside a numerical box that is large enough to contain the large-scale superbubble (501005010050-10050 - 100 pc) over millions of years.

The goal of the present work is two-fold: Firstly, we aim at highlighting astrometric results which are often oversimplified by the high-energy astrophysics community. Secondly, we probe the presence of shocks inside the superbubble by running a high-resolution 3D hydrodynamic simulation of a cluster whose distribution of massive stars statistically matches that of Cygnus OB2.

Section 2 is devoted to a detailed description of the massive star population of Cygnus OB2 which is then used as a basis to set-up our simulation. Section 3 described the results of the simulation, highlighting in particular the absence of a large-scale cluster WTS around the cluster, and discussing consequences on scenarios of particle acceleration. We conclude in Section 4.

2 Simulating the massive star population in Cygnus OB2

The aim of this section is to model Cygnus OB2 as accurately as allowed by the available data. In order to get around uncertainties (e.g. in parallaxes measurements or in stellar properties), we will have to make assumptions and extrapolations. These will be systematically designed to maximise the mechanical feedback of the cluster, i.e. our model cluster will be likely more compact and more powerful than the real configuration of Cygnus OB2.

Refer to caption
Figure 1: Mechanical energy versus effective temperature for the considered OB and WR stars in Cygnus OB2. Orange dots are denoting the O stars with known wind strength parameters from Berlanas et al. (2020). Specific OB and WR stars are outlined in green or red respectively. For objects with effective temperatures outside the plot range, arrows are given instead of dots. The dashed orange line shows the integrated mechanical luminosity of all O stars with individually known wind strength parameters. The red dashed line shows the total mechanical luminosity including the three WR stars.

2.1 Stellar population from observations

Cygnus OB2 is one of the best examples of a young massive star cluster close to Earth, with a total stellar mass estimated around 16 5002800+3800subscriptsuperscript165003800280016\,500^{+3800}_{-2800}16 500 start_POSTSUPERSCRIPT + 3800 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT - 2800 end_POSTSUBSCRIPT M (Wright et al., 2015). Its massive star population has been continuously studied and extended over the years (e.g., Comerón et al., 2002; Comerón & Pasquali, 2012; Wright et al., 2015; Berlanas et al., 2018), although the high extinction in its direction (AV56similar-tosubscript𝐴𝑉56A_{V}\sim 5-6\,italic_A start_POSTSUBSCRIPT italic_V end_POSTSUBSCRIPT ∼ 5 - 6mag, Wright et al., 2015) hinders the observations. The OB population is still not considered to be complete, nevertheless this mainly limits the completeness with respect to faint, late-type O stars which are not expected to provide a significant contribution to the total mechanical luminosity of the cluster (Berlanas et al., 2018).

Cygnus OB2 is located in the northern sky, at a distance closer than 2 kpc (Rygl et al., 2012). The Gaia parallaxes reveal substructures along the line of sight (Berlanas et al., 2019; Orellana et al., 2021), namely a small complex in the foreground at similar-to\sim1.41.41.4\,1.4kpc and a main complex at similar-to\sim1.71.71.7\,1.7kpc. The latter is further divided into two groups of stars centred around the massive systems Cygnus OB2 #22 and Cygnus OB2 #8 (Bica et al., 2003), which are sometimes referred to as clusters themselves named “Bica 1” and “Bica 2” (e.g., Maíz Apellániz et al., 2020). Both systems contain several luminous, early-type O stars including one of the rare O3 spectral class in each of them, namely Cygnus OB2 #7 (O3f*) in Bica 2 and Cygnus OB2 #22 A (also O3f*) in Bica 1 (Sota et al., 2011). The two clusters form a double-cluster as they are located only similar-to\sim0.10.10.1\,0.1deg away from each other in the projected plane on the sky.

In the present work, we consider the area within 1deg21superscriptdegree21\deg^{2}1 roman_deg start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT centred on Cygnus OB2 #8 (i.e., the centre of the “Bica 2 cluster”) as in e.g. Wright et al. (2015). We assume a cluster distance of 1.651.651.65\,1.65kpc, motivated by the recent analysis from Maíz Apellániz et al. (2020), who obtained a Gaia-based distance of 1.640.11+0.13subscriptsuperscript1.640.130.111.64^{+0.13}_{-0.11}\,1.64 start_POSTSUPERSCRIPT + 0.13 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT - 0.11 end_POSTSUBSCRIPTkpc considering the Bica 2 group. For comparison, the derived distance by Maíz Apellániz et al. (2020) for Bica 1 (containing Cyg OB #22) is 1.720.12+0.14subscriptsuperscript1.720.140.121.72^{+0.14}_{-0.12}\,1.72 start_POSTSUPERSCRIPT + 0.14 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT - 0.12 end_POSTSUBSCRIPTkpc.

The age of Cygnus OB2 is ill-defined as the region is not coeval. Estimates range between 2±1plus-or-minus212\pm 12 ± 1 Myr for a population of mid-type O and B supergiant (BSG) stars (Hanson, 2003) up to 5–7 Myr for a population of A-type stars (Drew et al., 2008). Star-formation has likely happened over the last 6 Myr, with two main starburst events at approximately 3 and 5 Myr (Berlanas et al., 2020). It is believed that the current OB stellar population formed likely later than most of the observed lower-mass stars, which has also been concluded for other OB clusters (e.g., Povich & Whitney, 2010; Ramachandran et al., 2018; Schneider et al., 2018).

Three Wolf-Rayet (WR) stars can be found in the vicinity of Cygnus OB2. These stars provide a significant fraction of the total stellar wind power, and need to be considered in detail due to their positions far from the centre of the cluster. WR 144 is a WC4-type star with logM˙=4.62˙𝑀4.62\log\dot{M}=-4.62roman_log over˙ start_ARG italic_M end_ARG = - 4.62 and \varv=3500km s1subscript\varv3500superscriptkm s1\varv_{\infty}=3500\,~{}\textrm{km\,s}^{-1}start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT = 3500 km s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT (Sander et al., 2019). WR 145 is a WR+O7V binary with a transition type (WN7/WC) primary (Muntean et al., 2009). We cannot resolve the individual stars in the system and thus use the values of logM˙=4.49˙𝑀4.49\log\dot{M}=-4.49roman_log over˙ start_ARG italic_M end_ARG = - 4.49 and \varv=1440km s1subscript\varv1440superscriptkm s1\varv_{\infty}=1440\,~{}\textrm{km\,s}^{-1}start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT = 1440 km s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT from the single-star analysis by Sander et al. (2012); Sander et al. (2019) which provides a conservative estimate of the mass-loss rate. WR 146 is a WC6+O8I binary, with \varv=2900km s1subscript\varv2900superscriptkm s1\varv_{\infty}=2900\,~{}\textrm{km\,s}^{-1}start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT = 2900 km s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT (Eenens & Williams, 1994) for the WR component, though the mass-loss rate is uncertain, both for the individual components and for the combined spectrum. We therefore make an assumption based on the average value for Galactic WC stars from Table 5 in Sander et al. (2019) and scale this to the luminosity of WR 146, giving us logM˙4.9˙𝑀4.9\log\dot{M}\approx-4.9roman_log over˙ start_ARG italic_M end_ARG ≈ - 4.9 (see, e.g., Sander et al., 2019, for the corresponding equations). This estimate is compatible with recent analyses showing that the mass-loss rate of WR146 is likely smaller than previously thought (Zhekov, 2017; Pittard et al., 2021a).

In addition to the WR stars, we also specifically model Cygnus OB2 #12, a blue hypergiant of spectral type B3Ia+ and one of the most luminous star in the Milky Way (Clark et al., 2012). The star is possibly a Colliding Wind Binary (Oskinova et al., 2017), but similar to other binary systems the components cannot be treated individually. We therefore adopt \varv=400km s1subscript\varv400superscriptkm s1\varv_{\infty}=400\,~{}\textrm{km\,s}^{-1}start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT = 400 km s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT and logM˙=5.52˙𝑀5.52\log\dot{M}=-5.52roman_log over˙ start_ARG italic_M end_ARG = - 5.52 from Clark et al. (2012).

2.2 Simulated population of massive stars

The main data source for the stars in our simulation is the spectroscopic sample of 78 O stars from Berlanas et al. (2020), which marks the most complete spectral sample for the O-star population of Cygnus OB2 so far. Using a grid of FASTWIND (Puls et al., 2005; Rivero González et al., 2012) stellar atmosphere models and the iacob-gbat tool (Simón-Díaz et al., 2011; Holgado et al., 2018), Berlanas et al. (2020) have provided the effective temperature Teffsubscript𝑇effT_{\mathrm{eff}}italic_T start_POSTSUBSCRIPT roman_eff end_POSTSUBSCRIPT, the stellar radius Rsubscript𝑅R_{\star}italic_R start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT and the wind-strength parameter

Q𝑄\displaystyle Qitalic_Q =M˙M/yr(km s1\varvRR)3/2,absent˙𝑀subscript𝑀direct-productyrsuperscriptsuperscriptkm s1subscript\varvsubscript𝑅direct-productsubscript𝑅32\displaystyle=\frac{\dot{M}}{M_{\odot}/\mathrm{yr}}\left(\frac{\,~{}\textrm{km% \,s}^{-1}}{\varv_{\infty}}\frac{R_{\odot}}{R_{\star}}\right)^{3/2}\,,= divide start_ARG over˙ start_ARG italic_M end_ARG end_ARG start_ARG italic_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT / roman_yr end_ARG ( divide start_ARG km s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT end_ARG start_ARG start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT end_ARG divide start_ARG italic_R start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT end_ARG start_ARG italic_R start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT end_ARG ) start_POSTSUPERSCRIPT 3 / 2 end_POSTSUPERSCRIPT , (1)

for the 52 O stars in Cygnus OB2 that have reliable Gaia DR2 astrometry. Given the considerable reddening of Cygnus OB2, the analysis of Berlanas et al. (2020) is limited to the optical wavelength range. Without additional constraints, e.g., UV measurements of the terminal wind velocity, \varvsubscript\varv\varv_{\infty}start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT, the combined wind strength parameter Q𝑄Qitalic_Q can be robustly constrained from the spectral fit but not the individual values of \varvsubscript\varv\varv_{\infty}start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT and the mass-loss rate M˙˙𝑀\dot{M}over˙ start_ARG italic_M end_ARG. As we wish to know the mechanical luminosity

Lmech=12M˙\varv2,subscript𝐿mech12˙𝑀superscriptsubscript\varv2L_{\text{mech}}=\frac{1}{2}\dot{M}\varv_{\infty}^{2}\,,italic_L start_POSTSUBSCRIPT mech end_POSTSUBSCRIPT = divide start_ARG 1 end_ARG start_ARG 2 end_ARG over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT , (2)

for our simulation, we need an additional estimate for \varvsubscript\varv\varv_{\infty}start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT to calculate M˙˙𝑀\dot{M}over˙ start_ARG italic_M end_ARG via Eq. (1). In principle, we could estimate the terminal wind velocity from the (effective) escape velocity using the modified-CAK relation between them (Castor et al., 1975; Kudritzki & Puls, 2000), but given the considerable uncertainties of spectroscopic masses and the available quantities from Berlanas et al. (2020), we instead use the empirically well-constrained relation between \varvsubscript\varv\varv_{\infty}start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT and Teffsubscript𝑇effT_{\mathrm{eff}}italic_T start_POSTSUBSCRIPT roman_eff end_POSTSUBSCRIPT. The recent determination by Hawcroft et al. (2023) for O stars at Galactic metallicity yields the formula

\varvsubscript\varv\displaystyle\varv_{\infty}start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT =(0.102TeffK1300)km s1,absent0.102subscript𝑇effK1300superscriptkm s1\displaystyle=\left(0.102\frac{T_{\mathrm{eff}}}{\mathrm{K}}-1300\right)\,~{}% \textrm{km\,s}^{-1}\,,= ( 0.102 divide start_ARG italic_T start_POSTSUBSCRIPT roman_eff end_POSTSUBSCRIPT end_ARG start_ARG roman_K end_ARG - 1300 ) km s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT , (3)

which we use to estimate values of \varvsubscript\varv\varv_{\infty}start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT and thus obtain both M˙˙𝑀\dot{M}over˙ start_ARG italic_M end_ARG and Lmechsubscript𝐿mechL_{\text{mech}}italic_L start_POSTSUBSCRIPT mech end_POSTSUBSCRIPT, avoiding any explicit assumptions about the spectroscopic mass. We note that the values for Q𝑄Qitalic_Q derived by Berlanas et al. (2020) and also the inferred M˙˙𝑀\dot{M}over˙ start_ARG italic_M end_ARG values assume smooth winds and thus the real mass-loss rates could even be a bit lower (see Puls et al., 2008; Hamann et al., 2008; Oskinova et al., 2016, for recent reviews). Moreover, Eq. (3) only reproduces a statistical trend with a non-negligible spread, meaning that the values for individual objects could be off (see Table 1), but given the aims of our simulations, these are expected to average out. In any case, as we will discuss further below, we do not aim at an exact reproduction of Cygnus OB2, but instead a statistically equivalent cluster.

There remain 26 O stars for which the stellar radius or wind parameters are not known. In order to complete our sample, we generate a synthetic stellar population by exactly sampling the initial mass function of Cygnus OB2, with a slope of Γ=1.39Γ1.39\Gamma=1.39roman_Γ = 1.39 (Wright et al., 2015). We pick stars with masses between 20 and 50 Msubscript𝑀direct-product\,M_{\odot}italic_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT in order to avoid introducing unusually powerful stars which are not seen in the observations. On the other hand, starting at 20 Msubscript𝑀direct-product\,M_{\odot}italic_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT is expected to produce more powerful stars than in reality, since the remaining sample is likely biased toward lower masses. Overestimating the power provides an optimistic case for the formation of a cluster WTS around Cygnus OB2.

To infer plausible values for the wind parameters from the initial mass M𝑀Mitalic_M, we use the prescriptions of Seo et al. (2018) (in the case of non-rotating single stars):

log10M˙Myr1=3.38(log10MM)2+14.59log10MM20.84,subscript10˙𝑀subscript𝑀direct-productsuperscriptyr13.38superscriptsubscript10𝑀subscript𝑀direct-product214.59subscript10𝑀subscript𝑀direct-product20.84\displaystyle\log_{10}\frac{\dot{M}}{{M_{\odot}\,\mathrm{yr}^{-1}}}=-3.38\left% (\log_{10}\frac{M}{M_{\odot}}\right)^{2}+14.59\log_{10}\frac{M}{M_{\odot}}-20.% 84\,,roman_log start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT divide start_ARG over˙ start_ARG italic_M end_ARG end_ARG start_ARG italic_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT roman_yr start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT end_ARG = - 3.38 ( roman_log start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT divide start_ARG italic_M end_ARG start_ARG italic_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT end_ARG ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + 14.59 roman_log start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT divide start_ARG italic_M end_ARG start_ARG italic_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT end_ARG - 20.84 , (4)
log10\varvkm s1=0.08log10MM+3.28.subscript10subscript\varvsuperscriptkm s10.08subscript10𝑀subscript𝑀direct-product3.28\displaystyle\log_{10}\frac{\varv_{\infty}}{\textrm{km\,s}^{-1}}=0.08\log_{10}% \frac{M}{M_{\odot}}+3.28\,.roman_log start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT divide start_ARG start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT end_ARG start_ARG km s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT end_ARG = 0.08 roman_log start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT divide start_ARG italic_M end_ARG start_ARG italic_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT end_ARG + 3.28 . (5)

In the end we obtain a population of O stars with a mechanical power of 1038superscript103810^{38}10 start_POSTSUPERSCRIPT 38 end_POSTSUPERSCRIPT erg/s. The mechanical luminosities of this OB stellar population is shown in Fig. 1, including their total sum. When adding the three nearby WR stars, as well as the Blue Hypergiant Cygnus OB2 #12, it becomes evident that the latter, despite its high luminosity, has very little impact on the overall budget. This is due to the comparably low terminal velocity of the B3-4 hypergiant of only 400400400\,400km/s. When including the WR stars (cf. Fig. 1), we obtain a total mechanical power of 2.6×10382.6superscript10382.6\times 10^{38}2.6 × 10 start_POSTSUPERSCRIPT 38 end_POSTSUPERSCRIPT erg/s for the Cygnus OB2 cluster. The most powerful stars are listed in Table 1. Note that the whole O star population represents only about 40% of the total mechanical power in the cluster with the remaining 60% being generated by the three WR stars. This underlines the importance to properly identify WR stars in stellar populations as they can completely change the mechanical energy budget.

We also show estimated uncertainties in Table 1. These are derived by standard error propagation, assuming conservative standard values for ϵ(logM˙)=0.15 [dex]italic-ϵ˙𝑀0.15 [dex]\epsilon(\log\dot{M})=0.15\text{ [dex]}italic_ϵ ( roman_log over˙ start_ARG italic_M end_ARG ) = 0.15 [dex] and ϵ(\varv)=200km s1italic-ϵsubscript\varv200superscriptkm s1\epsilon(\varv_{\infty})=200\,~{}\textrm{km\,s}^{-1}italic_ϵ ( start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT ) = 200 km s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT.

WR144 WR145 WR146 Cygnus OB2 #7 Cygnus OB2 Cygnus OB2 #11 Cygnus OB2 #22 O stars O stars
#8B+#8C (total) (simulated)
Lwsubscript𝐿𝑤L_{w}italic_L start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT 9.3±2.7plus-or-minus9.32.79.3\pm 2.79.3 ± 2.7 2.1±0.95plus-or-minus2.10.952.1\pm 0.952.1 ± 0.95 3.4±1.0plus-or-minus3.41.03.4\pm 1.03.4 ± 1.0 2.4±0.71plus-or-minus2.40.712.4\pm 0.712.4 ± 0.71 0.99±0.32plus-or-minus0.990.320.99\pm 0.320.99 ± 0.32 1.2±0.39plus-or-minus1.20.391.2\pm 0.391.2 ± 0.39 1.1±0.32plus-or-minus1.10.321.1\pm 0.321.1 ± 0.32 10101010 7.77.77.77.7
Table 1: Distribution of the stellar power in Cygnus OB2, in units of 1037superscript103710^{37}10 start_POSTSUPERSCRIPT 37 end_POSTSUPERSCRIPT erg/s. Most of the feedback is blown by the WR stars and 5 very powerful O stars. The discrepancy between the total sample and the simulated sample is due to the stars located beyond 12 pc being discarded because of numerical limitations, as discussed in Section 2.3.

Although the simulated stars are not individually identical to that of Cygnus OB2, our procedure allows to generate a synthetic cluster that is statistically equivalent to the real distribution. In any case, the precise power of a standard O star is not expected to strongly affect the simulation results.

2.3 3D positioning

Refer to caption
Figure 2: Statistical properties of our sample. The black dashed line at 12 pc in each panel shows where we cut off our simulated population. The upper panel shows log10Qsubscript10𝑄\log_{10}Qroman_log start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT italic_Q as a function of projected radial distance Rprojsubscript𝑅projR_{\mathrm{proj}}italic_R start_POSTSUBSCRIPT roman_proj end_POSTSUBSCRIPT for the individual stars where we have Q𝑄Qitalic_Q. The lower panel is a histogram showing the number of O stars with Rsubscript𝑅R_{\star}italic_R start_POSTSUBSCRIPT ⋆ end_POSTSUBSCRIPT and Q𝑄Qitalic_Q available.

The most important factor that influences the feedback is the spatial distribution of the stars. Indeed, one should expect qualitatively different results if all the stars are gathered in a compact core, or rather loosely scattered in a large volume. The radial distribution of Cygnus OB2 as catalogued in Berlanas et al. (2020) is shown in Fig. 2. The histogram in the bottom panel shows that Cygnus OB2 is actually neither very compact nor very loose. Although most of the stars and most of the wind power are concentrated in the inner few parsecs, the distribution does extend up to large distances. In particular, the WR stars are notably off-centred, at projected distances of 6.8 pc (WR144), 13.8 pc (WR146) and 15.8 pc (WR145). In fact, most of the cluster mechanical power is blown by a few stars, as seen in Fig. 1 and summarised in Table 1: the three WR stars, the very powerful O star Cygnus OB2 #7 , and eventually the stars Cygnus OB2 #8B and Cygnus OB2 #8C, which are in such close vicinity that they can be considered as a single star from the point of view of their feedback. It is clear that such a scattered distribution of stars cannot be approximated by a continuous distribution. Massive stars need to be resolved individually in order to properly encapsulate the wind-wind interactions and a 3D positioning is required for all the simulated stars.

Unfortunately, reliable Gaia parallaxes only exist for just over half of the sample (Berlanas et al., 2020). We therefore need to extrapolate from the right ascension and declination in order to determine a position in 3D. To produce a list of 3D coordinates for each of our stars, we first define a pair of (x,y)𝑥𝑦(x,y)( italic_x , italic_y ) coordinates with respect to the centre of Cygnus OB2, which are calculated using Gaia right ascension and declination assuming a distance of 1.65 kpc to the centre of the cluster. We then extrapolate a 3D distance d=x2+y2+((|x|+|y|)/2)2𝑑superscript𝑥2superscript𝑦2superscript𝑥𝑦22d=\sqrt{x^{2}+y^{2}+\left((\lvert x\rvert+\lvert y\rvert)/2\right)^{2}}italic_d = square-root start_ARG italic_x start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_y start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + ( ( | italic_x | + | italic_y | ) / 2 ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG. This procedure provides a conservative estimate of the real distance, which could actually be much larger in 3D. Finally, in order to avoid a biased 3D reconstruction, in particular to produce a distribution that is statistically spherically symmetric, we randomly pick a point on the sphere located at the distance d𝑑ditalic_d from the cluster centre. By construction, this reproduces the radial distribution of stars in the cluster, while allowing us to make sensible estimates for the z-component. Fig. 3 shows that the projected distribution of stars in the projected plane resembles qualitatively that given in Wright et al. (2015).

Refer to caption
Figure 3: Comparison of the absolute x𝑥xitalic_x and y𝑦yitalic_y separations for the “True” values given in Berlanas et al. (2020) (blue markers, assuming a distance to the cluster of 1.65 kpc) compared with the separations of the 36 O stars in our simulation (orange). We also show the “True” positions of the WR stars with star markers in blue, and their shifted positions in the simulation in orange. A 3D visualisation of the simulated cluster is shown in Fig. 12.

To accurately resolve massive stars individually in a simulation requires a high resolution. Extending this resolution over the entire domain would be computationally too expensive. Therefore we need to define a “core region”, at |x|,|y|,|z|<Δc𝑥𝑦𝑧subscriptΔ𝑐|x|,|y|,|z|<\Delta_{c}| italic_x | , | italic_y | , | italic_z | < roman_Δ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT, such that stars located outside of the cube of half edge-length ΔcsubscriptΔ𝑐\Delta_{c}roman_Δ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT will be excluded from the simulation. Our computational resources allow us to set Δc=12subscriptΔ𝑐12\Delta_{c}=12roman_Δ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT = 12 pc. As seen in Fig. 2, this is enough to include most of the powerful stars, with the notable exceptions of WR145 and WR146, which, with a mechanical power of respectively 2.1×10372.1superscript10372.1\times 10^{37}2.1 × 10 start_POSTSUPERSCRIPT 37 end_POSTSUPERSCRIPT erg/s and 3.4×10373.4superscript10373.4\times 10^{37}3.4 × 10 start_POSTSUPERSCRIPT 37 end_POSTSUPERSCRIPT erg/s, are the second and third most powerful stars in Cygnus OB2, thus should not be discarded. We artificially bring them closer, shifting the y𝑦yitalic_y coordinate of WR145 from 14.5 pc to 11.9 pc, and the x𝑥xitalic_x coordinate of WR146 from 13.6 pc to 10.6 pc. For these two stars, we set z=0𝑧0z=0italic_z = 0 without loss of generality.

Stars which are not powerful enough to expand a supersonic wind against the ISM pressure beyond the injection cells cannot be included in the simulation (Pittard et al., 2021b): we must discard stars such that M˙\varv<(2000km/s)×(107M/yr)˙𝑀subscript\varv2000km/ssuperscript107subscript𝑀direct-product/yr\dot{M}\varv_{\infty}<(2000~{}\textrm{km/s})\times(10^{-7}~{}M_{\odot}\textrm{% /yr})over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT < ( 2000 km/s ) × ( 10 start_POSTSUPERSCRIPT - 7 end_POSTSUPERSCRIPT italic_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT /yr ). This represents a population of 20 O stars which has a negligible mechanical power (1.8×10361.8superscript10361.8\times 10^{36}1.8 × 10 start_POSTSUPERSCRIPT 36 end_POSTSUPERSCRIPT erg/s) and is therefore not expected to affect the cluster feedback at large scales.

Removing all O stars located outside of the cube of half edge-length 12 pc, we end up with a population of 36 O stars, which has a total mechanical power of 7.7×10377.7superscript10377.7\times 10^{37}7.7 × 10 start_POSTSUPERSCRIPT 37 end_POSTSUPERSCRIPT erg/s, 73% of which being blown by the most powerful O stars #7, #8, #11, #22. Cygnus OB2-8B and Cygnus OB2-8C are located too close to each other (0.12 pc) to be individually resolved in the inner parsec, so we merge them together by adding up their mass-loss rate and mechanical power.

Fig. 4 shows the velocity dispersion of the stars in the Berlanas et al. (2020) sample. We see that the vast majority of the stars are moving at no more than a few km s1superscriptkm s1\,~{}\textrm{km\,s}^{-1}km s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT, corresponding to a proper motion of order 4pc/Myrsimilar-toabsent4pcMyr\sim 4\,\mathrm{pc}/\mathrm{Myr}∼ 4 roman_pc / roman_Myr. Since it typically takes a few hundreds of kyrs to reach a quasi-stationary state in the simulation, it is reasonable to neglect the proper motion of individual stars. Besides, as discussed in detail in Wright et al. (2016), observations show no hint of a radial expansion of the cluster: it is unlikely that Cygnus OB2 was more compact in the past.

Finally we set up an additional O star near the centre of the cluster with fiducial wind parameters M˙=106M˙𝑀superscript106subscript𝑀direct-product\dot{M}=10^{-6}M_{\odot}over˙ start_ARG italic_M end_ARG = 10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT italic_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT/yr, \varv=2500subscript\varv2500\varv_{\infty}=2500start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT = 2500 km/s. This star is introduced as a putative supernova progenitor, which will be used in a follow-up work. In the present simulation, it is kept in the main-sequence and has a negligible impact on the results.

Refer to caption
Figure 4: Histogram of individual stars’ projected proper motion Vpropsubscript𝑉propV_{\text{prop}}italic_V start_POSTSUBSCRIPT prop end_POSTSUBSCRIPT minus the mean proper motion of the sample V¯propsubscript¯𝑉prop\bar{V}_{\text{prop}}over¯ start_ARG italic_V end_ARG start_POSTSUBSCRIPT prop end_POSTSUBSCRIPT based on Gaia DR3 (Gaia Collaboration et al., 2023).

The stars forming our simulated cluster are listed in Appendix A, and Fig. 12 displays the reconstructed 3D distribution. Although this population is not exactly that of Cygnus OB2, it is representative of the real statistical distribution of the wind mechanical power in 3D space and, as discussed earlier, was designed to maximise the feedback when the parameters were uncertain. This synthetic population can only be used to draw qualitative conclusions on the gas density and flow configuration across a few tens of parsecs. Our aim is not to perform a quantitative comparison between the simulation results and specific observations such as X-ray data or gas density maps, but rather to probe the presence of a large-scale cluster wind termination shock, which is key to understanding the gamma-ray emission from the region.

Refer to caption
Figure 5: Two different 2D slices at x=0𝑥0x=0italic_x = 0 (left) and y=0𝑦0y=0italic_y = 0 (right) showing the density map at 3 different times. The purple outlines highlight the Mach=1 contours.

2.4 Initial conditions and stellar wind injection

The 3D simulation domain extends in the three Cartesian axes from -70 pc to +70 pc. This is large enough to contain the superbubble blown by the cluster during about 2 Myr. We set up 7763 cells over the core region defined by -12.5 pc <x,y,z<12.5formulae-sequenceabsent𝑥𝑦𝑧12.5<x,y,z<12.5< italic_x , italic_y , italic_z < 12.5 pc, therefore using a resolution of 0.032 pc/cell in this region. The grid is stretched beyond 12.5 pc up to the boundaries of the domain, degrading progressively the resolution, eventually using a total number of 10003 cells.

The initial condition is an ideal background gas of number density 20 cm-3, mean molecular weight 0.61 and temperature 104superscript10410^{4}10 start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT K. These values are compatible with recent observational analyses of the molecular clouds in the region (e.g. Zhang et al., 2024; Astiasarain et al., 2023; LHAASO Collaboration, 2024), with typical HI column density measurements on the order of 1022superscript102210^{22}~{}10 start_POSTSUPERSCRIPT 22 end_POSTSUPERSCRIPTcm-2. Given the complex structure of the region, it is difficult to accurately isolate the contribution of the molecular gas near Cygnus OB2. We choose a density typical of a parent molecular cloud which is expected to form a massive star cluster.

Setting a somewhat high external density is expected to inhibit the expansion of the superbubble blown by the massive stars, which is necessary to maintain the entire structure inside the numerical box over the simulation time (2 Myr). On the other hand, one should ensure that the pressure inside the superbubble remains larger than the external pressure over the simulation time, which is always true with our assumed parameters (at t=2𝑡2t=2italic_t = 2 Myr the internal pressure is still more than twice the external pressure).

Within the high-resolution core, each massive star is individually modelled as a spherically symmetric stellar wind of terminal velocity \varv,isubscript\varv𝑖\varv_{\infty,i}start_POSTSUBSCRIPT ∞ , italic_i end_POSTSUBSCRIPT and mass-loss rate M˙isubscript˙𝑀𝑖\dot{M}_{i}over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT:

𝒖=\varv,i𝒆𝒓,𝒊,ρ=M˙i4πr2\varv,i,formulae-sequence𝒖subscript\varv𝑖subscript𝒆𝒓𝒊𝜌subscript˙𝑀𝑖4𝜋superscript𝑟2subscript\varv𝑖\displaystyle\boldsymbol{u}=\varv_{\infty,i}\boldsymbol{e_{r,i}}\,,\qquad\rho=% \frac{\dot{M}_{i}}{4\pi r^{2}\varv_{\infty,i}}\,,bold_italic_u = start_POSTSUBSCRIPT ∞ , italic_i end_POSTSUBSCRIPT bold_italic_e start_POSTSUBSCRIPT bold_italic_r bold_, bold_italic_i end_POSTSUBSCRIPT , italic_ρ = divide start_ARG over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_ARG start_ARG 4 italic_π italic_r start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT ∞ , italic_i end_POSTSUBSCRIPT end_ARG , (6)

where the parameters \varv,isubscript\varv𝑖\varv_{\infty,i}start_POSTSUBSCRIPT ∞ , italic_i end_POSTSUBSCRIPT and M˙isubscript˙𝑀𝑖\dot{M}_{i}over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT are summarised for each star in Table 2 and are kept constant in time. The pressure in the wind is prescribed assuming an isothermal wind: pgas=ρcs2subscript𝑝gas𝜌superscriptsubscript𝑐𝑠2p_{\rm gas}=\rho c_{s}^{2}italic_p start_POSTSUBSCRIPT roman_gas end_POSTSUBSCRIPT = italic_ρ italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT, with cs=20subscript𝑐𝑠20c_{s}=20italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT = 20 km/s for O stars and cs=30subscript𝑐𝑠30c_{s}=30italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT = 30 km/s for WR stars. During the simulation, the stellar winds are set up within spheres of radius RH=0.16subscript𝑅𝐻0.16R_{H}=0.16italic_R start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT = 0.16 pc (5 cells) around each star, which is almost five times smaller than the distance between the two closest stars, allowing to properly resolve the regions of stellar wind collisions (except that between Cygnus OB2 #8B and Cygnus OB2 #8C which have been merged). This also ensures that for all simulated stars, the injection radius is significantly smaller than rinj,max=(M˙\varv/(4πP0))1/2subscript𝑟injmaxsuperscript˙𝑀subscript\varv4𝜋subscript𝑃012r_{\rm inj,max}=\left(\dot{M}\varv_{\infty}/(4\pi P_{0})\right)^{1/2}italic_r start_POSTSUBSCRIPT roman_inj , roman_max end_POSTSUBSCRIPT = ( over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT / ( 4 italic_π italic_P start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) ) start_POSTSUPERSCRIPT 1 / 2 end_POSTSUPERSCRIPT, with P0subscript𝑃0P_{0}italic_P start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT the initial pressure. This condition is necessary to accurately expand wind-blown bubbles (Pittard et al., 2021b). For the stars with the lowest M˙\varv˙𝑀subscript\varv\dot{M}\varv_{\infty}over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT product, we get RH/rinj,maxsubscript𝑅𝐻subscript𝑟injmaxR_{H}/r_{\rm inj,max}italic_R start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT / italic_r start_POSTSUBSCRIPT roman_inj , roman_max end_POSTSUBSCRIPT = 0.3, while for the most powerful stars, which dominate the stellar feedback at large scales, we have RH/rinj,max<0.1subscript𝑅𝐻subscript𝑟injmax0.1R_{H}/r_{\rm inj,max}<0.1italic_R start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT / italic_r start_POSTSUBSCRIPT roman_inj , roman_max end_POSTSUBSCRIPT < 0.1.

The initial (t=0𝑡0t=0italic_t = 0) injection setup is slightly different. The stellar winds are defined over 0.39 pc (12 cells), which corresponds to half the distance between the two closest stars. In this extended injection region, we can define a smoother cut-off for the velocity and the density: 𝒖=\varv,i/(1+(rRH)/RH)2𝒆𝒓,𝒊𝒖subscript\varv𝑖superscript1𝑟subscript𝑅𝐻subscript𝑅𝐻2subscript𝒆𝒓𝒊\boldsymbol{u}=\varv_{\infty,i}/\left(1+(r-R_{H})/R_{H}\right)^{2}\boldsymbol{% e_{r,i}}bold_italic_u = start_POSTSUBSCRIPT ∞ , italic_i end_POSTSUBSCRIPT / ( 1 + ( italic_r - italic_R start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT ) / italic_R start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT bold_italic_e start_POSTSUBSCRIPT bold_italic_r bold_, bold_italic_i end_POSTSUBSCRIPT, ρ=M˙i/(4πr2𝒖)𝜌subscript˙𝑀𝑖4𝜋superscript𝑟2norm𝒖\rho=\dot{M}_{i}/(4\pi r^{2}\|\boldsymbol{u}\|)italic_ρ = over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT / ( 4 italic_π italic_r start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ∥ bold_italic_u ∥ ). This was deemed necessary in order to obtain spherical stellar winds with the Cartesian grid.

WR stars are only expected to appear at the end of the main-sequence phase. WR progenitors are assumed to be O stars with fiducial wind parameters M˙=106M˙𝑀superscript106subscript𝑀direct-product\dot{M}=10^{-6}M_{\odot}over˙ start_ARG italic_M end_ARG = 10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT italic_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT/yr and \varv=2500subscript\varv2500\varv_{\infty}=2500start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT = 2500 km/s. Since the WR stars are considerably off-centred, their feedback in the main-sequence is not expected to play a major role. In the simulation, the transition to the WR phase happens at t=1.6𝑡1.6t=1.6italic_t = 1.6 Myr. We inject the WR winds over 6 cells using Eq. 6.

The Euler equations of hydrodynamics are solved using the publicly available code PLUTO (Mignone et al., 2007) with a Lax-Friedrichs scheme (TVDLF). Thermal conduction and radiative cooling are neglected in order to minimise the computational overhead in the high-resolution core. Besides, a proper treatment of thermal conduction cannot be done without accounting for the magnetic field. Thermal conduction and cooling are known to impact on the properties and stability of the superbubble shell, to increase the density inside the cavity due to the evaporation at the shell interface and to lower the temperature inside the cavity (Weaver et al., 1977; El-Badry et al., 2019). On the one hand, not encapsulating these effects implies that our simulation is not expected to provide reliable results for the dynamics of the superbubble shell, nor to provide quantitative predictions for the density and temperature, which could be confronted to radio, optical or X-ray observations. On the other hand, neither thermal conduction nor cooling are expected to substantially change the pressure inside the superbubble, which is the main driver for the formation of the cluster WTS. For the purpose of this work, it is therefore reasonable to neglect these processes.

3 Results

3.1 Structure of the superbubble

Snapshots, using 2D slices through the simulation at 150 kyr (early stage), 1 Myr and 1.59 Myr (just before the onset of WR stars) are shown in Fig. 5 while averaged radial profiles of density and velocity are plotted in Fig. 6. At the beginning of the simulation, each massive star is surrounded by its own wind-blown bubble. The latter expand over several tens of kyrs until they start to percolate and eventually merge to form an approximately spherical superbubble. The matter excavated out accumulates in the superbubble shell, whose position matches the theoretical expectation (see the top panel of Fig. 6), although the shell width is much broader than theoretically expected due to the absence of radiative cooling in the simulation. The density inside the superbubble is low, of the order of 0.01 cm-3, although accounting for thermal conduction would trigger mass loading from the shell and increase this value.

Refer to caption
Figure 6: Radial profiles of number density (top panel) and velocity (bottom panel) at 0.2 Myr (early phase), 1 Myr (quasi-stationary state), 1.59 Myr (just before the onset of WR stars) and 2 Myr (end of the simulation). The vertical lines in the top panel show the theoretical position of the superbubble shell according to Weaver et al. (1977).
Refer to caption
Figure 7: Density maps at 2 Myr (400 kyr after the onset of WR stars) in two different slices at y=0𝑦0y=0italic_y = 0 (top) and z=0𝑧0z=0italic_z = 0 (bottom). Left panels: full view of the superbubble. Right panels: zoomed-in views. The red outlines show the Mach=1 contours.

From Fig. 6 we see that the radial profiles averaged over the azimuthal angles are relatively smooth, in contrast with theoretical expectation (see e.g. Fig. 1 of Gupta et al., 2018a). This is understood since we do not have a spherical wind around a compact central cluster. The winds blown by the powerful O stars in the inner few parsecs (namely Cygnus OB2 #7, Cygnus OB2 #8B+#8C and Cygnus OB2 #22) are strongly asymmetric. The supersonic outflows are shaped, through wind-wind collisions, into two-dimensional sheets, and then possibly, via subsequent collisions with other stellar winds, to one-dimensional jets that extend over large distances (up to 10 pc).

Overall, the level of wind-wind interaction is very low. At this stage, most of the stars are still surrounded by a small-scale (1absent1\approx 1≈ 1 pc) stellar wind termination shock, which slowly expands in the low-density cavity.

After 1 Myr, the simulation has already reached a quasi-stationary state, in the sense that we do not witness qualitative changes anymore. Due to the high external pressure, the forward shock of the superbubble has already become subsonic. The expansion of the superbubble follows the expected t3/5superscript𝑡35t^{3/5}italic_t start_POSTSUPERSCRIPT 3 / 5 end_POSTSUPERSCRIPT scaling (Weaver et al., 1977). Individual wind termination shocks, trans-sonic sheets and jets continue to expand slowly, in an almost self-similar manner.

Wolf-Rayet stars are introduced at 1.6 Myr and blow powerful winds until the end of the simulation. Fig. 7 shows the final snapshot at 2 Myr. At this point, the Wolf-Rayet stars are expected to either explode in supernovae or collapse into black-holes. The presence of WR stars has overall little impact on the dynamics of the superbubble. This is not surprinsing since the size of the superbubble shell theoretically scales as Lc1/5superscriptsubscript𝐿𝑐15L_{c}^{1/5}italic_L start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 1 / 5 end_POSTSUPERSCRIPT, where Lcsubscript𝐿𝑐L_{c}italic_L start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT is the total wind power of the cluster. Since the three WR stars are strongly off-centred, their feedback actually acts against the expansion of the wind termination shocks from the central region. Indeed, their onset provokes the shrinking of the individual O star WTSs. Finally, 400 kyr after they have appeared, they are still far away from interacting despite having blown their own individual WTSs over a couple of parsecs.

Refer to caption
Figure 8: 3D views of the cluster core at 2 Myr, showing the isosurfaces at Mach=1. The colorscale shows the density on these surfaces. See also the corresponding movie in the online Supplementary Material.

Fig. 8 provides a 3D view of the cluster at 2 Myr showing the isosurfaces of Mach number equal to 1 (see also the corresponding movie in the online Supplementary Material). The spheres correspond to individual stellar WTS around single stars. It is clearly seen that only a handful of powerful stars dominate the mechanical feedback. Most O stars are not powerful enough to expand a stellar WTS beyond one parsec: there are seen as isolated small spheres on the 3D views. These views provide a better understanding of the flow configuration and highlight in particular the absence of any large-scale spherical shock in the system, apart from individual WR shocks which extend over a couple of parsecs, as theoretically expected. Wind-wind collisions are clearly seen in particular between WR stars and O stars. These collisions produce trans-sonic planes, possibly curved inward the less powerful star. Further interaction with multiple stars sculpt the trans-sonic flows. We witness in particular the formation of a jet in the (x,y)𝑥𝑦(x,y)( italic_x , italic_y ) plane, as can also be seen in lower panels of Fig. 7.

3.2 Understanding the absence of a large-scale cluster wind termination shock

It is usually assumed that if the putative cluster WTS is larger than the cluster core, then it will actually form beyond the boundary of the core (e.g. Vieu & Reville, 2023; Menchiari et al., 2024). It appears that this argument is too naive. Firstly, a high number of massive stars organised in several layers is required to trap the kinetic energy of the winds within a rather small region of energy deposition, and eventually produce a collective, radial, outflow. If the number of stars is too small, the hydrodynamics is dominated by a few wind-wind collisions which create highly asymmetric structures. Secondly, the stars must be close enough to each other in order for their winds to interact. Indeed, the kinetic energy drops rapidly (typically as r4superscript𝑟4r^{-4}italic_r start_POSTSUPERSCRIPT - 4 end_POSTSUPERSCRIPT) beyond the WTS of a single star, which implies that the kinetic energy injected by a single star is contained within its own WTS. The latter expands in the shocked plasma which is nearly statistically homogeneous inside the superbubble. In order to enable efficient wind-wind interactions, the individual wind termination shocks must be close enough to one another, something which is not the case if the average distance between the stars greatly exceeds the typical size of a WTS blown by an O star. The textbook scenario of a continuous extended deposition of energy (Chevalier & Clegg, 1985) actually describes the extreme case where the wind-wind interactions are so efficient that all the stellar wind kinetic energy is converted to thermal pressure within the core. Although this might apply in some extremely compact and powerful clusters such as Westerlund 1 (e.g. Härer et al., 2023), it is not the case for Cygnus OB2.

The average pressure inside the superbubble (in the shocked plasma downstream of individual wind termination shocks) scales as follows (Weaver et al., 1977):

PSB0.16Lc2/5ρ03/5t4/5,subscript𝑃𝑆𝐵0.16superscriptsubscript𝐿𝑐25superscriptsubscript𝜌035superscript𝑡45P_{SB}\approx 0.16L_{c}^{2/5}\rho_{0}^{3/5}t^{-4/5}\,,italic_P start_POSTSUBSCRIPT italic_S italic_B end_POSTSUBSCRIPT ≈ 0.16 italic_L start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 / 5 end_POSTSUPERSCRIPT italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 / 5 end_POSTSUPERSCRIPT italic_t start_POSTSUPERSCRIPT - 4 / 5 end_POSTSUPERSCRIPT , (7)

where Lcsubscript𝐿𝑐L_{c}italic_L start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT is the total power of the cluster, ρ0subscript𝜌0\rho_{0}italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT is the initial density and t𝑡titalic_t the age of the cluster. This relation is derived by imposing energy and momentum conservation in the entire superbubble and should therefore hold independently of the spatial distribution of the stars inside the superbubble. In order to estimate the position of the stellar WTS of a single star isolated in the superbubble, one has to balance the ram pressure of the stellar wind, Pw=M˙\varv/(4πRwts2)subscript𝑃𝑤˙𝑀subscript\varv4𝜋superscriptsubscript𝑅wts2P_{w}=\dot{M}\varv_{\infty}/(4\pi R_{\text{wts}}^{2})italic_P start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT = over˙ start_ARG italic_M end_ARG start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT / ( 4 italic_π italic_R start_POSTSUBSCRIPT wts end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) with the pressure inside the superbubble. This provides:

Rwts1.5(M˙106M/yr)1/2(\varv2500km/s)1/2×(Lc2×1038erg/s)1/5(t5Myr)2/5(ρ020mpcm3)3/10pc,subscript𝑅wts1.5superscript˙𝑀superscript106subscript𝑀direct-productyr12superscriptsubscript\varv2500kms12superscriptsubscript𝐿𝑐2superscript1038ergs15superscript𝑡5Myr25superscriptsubscript𝜌020subscript𝑚𝑝superscriptcm3310pcR_{\text{wts}}\approx 1.5\left(\frac{\dot{M}}{10^{-6}M_{\odot}\rm/yr}\right)^{% 1/2}\left(\frac{\varv_{\infty}}{2500~{}\rm km/s}\right)^{1/2}\\ \times\left(\frac{L_{c}}{2\times 10^{38}{\rm~{}erg/s}}\right)^{-1/5}\left(% \frac{t}{5~{}\rm Myr}\right)^{2/5}\left(\frac{\rho_{0}}{20m_{p}~{}{\rm cm}^{-3% }}\right)^{-3/10}\rm pc\,,start_ROW start_CELL italic_R start_POSTSUBSCRIPT wts end_POSTSUBSCRIPT ≈ 1.5 ( divide start_ARG over˙ start_ARG italic_M end_ARG end_ARG start_ARG 10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT italic_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT / roman_yr end_ARG ) start_POSTSUPERSCRIPT 1 / 2 end_POSTSUPERSCRIPT ( divide start_ARG start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT end_ARG start_ARG 2500 roman_km / roman_s end_ARG ) start_POSTSUPERSCRIPT 1 / 2 end_POSTSUPERSCRIPT end_CELL end_ROW start_ROW start_CELL × ( divide start_ARG italic_L start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT end_ARG start_ARG 2 × 10 start_POSTSUPERSCRIPT 38 end_POSTSUPERSCRIPT roman_erg / roman_s end_ARG ) start_POSTSUPERSCRIPT - 1 / 5 end_POSTSUPERSCRIPT ( divide start_ARG italic_t end_ARG start_ARG 5 roman_Myr end_ARG ) start_POSTSUPERSCRIPT 2 / 5 end_POSTSUPERSCRIPT ( divide start_ARG italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG 20 italic_m start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT roman_cm start_POSTSUPERSCRIPT - 3 end_POSTSUPERSCRIPT end_ARG ) start_POSTSUPERSCRIPT - 3 / 10 end_POSTSUPERSCRIPT roman_pc , end_CELL end_ROW (8)

where M˙˙𝑀\dot{M}over˙ start_ARG italic_M end_ARG is the mass-loss rate of the isolated star and \varvsubscript\varv\varv_{\infty}start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT its wind velocity. We have normalised the parameters to values typical of Cygnus OB2. In fact, Rwtssubscript𝑅wtsR_{\text{wts}}italic_R start_POSTSUBSCRIPT wts end_POSTSUBSCRIPT depends only weakly on the parameters, which means it is quite general to state that a typical O star isolated in any superbubble will blow a supersonic wind over at most a few parsecs. In order for stellar winds to interact, the typical distance between the O stars must therefore be less than a few parsecs.

Refer to caption
Figure 9: Average distance between the stars within a given radius for stars with M˙>107M˙𝑀superscript107subscript𝑀direct-product\dot{M}>10^{-7}M_{\odot}over˙ start_ARG italic_M end_ARG > 10 start_POSTSUPERSCRIPT - 7 end_POSTSUPERSCRIPT italic_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT/yr (nearly all simulated stars) and stars with M˙>106M˙𝑀superscript106subscript𝑀direct-product\dot{M}>10^{-6}M_{\odot}over˙ start_ARG italic_M end_ARG > 10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT italic_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT/yr (8 stars). Wolf-Rayet stars are not included. The red dashed line shows the theoretical size of the WTS blown by a 5 Myr old star with M˙=106M˙𝑀superscript106subscript𝑀direct-product\dot{M}=10^{-6}M_{\odot}over˙ start_ARG italic_M end_ARG = 10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT italic_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT/yr and \varv=2500subscript\varv2500\varv_{\infty}=2500start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT = 2500 km/s. The blue dashed line shows the theoretical size of the WTS blown by a 5 Myr old star with M˙=107M˙𝑀superscript107subscript𝑀direct-product\dot{M}=10^{-7}M_{\odot}over˙ start_ARG italic_M end_ARG = 10 start_POSTSUPERSCRIPT - 7 end_POSTSUPERSCRIPT italic_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT/yr and \varv=2500subscript\varv2500\varv_{\infty}=2500start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT = 2500 km/s.

In Cygnus OB2, the distribution of stars is not homogeneous but peaked around the centre, as can be seen in the middle panel of Fig. 2. Therefore, the average distance between O stars increases as one moves away from the cluster centre. Fig. 9 shows the evolution of the average distance between the stars, d=(2/3)1/3N(r)1/3r𝑑superscript2313𝑁superscript𝑟13𝑟d=(2/3)^{1/3}N(r)^{-1/3}ritalic_d = ( 2 / 3 ) start_POSTSUPERSCRIPT 1 / 3 end_POSTSUPERSCRIPT italic_N ( italic_r ) start_POSTSUPERSCRIPT - 1 / 3 end_POSTSUPERSCRIPT italic_r, where N(r)𝑁𝑟N(r)italic_N ( italic_r ) is the number of O stars within the radius r𝑟ritalic_r in our simulated sample. This should be representative of Cygnus OB2. Efficient wind-wind interactions are expected only if a single star WTS is a sizeable fraction of the average distance between the stars. One concludes that, in our sample, only the most powerful O stars (M˙>106M/yr˙𝑀superscript106subscript𝑀direct-productyr\dot{M}>10^{-6}M_{\odot}/\mathrm{yr}over˙ start_ARG italic_M end_ARG > 10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT italic_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT / roman_yr) within a distance of a couple of parsecs from the cluster centre, can interact efficiently. The other, less powerful, stars can only interact within the inner few parsecs. This explains why wind-wind interactions are only efficient in the inner few parsecs in the simulation and dominated by the most powerful stars.

The argument exposed above is quite general and the conclusion depends only weakly on the parameters of the problem. In particular, even if the real configuration of the stellar population in Cygnus OB2 is slightly different than the one we simulate, or if it happens that the population was different in the past, or if the initial density is lower than the value assumed in the simulation, one would still not expect to see a large-scale WTS in the region.

Refer to caption
Figure 10: Result of a complementary simulation with 100 massive stars, from 24M24subscript𝑀direct-product24M_{\odot}24 italic_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT to 90M90subscript𝑀direct-product90M_{\odot}90 italic_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT, sampled from a Salpeter IMF and uniformly gathered within 3 pc. Top panel: 3D distribution of the stars. The size of the spheres scales as Lw1/5superscriptsubscript𝐿𝑤15L_{w}^{1/5}italic_L start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 1 / 5 end_POSTSUPERSCRIPT, with Lwsubscript𝐿𝑤L_{w}italic_L start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT the mechanical power of the star. Middle panel: Mach=1 isosurfaces at 2 Myr. Bottom panel: Mach number in the (x,y)𝑥𝑦(x,y)( italic_x , italic_y ) plane (z=0𝑧0z=0italic_z = 0) at 2 Myr.

As rule of thumb, a massive star cluster could create a large-scale WTS only if the cluster core is smaller than a few parsecs. However, this might not be a sufficient condition to generate a spherical wind termination shock. In order to investigate this last point, we ran a fiducial simulation of a more compact cluster, with a synthetic population of 100 massive stars generated from a Salpeter IMF and estimating the wind parameters using Eqs. 4,5. The spatial distribution is randomly generated assuming a uniform distribution of stars within 3 pc, in an ISM of initial density 10 cm-3. The total power of this fiducial cluster is 1.1×10381.1superscript10381.1\times 10^{38}1.1 × 10 start_POSTSUPERSCRIPT 38 end_POSTSUPERSCRIPT erg/s. If such a power would be injected as a point-like source, one would expect a spherical cluster WTS to develop up to 13 pc after 2 Myr. Fig. 10 shows the simulation result after 2 Myr. Clearly, wind-wind interactions are much more efficient in such a compact case than it is in the case of Cygnus OB2. These interactions create a number of trans-sonic sheets and jets, which are still able to escape from the core region. Therefore, even in such a compact case, the wind kinetic energy escapes the region of energy deposition in a very inhomogeneous manner and the shock surface shape is irregular. This highlights the difficulty to generate a spherical cluster wind termination shock from a population of O stars that is only statistically spherically symmetric.

A detailed study of the structure of the shocks and jets in such a compact cluster is left for future work.

Refer to caption
Figure 11: Wind-wind interactions in the inner region at 2 Myr. Top panel: magnitude of the velocity field. Middle panel: magnitude of the vorticity. Bottom panel: Mach number. Dark red shows regions of Mach number larger than 2. The largest WTS on the left and right are respectively blown by WR146 and WR144. The slices are taken at y=0𝑦0y=0italic_y = 0.

3.3 Consequences for high-energy astrophysics

The aim of the present study is to diagnose the shocks in the Cygnus OB2 star cluster, which are key to understanding the high-energy emission observed from the region. Our simulation shows that a scenario based on a large-scale spherical cluster WTS is excluded. On the other hand, the powerful Wolf-Rayet stars do blow strong large-scale spherical shocks which could in principle be sites of particle acceleration (e.g., Casse & Paul, 1980; Zirakashvili & Ptuskin, 2018). In the simulation, we ran the Wolf-Rayet phase during 400 kyr, which corresponds to the typical lifetime of central Helium burning in a massive star (Woosley, 2019; Higgins et al., 2021, e.g.). Given that some stars might only reach the WR stage later in their central He burning (e.g., Josiek et al., 2024), this can be considered a upper limit for the typical WR lifetime. At this point, the most powerful WR star has expanded a shock up to a radius of 5 pc.

We can work out an upper bound on the maximum achievable energy of a proton accelerated by the WTS around WR144. Efficient particle acceleration requires a super-Alfvénic shock, which limits the upstream magnetic field as BRwts<\varvM˙𝐵subscript𝑅wtssubscript\varv˙𝑀BR_{\text{wts}}<\sqrt{\varv_{\infty}\dot{M}}italic_B italic_R start_POSTSUBSCRIPT wts end_POSTSUBSCRIPT < square-root start_ARG start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT over˙ start_ARG italic_M end_ARG end_ARG. Besides, the Hillas criterion (Hillas, 1984) sets an upper bound on the maximum energy as Emax<qBRwts\varv/csubscript𝐸max𝑞𝐵subscript𝑅wtssubscript\varv𝑐E_{\text{max}}<qBR_{\text{wts}}\varv_{\infty}/citalic_E start_POSTSUBSCRIPT max end_POSTSUBSCRIPT < italic_q italic_B italic_R start_POSTSUBSCRIPT wts end_POSTSUBSCRIPT start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT / italic_c, which is actually the limit imposed by adiabatic losses in the radial wind and therefore cannot be overcome by enhanced magnetic field amplification downstream of the shock e.g. due to the Cranfill effect or nontrivial interactions with the nearby O star winds. We eventually get Emax<1.5Zsubscript𝐸max1.5𝑍E_{\text{max}}<1.5Zitalic_E start_POSTSUBSCRIPT max end_POSTSUBSCRIPT < 1.5 italic_Z PeV for the wind parameters of WR144 (M˙=2.4×105M˙𝑀2.4superscript105subscript𝑀direct-product\dot{M}=2.4\times 10^{-5}M_{\odot}over˙ start_ARG italic_M end_ARG = 2.4 × 10 start_POSTSUPERSCRIPT - 5 end_POSTSUPERSCRIPT italic_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT/yr, \varv=3500subscript\varv3500\varv_{\infty}=3500start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT = 3500 km/s). This should be taken as a strict upper bound and not an estimate, for we know that the limiting factors on the maximum achievable non-thermal energy around a spherical shock are more stringent than the Hillas limit (Morlino et al., 2021), and the efficiency of particle acceleration at quasi-perpendicular shocks is debated (Caprioli & Spitkovsky, 2014; Xu et al., 2020; Kumar & Reville, 2021). On the other hand, a regime of perpendicular diffusion around the WTS might enhance the maximum achievable energy (Jokipii & Morfill, 1987), although this requires a detailed treatment of the transport along the shock surface in a Parker-spiral magnetic field (Kamijima & Ohira, 2022). In the end, even if the WR stars in Cygnus OB2 are unlikely to account alone for the ultra-high energy gamma-rays detected from the Cygnus region, they might provide a non-negligible flux of high-energy particles at the level of a fraction of the wind power up to hundreds of TeV.

In addition to the individual stellar wind termination shocks, there is still some level of collective interaction in the inner few parsecs, mostly due to the stars WR144, Cyg OB2 #7, Cyg OB2 #8B+#8C, Cyg OB2 #22 and a couple of less powerful O stars. As can be seen in Fig. 11, these stars are close enough to interact efficiently, creating structured trans-sonic sheets and jets. The shocked plasma is characterised by a high velocity (hundreds of km/s) and high vorticity, which suggests a high level of hydrodynamic turbulence with the formation of large-scale eddies. Particle acceleration in a collection of strong shocks and strong turbulence in Cygnus OB2 was investigated by Bykov & Kalyashova (2022). The authors claimed that this model was able to reproduce the observations up to the highest energies, assuming a cluster size of 55 pc, a turbulent velocity of 1500 km/s in the entire region. Indeed, these parameters together with the assumption of a 15 µG magnetic field provide a maximum energy of about 4 PeV, assuming that the Hillas limit applies, which is not unlikely in this acceleration scenario (Vieu et al., 2022b). However, Fig. 11 shows that this choice of parameters is not realistic. The turbulent region extends at most over 10 pc and the mean velocity is only a few hundreds of km/s, which would provide a maximum energy of at most 100 TeV. Furthermore, the purely stochastic model of Bykov (2001) does not properly describe the plasma in star cluster environments, because the flow remains laminar inside the stellar wind cavities, which inhibits re-acceleration effects (Vieu et al., 2024). Finally, we stress again that direct wind-wind collisions are not expected to enhance the maximum energy achieved by accelerated particles when a realistic flow configuration in multi-dimensions is taken into account (Vieu et al., 2020, 2024).

4 Conclusions

We analysed the stellar population of Cygnus OB2 as provided by observations, highlighting that i) most of the cluster power is actually provided by 3 WR stars and 5 O stars, and ii) the distribution of the massive stars extends up to more than 10 pc without a strong mass segregation. Cygnus OB2 must therefore be considered as a “loose association” rather than a “compact cluster”. In this case, one-dimensional modelling of the cluster feedback is not expected to provide reliable results. In particular, assuming a continuous deposition of thermal energy within a spherical region is an oversimplification. Detailed hydrodynamic simulations are then necessary to properly understand the stellar feedback on interstellar scales.

To this aim, we simulated the 40 most powerful stars in Cygnus OB2. We were able to resolve each massive star individually in the core, while allowing their feedback to develop features up to several tens of parsecs. Since the uncertainties on the measurements of parallaxes and stellar wind parameters prevent us from modelling the exact stellar population, we set up a cluster whose statistical properties match that of Cygnus OB2. Our 3D reconstruction aims to reproduce the radial distribution of massive stars, although it is expected to provide a slightly more compact configuration.

Despite these shortcomings, the simulation results remain reliable on the qualitative level. In particular, the absence of a large-scale cluster WTS in the simulation implies the absence of such a cluster WTS in Cygnus OB2. This result is a natural consequence of the loose distribution of the most powerful stars in space. Since the average interstellar separation is typically much larger than the size of an individual stellar WTS, there is a relatively low level of wind-wind interactions. The kinetic energy of the winds is found to thermalise in a highly inhomogeneous manner which cannot be accounted for using standard spherically symmetric models.

In order to further probe the formation of a cluster WTS, we ran a complementary simulation of a synthetic population of 100 massive stars, clustered in a core with a radius 4 times smaller than that of Cygnus OB2. We found that even in this case, non-trivial wind-wind interactions on sub-parsec scales prevent the formation of a spherical WTS. Wind-wind collisions rather produce trans-sonic sheets and jets which disperse the kinetic energy to large distances, hindering the expansion of a collective radial outflow beyond the cluster core.

The physical reason underlying the absence of a cluster WTS in Cygnus OB2 is thus well understood and easily generalisable. In principle it would equally apply in a more refined simulation, e.g. including magnetic fields, thermal conduction and cooling, or for a different cluster, provided the distribution of stars in space is similar. This is the only relevant parameter of the problem. The stars in Cygnus OB2 are – and most likely always were – too dispersed to enable strong collective effects.

Some level of collective interactions is nevertheless found in the inner few parsecs of our simulation, in particular after the onset of WR stars. Similarly to the case of the more compact cluster, collisions between a few individual stellar WTS produce focused trans-sonic outflows. This non-trivial feedback creates an intricate environment with multiple shocks. Although these shocks do provide channels for particle acceleration, they are not expected to contribute beyond energies of a few hundred TeV and are therefore unlikely to be the origin of the PeV gamma-rays detected from the Cygnus region by the LHAASO observatory (LHAASO Collaboration, 2024).

Acknowledgements

We thank the referee and the editor for useful suggestions. This work made use of the MHD code PLUTO. The simulation was performed on the HPC system Raven at the Max Planck Computing and Data Facility. We acknowledge J.S. Wang for helpful discussions on the implementation of the simulation. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. 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. CJKL, AS, and VR gratefully acknowledge support by the German Deutsche Forschungsgemeinschaft, DFG in the form of an Emmy Noether Research Group – Project-ID 445674056 (SA4064/1-1, PI Sander) and the Federal Ministry of Education and Research (BMBF) and the Baden-Württemberg Ministry of Science as part of the Excellence Strategy of the German Federal and State Governments.

Data Availability

The initial conditions and output of the simulation may be shared on reasonable request to the corresponding author.

Appendix A List of simulated stars

The stars set-up in our simulation are listed in Table 12. We stress that this only partly reproduces the observations, as measurement uncertainties on parallaxes and wind parameters called for some extrapolations which are described in Section 2.2 and Section 2.3. Fig. 12 provides a 3D view of the simulated cluster.

x𝑥xitalic_x [pc] y𝑦yitalic_y [pc] z𝑧zitalic_z [pc] M˙˙𝑀\dot{M}over˙ start_ARG italic_M end_ARG [Msubscript𝑀direct-productM_{\odot}italic_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT/yr] \varvsubscript\varv\varv_{\infty}start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT [km/s] Method Stellar type Comment
6.67 11.89 0.00 3.24E-05 1440 Inferred WR (after 1.6 Myr) Y coordinate shifted from 14.5 pc to 11.9 pc
-10.63 -1.93 0.00 1.26E-05 2900 Inferred WR (after 1.6 Myr) X coordinate shifted from 13.6 pc to 10.6 pc
6.48 1.61 1.29 2.40E-05 3500 Inferred WR (after 1.6 Myr)
2.16 2.03 -2.16 3.02E-06 400 Inferred BHG Cygnus OB2 #12
0.19 0.77 -0.45 7.41E-06 3188 Inferred O Cygnus OB2 #7
1.64 2.32 11.66 6.34E-06 2403 Inferred O Cygnus OB2 #11
0.48 1.13 2.87 3.30E-06 3229 Inferred O Cygnus OB2 #22
-0.06 -0.10 0.13 4.82E-06 2541 Inferred O Merging of OB2 #8B and OB2 #8C
-1.03 -1.39 5.25 1.28E-06 2352 Inferred O
-2.45 7.76 -0.03 5.92E-07 2372 Inferred O
-1.19 -6.96 4.54 3.93E-07 2637 Inferred O
5.70 -5.41 3.03 5.21E-07 2280 Inferred O
0.06 0.48 0.74 4.16E-07 2331 Inferred O
9.02 3.19 0.81 2.97E-07 2362 Inferred O
2.84 -0.42 -0.45 2.24E-07 2474 Inferred O
-2.80 0.97 6.48 2.18E-07 2474 Inferred O
-2.58 -1.42 -2.42 1.79E-07 2352 Inferred O
-7.60 2.48 10.89 1.51E-07 2484 Inferred O
-0.13 -0.06 2.74 2.35E-07 1862 Inferred O
10.21 -1.87 -2.87 1.62E-07 2158 Inferred O
0.10 0.55 -3.06 1.27E-07 2352 Inferred O
0.39 1.22 11.57 1.92E-07 1882 Inferred O
-0.77 -3.35 -10.92 1.03E-07 2474 Inferred O
-4.96 -1.77 0.32 1.08E-07 2270 Inferred O
4.48 -0.35 -0.35 9.54E-08 2372 Inferred O
-0.03 -0.06 5.61 9.41E-08 2352 Inferred O
1.39 -0.68 -0.35 8.14E-08 2505 Inferred O
0.48 1.51 2.19 1.31E-07 1791 Inferred O
3.00 2.84 -3.54 1.42E-06 2599 IMF, M=48.5𝑀48.5M=48.5italic_M = 48.5 O
10.31 2.58 -2.32 1.04E-06 2580 IMF, M=44.2𝑀44.2M=44.2italic_M = 44.2 O
-2.32 8.31 3.06 7.70E-07 2563 IMF, M=40.7𝑀40.7M=40.7italic_M = 40.7 O
6.73 -1.19 -0.58 5.78E-07 2548 IMF, M=37.7𝑀37.7M=37.7italic_M = 37.7 O
-2.58 6.41 5.38 4.38E-07 2534 IMF, M=35.2𝑀35.2M=35.2italic_M = 35.2 O
-4.41 -2.45 1.74 3.35E-07 2521 IMF, M=33.1𝑀33.1M=33.1italic_M = 33.1 O
-2.64 -3.25 3.51 2.02E-07 2498 IMF, M=29.6𝑀29.6M=29.6italic_M = 29.6 O
1.48 1.55 2.29 1.59E-07 2488 IMF, M=28.1𝑀28.1M=28.1italic_M = 28.1 O
-2.93 1.39 -0.71 1.26E-07 2479 IMF, M=26.8𝑀26.8M=26.8italic_M = 26.8 O
-0.77 3.32 -7.60 1.01E-07 2470 IMF, M=25.6𝑀25.6M=25.6italic_M = 25.6 O
-0.39 -2.13 3.99 8.14E-08 2462 IMF, M=24.6𝑀24.6M=24.6italic_M = 24.6 O
0.19 0.00 -0.61 1.00E-06 2500 Fiducial O Putative SN progenitor
Table 2: List of the simulated stars. In the Method column, “Inferred” means that the wind parameters have been obtained from measured values of the wind-strength parameter and effective temperature. “IMF” means that the wind parameters have been generated from a synthetic population following an initial mass function of index 1.39, and we provide the initial mass in the Table. See Sections 2.22.3 for details.
Refer to caption
Figure 12: Simulated cluster. Red: Wolf-Rayet stars. Yellow: O stars above 1037superscript103710^{37}10 start_POSTSUPERSCRIPT 37 end_POSTSUPERSCRIPT erg/s. Black: Remaining O stars from the 52 with inferred parameters. Blue: O stars from the synthetic population. The size of the spheres scales as Lw1/5superscriptsubscript𝐿𝑤15L_{w}^{1/5}italic_L start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 1 / 5 end_POSTSUPERSCRIPT, with Lwsubscript𝐿𝑤L_{w}italic_L start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT the mechanical power of the star.

References

  • Abeysekara et al. (2021) Abeysekara A. U., et al., 2021, Nature Astronomy, 5, 465
  • Ackermann et al. (2011) Ackermann M., et al., 2011, Science, 334, 1103
  • Agertz et al. (2013) Agertz O., Kravtsov A. V., Leitner S. N., Gnedin N. Y., 2013, ApJ, 770, 25
  • Albacete-Colombo et al. (2023) Albacete-Colombo J. F., et al., 2023, ApJS, 269, 14
  • Astiasarain et al. (2023) Astiasarain X., Tibaldo L., Martin P., Knödlseder J., Remy Q., 2023, A&A, 671, A47
  • Badmaev et al. (2022) Badmaev D. V., Bykov A. M., Kalyashova M. E., 2022, MNRAS, 517, 2818
  • Berlanas et al. (2018) Berlanas S. R., Herrero A., Comerón F., Pasquali A., Bertelli Motta C., Sota A., 2018, A&A, 612, A50
  • Berlanas et al. (2019) Berlanas S. R., Wright N. J., Herrero A., Drew J. E., Lennon D. J., 2019, MNRAS, 484, 1838
  • Berlanas et al. (2020) Berlanas S. R., et al., 2020, A&A, 642, A168
  • Bica et al. (2003) Bica E., Bonatto C., Dutra C. M., 2003, A&A, 405, 991
  • Bykov (2001) Bykov A. M., 2001, Space Sci. Rev., 99, 317
  • Bykov & Fleishman (1992) Bykov A. M., Fleishman G. D., 1992, MNRAS, 255, 269
  • Bykov & Kalyashova (2022) Bykov A. M., Kalyashova M. E., 2022, Advances in Space Research, 70, 2685
  • Bykov et al. (2013) Bykov A. M., Gladilin P. E., Osipov S. M., 2013, MNRAS, 429, 2755
  • Cantó et al. (2000) Cantó J., Raga A. C., Rodríguez L. F., 2000, ApJ, 536, 896
  • Cao et al. (2021) Cao Z., et al., 2021, Nature, 594, 33
  • Caprioli & Spitkovsky (2014) Caprioli D., Spitkovsky A., 2014, ApJ, 783, 91
  • Casse & Paul (1980) Casse M., Paul J. A., 1980, ApJ, 237, 236
  • Castor et al. (1975) Castor J., McCray R., Weaver R., 1975, ApJ, 200, L107
  • Chevalier & Clegg (1985) Chevalier R. A., Clegg A. W., 1985, Nature, 317, 44
  • Clark et al. (2005) Clark J. S., Negueruela I., Crowther P. A., Goodwin S. P., 2005, A&A, 434, 949
  • Clark et al. (2012) Clark J. S., Najarro F., Negueruela I., Ritchie B. W., Urbaneja M. A., Howarth I. D., 2012, A&A, 541, A145
  • Comerón & Pasquali (2012) Comerón F., Pasquali A., 2012, A&A, 543, A101
  • Comerón et al. (2002) Comerón F., et al., 2002, A&A, 389, 874
  • Drew et al. (2008) Drew J. E., Greimel R., Irwin M. J., Sale S. E., 2008, MNRAS, 386, 1761
  • Eenens & Williams (1994) Eenens P. R. J., Williams P. M., 1994, MNRAS, 269, 1082
  • El-Badry et al. (2019) El-Badry K., Ostriker E. C., Kim C.-G., Quataert E., Weisz D. R., 2019, MNRAS, 490, 1961
  • Ferrand & Marcowith (2010) Ferrand G., Marcowith A., 2010, A&A, 510, A101
  • Gabici et al. (2019) Gabici S., Evoli C., Gaggero D., Lipari P., Mertsch P., Orlando E., Strong A., Vittino A., 2019, International Journal of Modern Physics D, 28, 1930022
  • Gaia Collaboration et al. (2023) Gaia Collaboration et al., 2023, A&A, 674, A1
  • Gatto et al. (2017) Gatto A., et al., 2017, MNRAS, 466, 1903
  • Geen et al. (2021) Geen S., Bieri R., Rosdahl J., de Koter A., 2021, MNRAS, 501, 1352
  • Gupta et al. (2018a) Gupta S., Nath B. B., Sharma P., Eichler D., 2018a, MNRAS, 473, 1537
  • Gupta et al. (2018b) Gupta S., Nath B. B., Sharma P., 2018b, MNRAS, 479, 5220
  • Gupta et al. (2020) Gupta S., Nath B. B., Sharma P., Eichler D., 2020, MNRAS, 493, 3159
  • Hamann et al. (2008) Hamann W.-R., Oskinova L. M., Feldmeier A., 2008, in Hamann W.-R., Feldmeier A., Oskinova L. M., eds, Clumping in Hot-Star Winds. p. 75
  • Hanson (2003) Hanson M. M., 2003, ApJ, 597, 957
  • Härer et al. (2023) Härer L. K., Reville B., Hinton J., Mohrmann L., Vieu T., 2023, A&A, 671, A4
  • Hawcroft et al. (2023) Hawcroft C., et al., 2023, arXiv e-prints, p. arXiv:2303.12165
  • Higgins et al. (2021) Higgins E. R., Sander A. A. C., Vink J. S., Hirschi R., 2021, MNRAS, 505, 4874
  • Hillas (1984) Hillas A. M., 1984, ARA&A, 22, 425
  • Holgado et al. (2018) Holgado G., et al., 2018, A&A, 613, A65
  • Jokipii & Morfill (1987) Jokipii J. R., Morfill G., 1987, ApJ, 312, 170
  • Josiek et al. (2024) Josiek J., Ekström S., Sander A. A. C., 2024, arXiv e-prints, p. arXiv:2404.14488
  • Kamijima & Ohira (2022) Kamijima S. F., Ohira Y., 2022, Phys. Rev. D, 106, 123025
  • Kim et al. (2017) Kim C.-G., Ostriker E. C., Raileanu R., 2017, ApJ, 834, 25
  • Klepach et al. (2000) Klepach E. G., Ptuskin V. S., Zirakashvili V. N., 2000, Astroparticle Physics, 13, 161
  • Knödlseder (2000) Knödlseder J., 2000, A&A, 360, 539
  • Krause et al. (2013) Krause M., Fierlinger K., Diehl R., Burkert A., Voss R., Ziegler U., 2013, A&A, 550, A49
  • Kudritzki & Puls (2000) Kudritzki R.-P., Puls J., 2000, ARA&A, 38, 613
  • Kumar & Reville (2021) Kumar N., Reville B., 2021, ApJ, 921, L14
  • LHAASO Collaboration (2024) LHAASO Collaboration 2024, Science Bulletin, 69, 449
  • Lancaster et al. (2021) Lancaster L., Ostriker E. C., Kim J.-G., Kim C.-G., 2021, ApJ, 914, 90
  • Maíz Apellániz et al. (2020) Maíz Apellániz J., Crespo Bellido P., Barbá R. H., Fernández Aranda R., Sota A., 2020, A&A, 643, A138
  • Massey & Thompson (1991) Massey P., Thompson A. B., 1991, AJ, 101, 1408
  • Menchiari et al. (2024) Menchiari S., Morlino G., Amato E., Bucciantini N., Beltrán M. T., 2024, arXiv e-prints, p. arXiv:2402.07784
  • Mignone et al. (2007) Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, ApJS, 170, 228
  • Morlino et al. (2021) Morlino G., Blasi P., Peretti E., Cristofari P., 2021, MNRAS, 504, 6096
  • Muntean et al. (2009) Muntean V., Moffat A. F. J., Chené A. N., de La Chevrotière A., 2009, MNRAS, 399, 1977
  • Orellana et al. (2021) Orellana R. B., De Biasi M. S., Paíz L. G., 2021, MNRAS, 502, 6080
  • Oskinova et al. (2016) Oskinova L. M., Kubátová B., Hamann W.-R., 2016, J. Quant. Spectrosc. Radiative Transfer, 183, 100
  • Oskinova et al. (2017) Oskinova L. M., Huenemoerder D. P., Hamann W. R., Shenar T., Sander A. A. C., Ignace R., Todt H., Hainich R., 2017, ApJ, 845, 39
  • Pittard et al. (2021a) Pittard J. M., Romero G. E., Vila G. S., 2021a, MNRAS, 504, 4204
  • Pittard et al. (2021b) Pittard J. M., Wareing C. J., Kupilas M. M., 2021b, MNRAS, 508, 1768
  • Povich & Whitney (2010) Povich M. S., Whitney B. A., 2010, ApJ, 714, L285
  • Puls et al. (2005) Puls J., Urbaneja M. A., Venero R., Repolust T., Springmann U., Jokuthy A., Mokiem M. R., 2005, A&A, 435, 669
  • Puls et al. (2008) Puls J., Vink J. S., Najarro F., 2008, A&ARv, 16, 209
  • Ramachandran et al. (2018) Ramachandran V., Hamann W. R., Hainich R., Oskinova L. M., Shenar T., Sander A. A. C., Todt H., Gallagher J. S., 2018, A&A, 615, A40
  • Rivero González et al. (2012) Rivero González J. G., Puls J., Massey P., Najarro F., 2012, A&A, 543, A95
  • Rogers & Pittard (2013) Rogers H., Pittard J. M., 2013, MNRAS, 431, 1337
  • Rygl et al. (2012) Rygl K. L. J., et al., 2012, A&A, 539, A79
  • Sander et al. (2012) Sander A., Hamann W. R., Todt H., 2012, A&A, 540, A144
  • Sander et al. (2019) Sander A. A. C., Hamann W. R., Todt H., Hainich R., Shenar T., Ramachandran V., Oskinova L. M., 2019, A&A, 621, A92
  • Schneider et al. (2018) Schneider F. R. N., et al., 2018, A&A, 618, A73
  • Seo et al. (2018) Seo J., Kang H., Ryu D., 2018, Journal of Korean Astronomical Society, 51, 37
  • Simón-Díaz et al. (2011) Simón-Díaz S., Castro N., Garcia M., Herrero A., Markova N., 2011, Bulletin de la Societe Royale des Sciences de Liege, 80, 514
  • Sota et al. (2011) Sota A., Maíz Apellániz J., Walborn N. R., Alfaro E. J., Barbá R. H., Morrell N. I., Gamen R. C., Arias J. I., 2011, ApJS, 193, 24
  • Uyanıker, B. et al. (2001) Uyanıker, B. Fürst, E. Reich, W. Aschenbach, B. Wielebinski, R. 2001, A&A, 371, 675
  • Vieu & Reville (2023) Vieu T., Reville B., 2023, MNRAS, 519, 136
  • Vieu et al. (2020) Vieu T., Gabici S., Tatischeff V., 2020, MNRAS, 494, 3166
  • Vieu et al. (2022a) Vieu T., Gabici S., Tatischeff V., Ravikularaman S., 2022a, MNRAS, 512, 1275
  • Vieu et al. (2022b) Vieu T., Reville B., Aharonian F., 2022b, MNRAS, 515, 2256
  • Vieu et al. (2024) Vieu T., Härer L., Reville B., 2024, MNRAS, p. stae1039
  • Weaver et al. (1977) Weaver R., McCray R., Castor J., Shapiro P., Moore R., 1977, ApJ, 218, 377
  • Woosley (2019) Woosley S. E., 2019, ApJ, 878, 49
  • Wright et al. (2015) Wright N. J., Drew J. E., Mohr-Smith M., 2015, MNRAS, 449, 741
  • Wright et al. (2016) Wright N. J., Bouy H., Drew J. E., Sarro L. M., Bertin E., Cuillandre J.-C., Barrado D., 2016, MNRAS, 460, 2593
  • Xu et al. (2020) Xu R., Spitkovsky A., Caprioli D., 2020, ApJ, 897, L41
  • Zhang et al. (2024) Zhang S., et al., 2024, AJ, 167, 220
  • Zhekov (2017) Zhekov S. A., 2017, MNRAS, 472, 4374
  • Zirakashvili & Ptuskin (2018) Zirakashvili V. N., Ptuskin V. S., 2018, Astroparticle Physics, 98, 21