Hydrodynamic simulation of Cygnus OB2: the absence of a cluster wind termination shock
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 rays1 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 ( 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 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 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 () than what is required to explain the gamma-ray observations ( 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 ( 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 ( pc) inside a numerical box that is large enough to contain the large-scale superbubble ( 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.
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 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 (mag, 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 kpc and a main complex at kpc. 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 deg away from each other in the projected plane on the sky.
In the present work, we consider the area within 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 kpc, motivated by the recent analysis from Maíz Apellániz et al. (2020), who obtained a Gaia-based distance of kpc 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 kpc.
The age of Cygnus OB2 is ill-defined as the region is not coeval. Estimates range between 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 and (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 and 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 (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 (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 and 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 , the stellar radius and the wind-strength parameter
(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, , the combined wind strength parameter can be robustly constrained from the spectral fit but not the individual values of and the mass-loss rate . As we wish to know the mechanical luminosity
(2) |
for our simulation, we need an additional estimate for to calculate 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 and . The recent determination by Hawcroft et al. (2023) for O stars at Galactic metallicity yields the formula
(3) |
which we use to estimate values of and thus obtain both and , avoiding any explicit assumptions about the spectroscopic mass. We note that the values for derived by Berlanas et al. (2020) and also the inferred 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 (Wright et al., 2015). We pick stars with masses between 20 and 50 in order to avoid introducing unusually powerful stars which are not seen in the observations. On the other hand, starting at 20 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 , we use the prescriptions of Seo et al. (2018) (in the case of non-rotating single stars):
(4) | |||
(5) |
In the end we obtain a population of O stars with a mechanical power of 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 km/s. When including the WR stars (cf. Fig. 1), we obtain a total mechanical power of 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 and .
WR144 | WR145 | WR146 | Cygnus OB2 #7 | Cygnus OB2 | Cygnus OB2 #11 | Cygnus OB2 #22 | O stars | O stars | |
#8B+#8C | (total) | (simulated) | |||||||
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
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 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 . 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 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).
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 , such that stars located outside of the cube of half edge-length will be excluded from the simulation. Our computational resources allow us to set 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 erg/s and 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 coordinate of WR145 from 14.5 pc to 11.9 pc, and the coordinate of WR146 from 13.6 pc to 10.6 pc. For these two stars, we set 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 . This represents a population of 20 O stars which has a negligible mechanical power ( 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 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 , corresponding to a proper motion of order . 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 /yr, 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.
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.
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 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 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 cm-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 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 and mass-loss rate :
(6) |
where the parameters and 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: , with km/s for O stars and km/s for WR stars. During the simulation, the stellar winds are set up within spheres of radius 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 , with the initial pressure. This condition is necessary to accurately expand wind-blown bubbles (Pittard et al., 2021b). For the stars with the lowest product, we get = 0.3, while for the most powerful stars, which dominate the stellar feedback at large scales, we have .
The initial () 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: , . 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 /yr and 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 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.
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 ( 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 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 , where 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.
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 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 ) 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):
(7) |
where is the total power of the cluster, is the initial density and 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, with the pressure inside the superbubble. This provides:
(8) |
where is the mass-loss rate of the isolated star and its wind velocity. We have normalised the parameters to values typical of Cygnus OB2. In fact, 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.
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, , where is the number of O stars within the radius 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 () 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.
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 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.
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 . Besides, the Hillas criterion (Hillas, 1984) sets an upper bound on the maximum energy as , 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 PeV for the wind parameters of WR144 (/yr, 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.
[pc] | [pc] | [pc] | [/yr] | [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, | O | |
10.31 | 2.58 | -2.32 | 1.04E-06 | 2580 | IMF, | O | |
-2.32 | 8.31 | 3.06 | 7.70E-07 | 2563 | IMF, | O | |
6.73 | -1.19 | -0.58 | 5.78E-07 | 2548 | IMF, | O | |
-2.58 | 6.41 | 5.38 | 4.38E-07 | 2534 | IMF, | O | |
-4.41 | -2.45 | 1.74 | 3.35E-07 | 2521 | IMF, | O | |
-2.64 | -3.25 | 3.51 | 2.02E-07 | 2498 | IMF, | O | |
1.48 | 1.55 | 2.29 | 1.59E-07 | 2488 | IMF, | O | |
-2.93 | 1.39 | -0.71 | 1.26E-07 | 2479 | IMF, | O | |
-0.77 | 3.32 | -7.60 | 1.01E-07 | 2470 | IMF, | O | |
-0.39 | -2.13 | 3.99 | 8.14E-08 | 2462 | IMF, | O | |
0.19 | 0.00 | -0.61 | 1.00E-06 | 2500 | Fiducial | O | Putative SN progenitor |
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