# Cosmic ray production in superbubbles

T. Vieu,1 S. Gabici,1 V. Tatischeff,2 S. Ravikularaman1
1Université de Paris, CNRS, Astroparticule et Cosmologie, F-75013 Paris,
France
2Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We compute the production of cosmic rays in the dynamical superbubble produced
by a cluster of massive stars. Stellar winds, supernova remnants and
turbulence are found to accelerate particles so efficiently that the nonlinear
feedback of the particles must be taken into account in order to ensure that
the energy balance is not violated. High energy particles do not scatter
efficiently on the turbulence and escape quickly after each supernova
explosion, which makes both their intensity inside the bubble and injection in
the interstellar medium intermittent. On the other hand, the stochastic
acceleration of low energy particles hardens the spectra at GeV energies.
Because cosmic rays damp the turbulence cascade, this hardening is less
pronounced when nonlinearities are taken into account. Nevertheless, spectra
with hard components extending up to 1 to 10 GeV and normalised to an energy
density of 1 to 100 eV cm-3 are found to be typical signatures of cosmic rays
produced in superbubbles. Efficient shock reacceleration within compact
clusters is further shown to produce hard, slightly concave spectra, while the
presence of a magnetised shell is shown to enhance the confinement of cosmic
rays in the bubble and therefore the collective plasma effects acting on them.
We eventually estimate the overall contribution of superbubbles to the
galactic cosmic ray content and show typical gamma-ray spectra expected from
hadronic interactions in superbuble shells. In both cases, a qualitative
agreement with observations is obtained.

###### keywords:

acceleration of particles – shock waves – cosmic rays

††pubyear: 2019††pagerange: Cosmic ray production in superbubbles–A

##  1 Introduction

Superbubbles (SBs) are expanding cavities carved inside dense molecular clouds
by the feedback of clustered massive stars. It is thought that most massive
stars are born in either compact or loose associations, in which they spend
most of their lives until they explode as supernovae (Higdon & Lingenfelter,
2005). On the other hand, massive stars have long been thought to be the major
sources of galactic cosmic rays (CR), in particular via diffusive shock
acceleration (DSA) at wind termination shocks (WTS) and supernova remnants
(SNR) (e.g. Cesarsky & Montmerle, 1983; Gabici et al., 2019).

Inside SBs, star winds and supernovae combine to inject mechanical energy in
the surrounding medium. It was early realised (e.g. Montmerle, 1979) that this
energy could be available for efficient particle acceleration, in particular
due to collective plasma effects (Bykov & Fleishman, 1992b). This scenario is
an alternative to the standard paradigm of CR acceleration in isolated SNR. A
number of phenomenological studies (e.g. Higdon et al., 1998; Parizot & Drury,
1999; Higdon & Lingenfelter, 2006; Tatischeff & Gabici, 2018) support the SB
origin of CRs. The composition of galactic CRS also points to an origin of
these particles in SBs (Tatischeff et al., 2021). Besides, particle
reacceleration at multiple shocks and strong turbulence inside large
superbubbles is expected to be more efficient than the acceleration around
single supernovae, such that one can hope to push the maximum energy above the
spectral “knee”, as required to explain the transition between the galactic
and extragalactic components of the CR spectrum measured in the vicinity of
the Earth. Although this scenario for the origin of CRs has been explored over
the years (see the reviews by Parizot et al., 2004; Bykov, 2014; Lingenfelter,
2018; Bykov et al., 2020, and references therein), a renewed interest in it
rised after the recent detections in gamma-rays of star clusters and SBs, both
in the Milky Way (MW) and in the Large Magellanic Cloud (LMC) (Aharonian et
al., 2007; Abramowski et al., 2012; Ackermann et al., 2011; Abramowski et al.,
2015; Katsuta et al., 2017; Aharonian et al., 2019). This confirms that
stellar clusters and SBs are important CR sources. However, few attempts have
been made to model them in a self-consistent way.

Klepach et al. (2000) considered CR acceleration either at the collective
termination shock produced by the star winds in a compact cluster, or around
the multiple individual winds in the case of a loose cluster. They found that
the maximum CR energy could be increased by up to two orders of magnitude in
these systems compared to the acceleration at a single shock, and that the
more efficient acceleration would lead to hard spectra. Efforts to model the
acceleration of particles in interacting shocks, such as colliding winds or
SNRs, have been made in a series of paper by Bykov and collaborators (Bykov et
al., 2013; Bykov et al., 2018). It was shown that very hard spectra could be
expected in the idealised case of stationary plane infinite shocks. An attempt
to include geometrical effects has been performed in Bykov et al. (2015), and
an insight in the time-dependent aspects was discussed in Vieu et al. (2020).
Models of CR acceleration at the collective wind of a compact cluster were
developed recently by Gupta et al. (2020) and Morlino et al. (2021). Both
studies showed that the maximum energy of accelerated protons could exceed 1
PeV if the magnetic field close to the termination shock is efficiently
amplified, which is required to better confine the particles around the shock.
The stochastic acceleration in the turbulent medium was not included in these
models, even though a substantial fraction of the stellar energetic input is
believed to be stored in turbulent MHD waves (e.g. Butt & Bykov, 2008;
Gallegos-Garcia et al., 2020).

In the 90’s, Bykov et al. used renormalisation methods (Bykov & Toptygin,
1993) to model the acceleration of particles inside young and turbulent SBs
obtaining interesting results regarding the shape of the spectrum and break
energies (Bykov et al., 1995; Bykov, 2001). This stochastic model includes the
effect of the strong, possibly supersonic, turbulence, and also accounts for
the particle feedback onto turbulence, making it the first - and up to now, to
our knowledge, only - self-consistent modelling of CR acceleration in SB
environments. A time-dependent model was later developed by Ferrand &
Marcowith (2010), who considered the reacceleration of the particles by
successive SNR shocks sweeping the entire SB volume, as well as the modulation
of the spectrum in weak turbulence. However, the model of Ferrand & Marcowith
(2010) neglected the backreaction of the accelerated particles on the plasma
and assumed that stellar winds and energy losses would provide negligible
contributions to the nonthermal particle population. More recently, Tolksdorf
et al. (2019) solved the CR transport around SB cavities. They showed that
interstellar particles could be efficiently reaccelerated in the turbulent
interior. The primary acceleration at shock waves was however disregarded.

In the present work, we seek to develop a self-consistent model of CR
acceleration in SBs, including all the relevant physical ingredients, namely
stellar winds, SNRs, turbulence, energy losses, and also accounting for the
backreaction of the particles on the shocks and the turbulence. Section 2
discusses the SB environment. The acceleration mechanisms are modelled in
Section 3. Section 4 details the main results regarding the CR energetics,
intermittency and spectra inside a typical SB. The model is then refined in
Section 5 by considering two different regions, which can either account for a
compact cluster inside the SB or for the presence of a magnetised shell around
the SB. We estimate the SB contribution to the galactic CRs in Section 6 and
compute typical gamma-ray signatures in Section 7. We eventually conclude in
Section 8.

In order to improve the clarity of the paper, a summary of the notations used
throughout this work is provided in Table 2.

##  2 Superbubble properties

###  2.1 Bubble structure and stellar power

We consider a cluster of hundreds of massive ($8-150~{}M_{\odot}$) stars
distributed according to an initial mass function of index 2.3 (Salpeter,
1955). The mass of a star can be connected to its lifetime by using the
numerical fit from the stellar evolution simulations performed by Limongi &
Chieffi (2006):
$\log_{10}(t~{}\textrm{[yr]})=9.598-2.879\log_{10}M+0.6679(\log_{10}M)^{2}$,
where $M$ is the initial mass of the star in solar masses. Assuming that all
stars in a cluster are born at the same time $t=0$, it is possible to
implement a Monte Carlo sampling to simulate the times at which they explode
into supernovae. One of such realisations for a cluster of 100 massive stars
is shown in Figure 1. Such a number of massive stars represents the typical
content of superbubbles (Lingenfelter, 2018). The most massive stellar
clusters located near the galactic centre contain of the order of $10^{3}$
massive stars (Krumholz et al., 2019).

It should be noted that supernova explosions are still poorly understood. In
particular, simulations suggest that very massive stars (> 20-40 $M_{\odot}$)
may only produce a weak, sometimes called failed supernova Fryer (1999, 2003),
while stars with masses higher than 40 $M_{\odot}$ may not explode at all and
directly collapse onto a black hole. There is still no consensus nowadays
about the fate of massive stars. It may well be the case that only a narrow
mass interval allows for a powerful supernova explosion (Sukhbold et al.,
2016; Ebinger et al., 2020). Disregarding these debated issues, we assume in
the following that all massive stars, i.e. with mass higher than 8
$M_{\odot}$, explode in a powerful supernova. In a typical cluster
characterised by a Salpeter initial mass function, only 25% of the stars have
a mass higher than 20 $M_{\odot}$, so the controversial fate of these stars
does not affect significantly the mechanisms of CR acceleration.

Figure 1: Supernova explosion times as a function of stellar mass for a
cluster of 100 stars in the range 8-150 $M_{\odot}$.

Stars inject kinetic energy in the surrounding medium not only when they
explode as supernovae, but also during their entire lives, mainly due to the
presence of stellar winds (Cesarsky & Montmerle, 1983; Seo et al., 2018). A
cavity filled by hot and diluted gas, called SB, forms as the result of the
combined effect of stellar winds and supernova explosions. A simple
dimensional argument can be invoked to show that the radius of the bubble
increases with time as (Castor et al., 1975; Weaver et al., 1977; McCray &
Kafatos, 1987):

$\displaystyle R_{b}(t)$ $\displaystyle\approx
68~{}\xi_{b}^{1/5}\left(\frac{L_{*}}{10^{38}~{}{\rm
erg/s}}\right)^{1/5}\left(\frac{n_{\rm ISM}}{\rm
cm^{-3}}\right)^{-1/5}\left(\frac{t}{\rm Myr}\right)^{3/5}\rm pc\,,$ (1)

where $L_{*}$ is the total average mechanical power injected in the SB by
stellar winds and supernova explosions, $n_{\rm ISM}$ is the hydrogen number
density of the surrounding medium, and $\xi_{b}$ a numerical factor to be
determined from more accurate studies (e.g. Yadav et al., 2017). In the
expression above, physical quantities have been normalised to typical values,
and the implicit assumptions have been made that $R_{b}$ is much larger than
the radius of the star cluster and that SNRs become subsonic before reaching
the edge of the SB.

To estimate the wind mechanical power, denoted $\mathcal{P}_{w}$ in the
following, we account for both the main sequence (MS) and Wolf-Rayet (WR) wind
phases, while we neglect the red supergiant phase which gives a negligible
contribution (Seo et al., 2018). Stellar evolution models provide an estimate
of the wind power in either phase as (Seo et al., 2018):

$\displaystyle\begin{aligned}
\log_{10}\left(\frac{\mathcal{P}_{w,MS}(M)}{\mathrm{~{}erg~{}Myr}^{-1}}\right)&\approx-3.4\left(\log_{10}M\right)^{2}+15\log_{10}M+34\,,\\\
\mathcal{P}_{w,WR}(M)&\approx 6\times
10^{48}M^{1.2}\mathrm{~{}erg~{}Myr}^{-1}\,.\end{aligned}$ (2)

While the duration of the MS is computed from the fits provided by Limongi &
Chieffi (2006), the duration of the WR phase is assumed to be 320 kyr,
independently on the initial mass of the star provided it is higher than 20
$M_{\odot}$ (Seo et al., 2018).

As an example, the lifetime of a star of initial mass $M\sim 20~{}M_{\odot}$
is around 20 Myr, its total MS wind energy is of the order 1049 erg, and its
total wind energy in the WR phase is around 1050 erg. For the most massive O
stars, $M\sim 100~{}M_{\odot}$, the total wind energy reaches 1051 erg, which
is comparable to the mechanical energy of a supernova explosion.

Figure 2: Evolution of the mechanical power of the cluster assuming a
relaxation time $\tau_{T}=1$ Myr for the SNR perturbations (see Equations 2
and 7). “Massive” stars are stars with $M>8M_{\odot}$. The dash-dotted lines
indicate the contribution of the winds, $\mathcal{P}_{w}$. The dashed lines
indicate the average mechanical power.

The mechanical power of two randomly generated clusters containing
respectively 100 and 500 massive stars is plotted in Figure 2. In the early
stage of the cluster ($t\lesssim 5$ Myr), the contribution of the WR stars
dominates, while in older SBs the supernovae are the main sources of
mechanical power. Because the WR and SNR phases last for a short time, the
injection of mechanical energy in the bubble is intermittent, although the
mean value is roughly constant in time. The wind power quickly drops after 10
Myr, which is the time at which all stars of mass higher than $20~{}M_{\odot}$
have exploded. Then there only remain MS stars. Eventually, the average power
of supernovae over the lifetime of the cluster is $9\times 10^{35}$ erg/s/star
while the average power of the winds of massive ($M>8M_{\odot}$) stars is
$9\times 10^{34}$ erg/s/star, for a total average power of $10^{36}$
erg/s/star. The winds are therefore subdominant, yet non negligible, in
particular because the input from supernovae is strongly intermittent. The
total production of energy by a cluster of 100 massive stars is of the order
$10^{53}$ erg. Note that the average power injected into the SB,
$\int_{0}^{\mathcal{T}_{\rm SB}}\mathrm{d}t\mathcal{P}_{\rm
tot}(t)/\mathcal{T}_{\rm SB}$, where $\mathcal{T}_{\rm SB}$ is the SB
lifetime, is equal to the average luminosity (or mechanical power) of the
cluster $L_{*}$ which drives the dynamics of the bubble. Eq. 1 can therefore
be rewritten as:

$R_{b}(t)\approx R_{0}N_{*}^{1/5}\left(\frac{t}{\rm Myr}\right)^{3/5}\,,$ (3)

where $N_{*}$ is the number of massive ($M>8M_{\odot}$) stars in the cluster
and $R_{0}=27~{}\xi^{1/5}\left(n_{\rm ISM}/\rm cm^{-3}\right)^{-1/5}\rm pc$.
Assigning a numerical value to $R_{0}$ is sufficient to determine the entire
dynamical evolution of the SB.

Figure 3 displays the radius of several known SBs as function of their age.
These observations suggest the scaling $n_{\rm ISM}/\xi_{b}\approx 450$ cm-3.
In the following we will assume $n_{\rm ISM}=100$ cm-3, which is a typical
density of giant molecular clouds where young clusters are expected to be born
and expand (Parizot et al., 2004). This implies $\xi_{b}\approx 22$% and
$R_{0}\approx 10$ pc.

Figure 3: Sizes and ages of several observed clusters. The error bars reflect
the variations found in the literature rather than the physical uncertainties.
In particular, the number of massive stars in the G25 region is only roughly
estimated, and Westerlund 2 is poorly constrained. Data are found in Ackermann
et al. (2011, and references therein) (Cygnus); Joubaud et al. (2019, and
references therein),Voss et al. (2010) (Orion-Eridani); Lopez et al. (2020,
and references therein) (30 Doradus C); Katsuta et al. (2017, and references
therein) (G25); Maddox et al. (2009) (N11); Dunne et al. (2003) (Omega
nebula); Smith (2006); Smith & Brooks (2007) ($\eta$-Carina); Rauw et al.
(2007) (Westerlund 2); Cooper et al. (2004) (DEM L192); Jaskot et al. (2011,
and references therein) (DEM L152 and DEM L50); Dunne et al. (2001); Massey et
al. (1995) (DEM 199). See also the compilation in Ferrand & Marcowith (2010,
and references therein).

As the SB expands, the medium is diluted and cooled. The seminal SB model of
Weaver et al. (1977) provides the following self-similar scalings for the
average density and temperature inside the bubble (Mac Low & McCray, 1988):

$\displaystyle n(t)$
$\displaystyle=0.34\left(\frac{\xi_{b}N_{*}}{100}\right)^{\frac{6}{35}}\left(\frac{n_{\rm
ISM}}{100~{}{\rm cm}^{-3}}\right)^{\frac{19}{35}}\left(\frac{t}{{\rm
Myr}}\right)^{-\frac{22}{35}}~{}{\rm cm}^{-3}\,,$ (4) $\displaystyle T(t)$
$\displaystyle=4.8\times
10^{6}\left(\frac{\xi_{b}N_{*}}{100}\right)^{\frac{8}{35}}\left(\frac{n_{\rm
ISM}}{100~{}{\rm cm}^{-3}}\right)^{\frac{2}{35}}\left(\frac{t}{{\rm
Myr}}\right)^{-\frac{6}{35}}~{}{\rm K}\,.$ (5)

Together with the expansion of the outer shock, the evolution of the density
and the temperature describes the dynamical SB environment.

###  2.2 Turbulence generation

Inside the SB, the stellar mechanical energy is converted into turbulence
which cascades from the largest scale to the dissipation scale. We assume
equipartition of the kinetic and magnetic energies in the waves and denote $W$
the spectral energy density such that $\int\mathrm{d}kW(k)=\delta B^{2}/4\pi$.
The dynamics of the turbulence can be described by the following nonlinear
local energy transfer equation (Zhou & Matthaeus, 1990; Norman & Ferrara,
1996; Miller & Roberts, 1995; Ptuskin et al., 2005):

$\partial_{t}W+\partial_{k}\left(\frac{a}{\sqrt{\rho}}k^{5/2}W^{3/2}\right)=-\Gamma(k)W+S\delta(k-k_{0})\,,$
(6)

where $a\approx 0.8$ is a numerical constant determined from experiments or
simulations (e.g. Verma et al., 1996), $\rho$ is the gas density, $2\pi/k_{0}$
is the largest turbulent scale, $S$ is the injection of energy at the largest
scale and $\Gamma$ is the damping rate.

We assume that both winds and supernovae convert a fraction $\eta_{T}$ of
their mechanical energy into magnetised waves, and that the conversion is
exponentially suppressed for SNRs, such that the source term of the turbulence
cascade is phenomenologically written as:

$\displaystyle\begin{aligned} S&=\frac{\eta_{T}}{V_{\rm
SB}(t)}\mathcal{P}_{\rm tot}\,,\\\ \mathcal{P}_{\rm
tot}&=\sum_{i=1}^{N(t)}\mathcal{P}_{w}(M_{i})+\sum_{N(t)+1}^{N_{*}}\frac{E_{\rm
SN}}{\tau_{T}}e^{-(t-t_{i})/\tau_{T}}\,,\end{aligned}$ (7)

where $N(t)$ is the number of stars remaining in the cluster at time $t$,
$t_{i}$ is the time at which the star number $i$ is expected to explode,
$E_{\rm SN}=10^{51}$ erg is the energy released per supernova and $\tau_{T}$
is the relaxation time of the turbulence after a supernova explosion. This
recipe is similar to that adopted in Fang et al. (2019) for the conversion of
the wind power. The exponential decay of the turbulence excited by SNRs is
motivated by numerical simulations (e.g. Krause et al., 2013), which suggest a
relaxation time $\tau_{T}\approx 1$ Myr.

In the stationary regime, the general solution of Eq. 6 reads:

$\displaystyle W(k)$
$\displaystyle=k^{-5/3}\left(\frac{\sqrt{\rho}S}{a}\right)^{2/3}\left(1-\frac{1}{3}\left(\frac{\rho}{Sa^{2}}\right)^{1/3}\int_{k_{0}}^{k}\mathrm{d}k^{\prime}\frac{\Gamma(k^{\prime})}{k^{\prime
5/3}}\right)^{2}\,.$ (8)

In the absence of damping, one retrieves the Kolmogorov scaling $W\propto
k^{-5/3}$. The generalisation to another turbulence regime (e.g. Kraichnan) is
straightforward. The test-particle solution provides an estimate of the energy
density contained in the turbulence:

$\frac{\delta B^{2}}{4\pi}\approx\rho\delta
u^{2}\approx\frac{3}{2}\left(\frac{\sqrt{\rho}S}{ak_{0}}\right)^{2/3}\,.$ (9)

Considering turbulence strengths above this limit would violate energy
conservation in the cluster. Eq. 9 allows to estimate the turbulent magnetic
field $\delta B$ and velocity $\delta u$:

$\displaystyle\delta B$ $\displaystyle\lesssim 5\left(\frac{n}{0.01~{}{\rm
cm}^{-3}}\right)^{\frac{1}{6}}\left(\frac{\eta_{T}\mathcal{P}_{\rm
tot}}{10^{51}~{}{\rm erg/Myr}}\frac{\lambda}{10~{}{\rm
pc}}\right)^{\frac{1}{3}}\left(\frac{R_{b}}{50~{}{\rm
pc}}\right)^{-1}~{}\text{\textmu G}\,,$ (10) $\displaystyle\delta u$
$\displaystyle\lesssim 110\left(\frac{n}{0.01~{}{\rm
cm}^{-3}}\right)^{-\frac{1}{3}}\left(\frac{\eta_{T}\mathcal{P}_{\rm
tot}}{10^{51}~{}{\rm erg/Myr}}\frac{\lambda}{10~{}{\rm
pc}}\right)^{\frac{1}{3}}\left(\frac{R_{b}}{50~{}{\rm
pc}}\right)^{-1}~{}\text{km/s}\,,$ (11)

where $\lambda=2\pi/k_{0}$ is the largest turbulent scale.

The density, power and radius of the SB depend on the total number of stars
$N_{*}$ as well as the time $t$. From the scalings derived above ($n\propto
N_{*}^{6/35}t^{-22/35}$; $\mathcal{P}_{\rm tot}\propto N_{*}$; $R_{b}\propto
N_{*}^{1/5}t^{3/5}$), we get to conclude that the random component of the
velocity $\delta u$ is weakly dependent on the total number of stars: $\delta
u\propto N_{*}^{8/105}$, while it decreases with time as $\delta u\propto
t^{-41/105}$. The sound speed being of the order of 100 km/s inside the
superbubble, the turbulence is expected to consist in an ensemble of weak
secondary shock waves in the early times of the cluster history, namely the
first few million years. The transport of charged particles in such strong
supersonic turbulence has been investigated by Bykov & Toptygin (1990, 1993),
who later showed that the acceleration of the particles could be very
efficient, with up to 30 % of the energy of the waves converted into CRs and
such that hard CR spectra should be expected up to TeV energies (Bykov &
Fleishman, 1992a, b; Bykov, 2001), while a steep component could remain up to
energies as high as 100 PeV (Bykov et al., 1995). In the following we shall
rather analyse the acceleration of particles in more evolved SBs, in which the
turbulence is expected to be subsonic. One should nevertheless keep in mind
that the early phase of efficient particle acceleration in supersonic
turbulence could last several million years and be an efficient CR accelerator
(e.g. Bykov, 2001).

Figure 4: Evolution of the turbulent magnetic field and turbulence level in a
typical SB with $N_{*}=100$, $B_{0}=10$ µG, $\eta_{T}=30\%$, as estimated from
Eq. 9.

The evolution of the random component of the magnetic field is plotted in
Figure 4 for a cluster containing 100 massive stars, where we assumed a
background field $B_{0}=10$ µG and an efficiency of turbulence generation
$\eta_{T}=30\%$. The quasi-linear theory is, strictly speaking, only valid
when $\delta B<B_{0}$, or $\eta<50\%$, where $\eta=\delta
B^{2}/(B_{0}^{2}+\delta B^{2})$ is the turbulence level. This is verified
after a few million years for a background field $B_{0}=10$ µG. On the other
hand, the turbulence may remain strong until the end of the cluster life if
e.g. $B_{0}=1$ µG (still assuming $\eta_{T}=30\%$).

The turbulence level decreases rapidly in the first million years as WR stars
explode and the wind power diminishes. Then it stabilises after about 10 Myr
around a few percent (for $B_{0}=10$ µG), when only MS stars remain and the
main sources of turbulence are the SN explosions.

##  3 Particle acceleration and transport in superbubbles

###  3.1 Particle acceleration at stellar wind termination shocks

The stellar outflows may form a collective wind termination shock (WTS) around
a star cluster. The position of the WTS $R_{s}$ is determined by the
equilibrium between the ram pressure of the collective wind and the pressure
in the interior of the SB, which provides (e.g. Gupta et al., 2020; Morlino et
al., 2021):

$R_{s}\approx\frac{(\xi_{b}L_{*})^{-1/5}\mathcal{P}_{w}^{1/2}}{\left(10^{36}~{}{\rm
erg/s}\right)^{3/10}}\left(\frac{n_{\rm ISM}}{100~{}{\rm
cm}^{-3}}\right)^{-3/10}\left(\frac{V_{w}}{2000~{}{\rm km/s}}\right)^{-1/2}\\\
\times\left(\frac{t}{1~{}{\rm Myr}}\right)^{2/5}~{}{\rm pc}\,,$ (12)

where $L_{*}\approx 10^{36}N_{*}$ erg/s is the average mechanical power
imparted by the winds and the supernovae (see the dashed lines in Fig. 2),
which determines the pressure in the SB interior, while $\mathcal{P}_{w}$ is
the mechanical power of the winds, which determines the ram pressure of the
WTS. Eq. 12 is only valid when the wind power is independent of time. Yet it
provides a reliable estimate of the radius of the termination shock in our
setup, where the wind power is constant most of the time and varies only
significantly during the brief WR phase of the stars just before the SN
explosions. In the early phase of the cluster life, the WTS radius may reach
about 10 pc. On the other hand, after 10 Myr only MS stars remain and the wind
power is subdominant, such that the collective WTS can barely be sustained and
shrinks to sub-parsec sizes. For simplicity, we will therefore consider loose
clusters in the following, assuming that the extent of the cluster (or that of
the region of energy deposition) is of the order of 10 pc. In this case, a
collective WTS cannot form and the winds produce isolated strong shocks around
each stars, with radii determined by Eq. 12 provided the total wind power
$\mathcal{P}_{w}$ is replaced by the power of the individual stars.

A fraction of the mechanical energy of the winds can be converted into CRs via
DSA operating at the individual WTS (Cesarsky & Montmerle, 1983; Seo et al.,
2018). Because winds do not produce a high energy density of particles, we
assume that this acceleration process takes place in the test-particle regime.
We can therefore describe the acceleration of particles at the WTS as a
continuous and time dependent injection of CRs characterised by an injection
rate of particles of momentum $p$ equal to:

$q_{w}(p)=\frac{\eta_{w}p^{-4}e^{-p/p_{w}}}{\mathcal{N}_{w}}\sum_{i=1}^{N(t)}\mathcal{P}_{w}(M_{i})\,,$
(13)

with $\eta_{w}\approx 10\%$ the fraction of the wind mechanical energy that
goes into CRs, $\mathcal{P}_{w}(M_{i})$ the mechanical power of the wind
ejected by the $i$th star in the cluster, $p_{w}$ the maximum momentum of the
particles, and $N(t)$ the number of stars remaining in the cluster (i.e. not
yet exploded) at time $t$. The maximum momentum achievable in a nearly
stationary WTS is expected to be limited by the size of the shock rather than
its lifetime: $p_{w}\sim ZeBV_{w}R_{s}\sim 10^{5}$ GeV (e.g. Voelk & Biermann,
1988; Morlino et al., 2021). Finally, $\mathcal{N}_{w}$ is a normalisation
constant given by the following integral:

$\displaystyle\mathcal{N}_{w}=4\pi\int_{p_{0}}^{p_{w}}\mathrm{d}p\,\epsilon(p)p^{-2}e^{-p/p_{w}},$
(14)

where $p_{0}$ is the injection momentum and $\epsilon(p)$ the kinetic energy
of protons of momentum $p$.

###  3.2 Particle reacceleration at supernova remnant shocks

When a star explodes into a supernova, a SNR is formed. The SNR not only
injects fresh CRs in the stellar cluster, but it also reaccelerates a fraction
of the already existing particle distribution. The lifetime of a SNR in the
medium of density $0.01-0.1$ cm-3 is of the order $10-100$ kyr, which is short
compared to the CR diffusion time away from the cluster (Myr), hence we assume
that the SNR are instantaneous and compute CR propagation only between them
(as in Ferrand & Marcowith, 2010). The CR spectrum after the passage of the
SNR, $n_{f}(p)$, is written as:

$n_{f}=\mathcal{V}_{\rm SNR}\,\mathcal{R}\left[\frac{n_{i}}{V_{\rm
SB}}\right]+\left(1-\frac{\mathcal{V}_{\rm SNR}}{V_{\rm SB}}\right)n_{i}\,,$
(15)

where $n_{i}(p)$ is the spectrum of CRs preexisting before the supernova
explosion, $\mathcal{V}_{\rm SNR}$ is the volume spanned by the SNR, $V_{\rm
SB}=4\pi/3R_{b}^{3}$ is the SB volume and $\mathcal{R}$ the “reacceleration
operator”. In the test-particle limit, this operator simply reads (e.g.
Melrose & Pope, 1993; Ferrand & Marcowith, 2010):

$\mathcal{R}[f]=4\int_{p_{0}}^{r^{\frac{1}{3}}p}\frac{\mathrm{d}p^{\prime}}{p^{\prime}}\left(\frac{r^{\frac{1}{3}}p}{p^{\prime}}\right)^{-4}\left(f(p^{\prime})+\eta_{\rm
SNR}\delta\left(p^{\prime}-p_{0}\right)\right)\,,$ (16)

which describes both the injection of particles from the thermal bath with an
efficiency $\eta_{\rm SNR}$ at the injection momentum $p_{0}$ and the
reacceleration of the preexisting distribution of particles. However, the
test-particle reacceleration operator cannot be applied in environments where
multiple shocks are expected over small timescales, because this would
systematically lead to the violation of the energy conservation as most SNR
would transfer more than $10^{51}$ erg into non-thermal particles (see Vieu et
al., 2021, for a detailed discussion of this issue). The linear reacceleration
operator must therefore be substituted by a nonlinear, self-consistent,
computation. In the following we use the procedure described in Vieu et al.
(2021). In order to take into account the feedback of the CR pressure on the
hydrodynamic profile and small-scale turbulence, two differential equations
must be solved together with suitable boundary conditions. On one hand, the
shock structure is modified by the non-thermal particles according to the
conservation of momentum, while on the other hand the acceleration of the
particles is driven by the transport equation integrated around the shock. The
injection efficiency $\eta$ can be related to the injection momentum by
imposing the continuity between the thermal and non thermal distributions of
particles at the injection momentum (Blasi et al., 2005). Finally, adiabatic
decompression is computed after the explosion of each supernova. As shown in
Vieu et al. (2021), the spectrum of reaccelerated particles tends towards a
universal shape after a few successive reaccelerations. Thus, this spectrum is
expected to be largely independent of the injection parameter $\xi$, which is
defined such that the injection momentum equals $\xi p_{\rm th}$, where
$p_{\rm th}$ is the typical momentum of the thermal particles (see Vieu et
al., 2021). In the following, we will take $\xi=3.5$, which is chosen such
that, in the absence of preexisting particles ahead of the shock, the remnant
transfers about $10^{50}$ erg into non-thermal CRs. Assuming that SNR are
strong shocks, the only remaining parameter is the maximum momentum. We
constrain the latter in a simple and conservative approach, without assuming
field amplification far upstream of the SNR.The most stringent limitation on
the maximum energy is therefore the acceleration time. The latter depends on
the dynamics of the SNR, which is not straightforward to compute when the
shock expands in an inhomogeneous medium. Disregarding these complications, we
assume that the velocity is constant in the free expansion phase, which is a
good approximation especially when the remnant expands in a wind density
profile (e.g. Finke & Dermer, 2012; Gaggero et al., 2018). SNRs expanding in a
SB spend a large fraction of their life into the free expansion phase, at the
end of which the mechanical power of the shock is maximal. The maximum
momentum is also expected to be achieved at the end of the free expansion
phase and can be estimated as:

$\displaystyle p_{\rm max}=1.4\frac{ZB}{10~{}\text{\textmu
G}}\left(\frac{V_{\rm
SN}}{3000~{}\text{km/s}}\right)\left(\frac{M_{e}}{10~{}M_{\odot}}\right)^{\frac{1}{3}}\left(\frac{n}{0.01{{\rm
cm}^{-3}}}\right)^{-\frac{1}{3}}\text{~{}PeV}\,,$ (17)

where $n$ and $B$ are respectively the density and the total magnetic field in
the SB interior and $Z$ is the atomic number of the particle. Although this is
admittedly a rough estimate, it allows to get rid of an otherwise free
parameter without solving the time-dependent nonlinear problem of shock
reacceleration accounting for the multiple environments the remnant
experiences before entering the adiabatic phase. Besides, it contains all the
relevant physical ingredients. The maximum momentum is indeed expected to
increase with enhanced magnetic field and in low density media.

As shown in Vieu et al. (2021), it is expected that a substantial fraction of
the non-thermal energy leaks upstream of the shock, carried by high-energy
particles beyond the maximum momentum. Because these particles will quickly
escape from the SB (within typically 10 kyr), they are not expected to affect
the overall CR content and we disregard them in the present analysis. However,
they should interestingly provide an intermittent input of very high energy
particles in the Galaxy, leaving a signature in the region of the knee which
is observed in the galactic cosmic ray spectrum. This additional component can
however only be treated accurately using a time-dependent nonlinear model of
CR reacceleration at SNR shocks, which is beyond the framework established in
Vieu et al. (2021) and will be considered in a future work.

The only remaining quantity to constrain is the volume spanned by the SNR,
$\mathcal{V}_{\rm SNR}$, which determines the fraction of CRs reaccelerated by
the expanding shock (Eq. 15). In the hot rarefied medium characterising the SB
interior, a SNR becomes subsonic before becoming radiative (Mac Low & McCray,
1988). In principle, $\mathcal{V}_{\rm SNR}$ should therefore correspond to
the volume spanned by the SNR before the forward shock weakens. However,
because DSA is solved in the stationary regime the Mach number at which the
shock “weakens” is not well defined. Instead of keeping it as a free
parameter, we rather impose that the kinetic energy of the stationary strong
shock multiplied by the volume $\mathcal{V}_{\rm SNR}$ equals the energy of
the supernova explosion:

$\rho\mathcal{V}_{\rm SNR}u_{0}^{2}/2=10^{51}~{}{\rm erg}\,,$ (18)

which ensures that the energy of the stationary shock will never exceed the
energy budget of the time-dependent SNR. The fraction of reaccelerated CRs
thus reads:

$\frac{\mathcal{V}_{\rm SNR}}{V_{\rm
SB}}\approx\min\left(0.59\left(\frac{R_{b}(t)}{20~{}\mathrm{pc}}\right)^{-3}\left(\frac{n(t)}{0.1~{}\mathrm{cm}^{-3}}\right)^{-1}\left(\frac{u_{0}}{1500~{}\mathrm{km/s}}\right)^{-2},1\right)\,,$
(19)

where $t$ is the age of the SB in Myr and $u_{0}\approx 1500$ km/s is the
velocity of the SNR shock at the beginning of the Sedov-Taylor phase, which is
independent of the density and thus stays nearly constant throughout the SB
evolution (Parizot et al., 2004). Eq. 18 could have been equivalently obtained
by imposing that the acceleration process stops when the Mach number of the
SNR shock falls below $\mathcal{M}_{\rm min}\approx
2.5\left(\frac{u_{0}}{1500~{}\mathrm{km/s}}\right)^{7/2}\left(\frac{T}{10^{6}~{}\rm
K}\right)^{-1/2}$, which is consistent with the intuitive definition of a
“weak” shock. Assuming $N_{*}=100$ and $n_{\rm ISM}=100$ cm-3, and using the
scaling described above, Eq. 19 becomes: $\mathcal{V}_{\rm SNR}/V_{\rm
SB}\approx 0.16(t/{\rm Myr})^{-1.17}$, which means that only a few percent of
the CR content will be reaccelerated by each SNR.

To summarise, in order to reduce an expanding supernova to a stationary shock,
we assume that (i) the stationary shock velocity is equal to the velocity of
the supernova at the beginning of the Sedov-Taylor phase ($u_{0}=1500$ km/s),
(ii) the stationary shock disappears when its integrated kinetic energy
$\rho\mathcal{V}_{\rm SNR}u_{0}^{2}/2$ is equal to the energy of the supernova
explosion ($10^{51}$ erg), (iii) particles are injected from the thermal bath
with a stationary efficiency such that in the absence of preexisting particles
$10^{50}$ erg are transferred into CRs ($\xi=3.5$), (iv) the time-independent
maximum energy achieved by the freshly injected particles as well as that
beyond which reaccelerated particles escape upstream is about that reached at
the end of the free expansion phase111In principle the energy of the
reaccelerated particles is not limited by the acceleration time but by the
finite size of the shock. However both criteria provide very similar estimates
at the beginning of the Sedov-Taylor phase., (v) we disregard the upstream
flux of very high energy particles assuming that they quickly leave the SB.

###  3.3 Stochastic acceleration

When they propagate inside the SB, particles experience resonant scattering
against the magnetised waves. Each scattering deflects the particles in such a
way that on average, an effective diffusion takes place both in space and
momentum. The momentum diffusion operator reads:

$\partial_{t}f=\frac{1}{p^{2}}\partial_{p}\left(p^{2}D_{p}\partial_{p}f(p)\right)\,,$
(20)

where $f(p)$ is the distribution function averaged over the spherical
coordinates $\mu$ and $\phi$:

$f(p)=\frac{1}{4\pi}\int_{0}^{2\pi}\mathrm{d}\phi\int_{-1}^{1}\mathrm{d}\mu\,f(p,\mu,\phi)\,,$
(21)

with $\mu$ the pitch-angle cosine and $\phi$ the gyrophase, and $D_{p}$ is the
pitch-angle averaged diffusion coefficient:
$D_{p}(p)=1/2\int_{-1}^{1}\mathrm{d}\mu D_{p}(p,\mu)$. These normalisations
ensure that the number density of the particles reads
$n=4\pi\int\mathrm{d}pp^{2}f(p)$.

In the high energy regime $c\gg v_{A}$, with $v_{A}$ the Alfven speed, the
quasi-linear theory of particle diffusion in turbulence provides the following
expression for the pitch-angle averaged diffusion coefficient $D_{p}$
(Berezinskii et al., 1990):

$\displaystyle D_{p}$
$\displaystyle=\frac{\pi^{2}Z^{2}e^{2}v_{A}^{2}}{v}\int_{0}^{1}\frac{\mathrm{d}\mu}{\mu}\left(1-\mu^{2}\right)W\left(\frac{ZeB}{\mu
p}\right)\,,$ (22)
$\displaystyle=\frac{\pi^{2}Z^{2}e^{2}v_{A}^{2}}{v}\int_{ZeB/p}^{\infty}\frac{\mathrm{d}k}{k}\left(1-\left(\frac{ZeB}{kp}\right)^{2}\right)W\left(k\right)\,,$
(23)

where $W$ describes both co- and counter-propagating waves of equal energy,
$Ze$ is the charge of the particle and $v$ its velocity. Here and in the
following the background magnetic field $B_{0}$ is systematically supplemented
by the perturbation $\delta B$ computed according to Eq. 10 and we write the
total field $B=B_{0}+\delta B$. This allows to consistently account for the
case of strong turbulence, $\delta B\gtrsim B_{0}$, where the quasi-linear
theory only makes sense locally, for particles follow the mean resulting
field. An accurate description of the strong turbulence regime is beyond the
scope of this work. We recall that for the parameters considered in this
analysis, the turbulence is strong only in the early phase of the SB
evolution.

Once again, in the case where the energy of the diffusing particles is similar
to that of the turbulence, one must account for the backreaction of the
particles on the turbulence, which is expected to damp the waves. To include
this nonlinear effect, we follow Ptuskin et al. (2003, 2005, 2006); Ptuskin et
al. (2017). Assuming an isotropic distribution of particles, the damping rate
can be estimated through the following flux conservation (Eilek, 1979;
Thornbury & Drury, 2014):

$\displaystyle\int\mathrm{d}k\Gamma(k)W(k)$
$\displaystyle=4\pi\int\mathrm{d}p\,p^{2}\frac{\epsilon(p)}{p^{2}}\partial_{p}(p^{2}D_{p}\partial_{p}f)\,,$
(24) $\displaystyle=4\pi\int\mathrm{d}p\,\partial_{p}(vp^{2}D_{p})f(p)\,.$
(25)

The quantity $\partial_{p}(vp^{2}D_{p})$ is the mean energy transferred per
unit time from the waves to the particles. Using Eq. 22, we work out the
right-hand side to identify:

$\Gamma(k)=\frac{8\pi^{3}Z^{2}e^{2}v_{A}^{2}}{k}\int_{ZeB/k}\mathrm{d}p\,pf(p)\,.$
(26)

This expression is identical to that written in e.g. Eilek (1979); Miller &
Roberts (1995); Brunetti et al. (2004). It differs from the expression
appearing in the works of Ptuskin et al. (e.g. 2006) by a factor $2\pi$ due to
the choice of normalisation of Eq. 6.

We assume that the turbulence cascade is faster than the other physical
processes occurring in the SB such that we can neglect the time-dependency of
the turbulence spectrum. We will comment a posteriori on the validity of this
hypothesis. From the stationary solution given by Eq. 8, and assuming that the
turbulence spectrum is steep enough for Eq. 22 to be approximated within a
factor of order unity as $D_{p}\approx
3\pi^{2}v_{A}^{2}Z^{2}e^{2}W(ZeB/p)/5v$, we obtain:

$D_{p}\approx\frac{3\pi
k_{0}^{1/3}p^{5/3}}{20v}\left(\frac{ZeB}{k_{0}}\right)^{1/3}\left(\frac{S}{a\rho}\right)^{2/3}\\\
\left(1-\frac{2\pi^{2}}{5}\left(\frac{ZeB}{\rho^{2}a^{2}S}\right)^{1/3}\int_{p}\mathrm{d}p^{\prime}f(p^{\prime})p^{\prime
8/3}\right)^{2}\,.$ (27)

The feedback of the particles on the turbulence is non trivial. Not only it
reduces the density of the magnetic waves, but it can also terminate the
turbulence cascade at a finite scale. Indeed, it may happen that for high
enough CR energy densities the quantity in the parenthesis in Eq. 27 vanishes
at small momenta. Low energy particles will not be able to resonate anymore
with the magnetised waves, and their momentum diffusion coefficient should
therefore be set to zero. This was already pointed out by Ptuskin et al.
(2006) in the context of CR diffusion in the interstellar medium.

Noticing that at $p=p_{0}$ the integral in the RHS of Eq. 27 is nearly the CR
energy density, the momentum diffusion coefficient will vanish at momenta
$p>p_{0}$ if the energy density of the particles $e_{\rm CR}$ satisfies:

$\displaystyle e_{\rm CR}$ $\displaystyle\gtrsim
3c\left(\frac{\rho^{2}Sp_{0}}{ZeB}\right)^{1/3}\,,$ (28)

i.e. if the CR energy density is higher than typically 10 eV cm-3 for standard
SB environments. This estimate is consistent with that obtained by Ptuskin et
al. (2006) who concluded that the feedback of the particles on the turbulence
spectrum of the ISM should be important below GeV energies. In young SBs, an
energy density of 1-10 eV/cm3 is already reached after the very first
supernova explosion: the nonlinear CR feedback on the turbulence cascade must
therefore be taken into account even for young clusters.

To obtain the turbulence spectrum, and thus the momentum diffusion
coefficient, we solved the stationary version of Eq. 6. The timescale of the
energy transfer at the scale $2\pi/k$ is $\tau=\rho/(a\sqrt{Wk^{3}})$. For a
Kolmogorov spectrum, this timescale is always much smaller than the scattering
time of the particles on the magnetic waves and the hypothesis of stationarity
is well justified.

Together with stationarity, other simplifications have been done in the above
derivation. We assumed that all modes of the turbulence followed a Kolmogorov-
type cascade, whereas in realistic cosmic environments two cascades may
coexist, one describing the dynamics of Alfvén waves and the other describing
the fast magnetosonic modes (Ptuskin et al., 2003).

###  3.4 Spatial diffusion and escape

When they propagate in the turbulent medium, not only particles are
reaccelerated by stochastic scatterings on magnetic waves, but they also
experience an effective spatial diffusion until they escape in the
interstellar medium (ISM). Under the quasilinear hypothesis, the following
approximate relation between the spatial and momentum diffusion coefficients
is obtained for a Kolmogorov turbulence222The “standard” relation
$D_{x}D_{p}=v_{A}^{2}p^{2}/9$ is only valid in Bohm’s regime (Thornbury &
Drury, 2014).:

$D_{x}D_{p}\approx 0.2v_{A}^{2}p^{2}\,.$ (29)

This relation, derived for a power law turbulence spectrum $W(k)\propto
k^{-5/3}$, is also valid if small scales are damped, providing the largest
scales do follow a Kolmogorov scaling. In the test-particle regime, Eq. 27
eventually provides the following estimate:

$D_{x}\approx
10^{28}\beta\left(\frac{p}{m_{p}c}\right)^{\frac{1}{3}}\left(\frac{B}{10~{}\text{\textmu
G}}\right)^{\frac{5}{3}}\left(\frac{\eta_{T}L_{*}}{10^{51}~{}\text{erg/Myr}}\right)^{-\frac{2}{3}}\\\
\times\left(\frac{n}{0.01~{}\text{cm}^{-3}}\right)^{-\frac{1}{3}}\left(\frac{R_{b}}{100~{}\text{pc}}\right)^{2}~{}\text{cm}^{2}/\text{s}\,.$
(30)

In the case where the CR energy density is so high that the turbulence cascade
terminates at a large scale such that $D_{p}$ vanishes (this happens when the
quantity in parenthesis in Eq. 27 vanishes), the above expression breaks down.
In this case, low energy particles do not scatter on waves anymore, but rather
follow the field lines of the background field. In a highly turbulent
medium333In this case, the quasi-linear approximation used to compute the
diffusion of high energy particles on the waves is understood as a local
approximation., the background field has a coherence length of the order of
$0.77/k_{0}$ (Casse et al., 2001). Low energy particles following the
background field do experience a random walk with mean displacement
$0.77/k_{0}$. From a macroscopic point of view, this gives rise to an
effective spatial diffusion with diffusion coefficient about $0.3v/k_{0}$. We
conveniently interpolate the coefficient between the low and high energy
regimes as:

$D_{x}\approx\left(\frac{5D_{p}}{v_{A}^{2}p^{2}}+\frac{3k_{0}}{v}\right)^{-1}\,.$
(31)

The quasilinear relation given by Eq. 29 is retrieved whenever the particles
scatter on magnetic waves, while in the absence of waves one gets
$D_{x}\approx 0.3v/k_{0}$.

In fact, the diffusion of CRs below the inertial turbulence cascade is a more
complicated process than the simple picture presented above, in particular
because the anisotropic streaming along the background field is expected to
trigger the so-called streaming instability, such that the particles will
scatter on the waves they produce. Accounting for the streaming instability
requires the knowledge of the spatial gradient of the distribution function
along the magnetic field lines, which is beyond the scope of this work.
Furthermore, it is not straightforward to solve the turbulence cascade
including the growth of the streaming instability. For simplicity, we assume
that the background field is chaotic enough to efficiently isotropise the low
energy particles, such that the streaming instability can be neglected.
Moreover, although it is crucial to account for the feedback of the energetic
particles on the hydromagnetic waves, the cascade is in fact very scarcely
suppressed for the parameters we consider.

We eventually average the CR transport equation over the SB volume such that
the spatial diffusion reduces to a leakage with characteristic escape time
$\tau=R_{b}^{2}/(6D_{x})$. The complete averaged CR transport equation in the
SB then reads:

$\partial_{t}n=-\frac{n}{\tau}+\frac{1}{p^{2}}\partial_{p}\left(p^{2}D_{p}\partial_{p}n\right)+\frac{1}{p^{2}}\partial_{p}\left(p^{2}\left[\left.\frac{\mathrm{d}p}{\mathrm{d}t}\right|_{\rm
Ad.}+\left.\frac{\mathrm{d}p}{\mathrm{d}t}\right|_{\rm
Int.}\right]n\right)+q_{w}\,,$ (32)

where $n=V_{\rm SB}f(p)$ is the CR spectrum in the SB. The first loss rate in
the right-hand side accounts for the adiabatic decompression of the CR
spectrum in the expanding medium,
$\left.\frac{\mathrm{d}p}{\mathrm{d}t}\right|_{\rm
Ad.}=(\epsilon(p)/v)\dot{R_{b}}/R_{b}$ (e.g. Finke & Dermer 2012 or Appendix B
of Aharonian 2004). The second loss rate
$\left.\frac{\mathrm{d}p}{\mathrm{d}t}\right|_{\rm Int.}$ encodes the losses
due to the Coulomb interactions and to the production of pions (Mannheim &
Schlickeiser, 1994).

The computation of the CR dynamics in the SB follows the procedure used by
Ferrand & Marcowith (2010). In between supernova explosions, the spectrum
evolves according to Eq. 32 implemented within a forward Euler scheme (a
Crank-Nicholson scheme can not be used since the momentum diffusion is
nonlinear). Whenever a supernova explodes, the reacceleration of the particles
is computed according to the formalism detailed in (Vieu et al., 2021). All
other physical processes are neglected during the supernova explosion.

The parameters used to model the fiducial SB considered in this work are
summarised in Table 1.

Initial number of massive ($M>8M_{\odot}$) stars | $N_{*}$ | $100-1000$
---|---|---
Superbubble radius normalised per star at 1 Myr | $R_{0}$ | 10 pc
ISM density around the SB | $n_{\rm ISM}$ |  $100$ cm-3
Large scale magnetic field | $B_{0}$ | 10 µG
Efficiency of turbulence generation | $\eta_{T}$ |  $1-30$ %
Turbulence relaxation time after SNR | $\tau_{T}$ | 1 Myr
Largest turbulent scale | $\lambda$ | 10 pc
Injection parameter of SNR shocks | $\xi$ | 3.5
Wind acceleration efficiency | $\eta_{w}$ | 10 %
Injection momentum at winds | $p_{0}$ | 5 MeV/c
Maximum momentum at winds | $p_{w}$ | 105 GeV/c
Table 1: Exhaustive list of the parameters characterising the fiducial
superbubble considered in this study.

##  4 Results

We now apply the model to the acceleration of protons in order to investigate
their energetics and spectra in SB.

###  4.1 Timescales

Figure 5: Typical timescales of the fiducial superbubble (see the parameters
in Table 1) for $N_{*}=100$, $\eta_{T}=30\%$, and a cluster age of 6 Myr. The
solid lines correspond to the nonlinear computation accounting for the CR
feedback on the turbulence spectrum, to be compared with the dotted lines
displaying the test-particle expectation. The horizontal black line indicates
the typical time interval between two supernovae in a cluster of 100 massive
stars (0.37 Myr).

Figure 5 displays the timescales of stochastic (Fermi II) acceleration
$p^{2}/(2D_{p})$ as well as the escape time $\tau$ for the fiducial SB. At
high energies ($p\gtrsim 10$ GeV/c), the standard Kolmogorov scalings are
always retrieved. Below 10 GeV, the cascade can be damped if there are enough
massive stars in the cluster. In this case, the stochastic acceleration is
less efficient than expected in the test-particle limit, while the escape is
faster. Sub-GeV energies are driven by the stochastic reacceleration while the
escape dominates above GeV energies. If the number of massive stars were
increased, supernovae would be more frequent and the horizontal black line in
Figure 5 would be lowered. Supernovae would then be frequent enough to sustain
the other processes in a broader energy range, between about 0.1 and 10 GeV,
leading to a nearly stationary production of CR in this range. On the other
hand, whenever the dominant process is faster than the SNR rate, the
production of CR is expected to be intermittent. This is the case at low and
high energies.

For completeness we added the timescale of the proton losses, including
Coulomb collisions, protons-protons interactions, as well as adiabatic
decompression. Proton losses are only relevant at the smallest energies.

###  4.2 Cosmic ray energetics in superbubbles

Figure 6 displays the CR energetics in a fiducial SB described by the
parameters given in Table 1. About 10% of the mechanical power is transferred
into CRs. For small clusters ($N_{*}<100$), the CR production becomes
intermittent, in particular if the turbulence is not efficiently generated, as
shown by the yellow curve where we set $\eta_{T}=1\%$. In this case the
particles are not confined and the SB is similar to a collection of individual
SNR.

As shown by the zoomed-in window in the left panel of Figure 6, CR
acceleration proceeds in two steps. First, DSA at SNR injects nonthermal
particles inside the SB. In the following few 10 kyr, the high energy
particles rapidly escape. Meanwhile, low energy particles are reaccelerated in
the magnetic turbulence. After around 100 kyr, this second order
reacceleration compensates the escape and the total energy increases
progressively. However, the process of stochastic reacceleration is only
effective at low energies, such that when a substantial fraction of the
particles has been reaccelerated, the escape dominates once again and the
total energy decreases, until the next star explodes.

Figure 6: Time evolution of the total CR energy (left) and CR energy density
(right) inside the SB for clusters containing initially 100, 500 and 1000
stars. All the results are obtained with $\eta_{T}=30\%$ except those shown by
the yellow curves for which $\eta_{T}=1\%$.

The right panel of Figure 6 displays the time evolution of the CR energy
density inside the SB. Densities of the order of 1 to 10 eV/cm3 are expected
in rather small clusters ($N_{*}<100$), while larger clusters gathering
hundreds of stars may contain a CR energy density higher than 100 eV/cm3.

Our computation also shows that even when accounting for nonlinear effects,
CRs are mainly injected by SNRs (e.g. $95$% of the particles come from SNRs
for $\eta_{w}=0.1$, or $78$% for $\eta_{w}=0.5$). The fraction is slightly
smaller for more massive clusters (e.g. 64% for a cluster of 500 stars with
$\eta_{w}=0.5$) because the nonlinear effects are stronger. Yet, the
supernovae are dominant in most cases. This result matches the requirement
needed to reproduce the composition of galactic CRs (Tatischeff et al., 2021).

In fact, the energy transfer between the winds and the CRs proceeds mainly
indirectly through the stochastic reacceleration in turbulence. As the low-
and intermediate-energy particles are confined in the SB, this process is very
efficient. Nearly all the energy injected from the stars to the turbulence is
transferred to CRs. This could explain the high gamma-ray fluxes detected in
e.g. G25.18+0.26 (Katsuta et al., 2017), without requiring an unrealistically
small diffusion coefficient. This statement has to be confirmed by further
studies including the inhomogeneity of the turbulence and distribution
function. Indeed, if the stars are gathered at the centre of the cluster, they
are expected to excite Alfvén waves only in a small region of space. On the
other hand, the CR density is expected to be higher close to the centre of the
bubble, such that most of the particles could still be efficiently
reaccelerated.

Because the main source of matter inside the SB is the evaporation of
interstellar material at the shell interface, most of the material accelerated
in the SB by SNRs is expected to be of galactic composition. A discrepancy is
however expected because of the efficient injection of WR-enriched material in
the early phase of the cluster life. Although the fraction of CRs accelerated
at WTS is only about 5-10% of that accelerated at SNRs, this may be enough to
account for the observed overabundance of 22Ne over 20Ne in the galactic CRs,
in particular if very massive stars collapse without exploding (Kalyashova et
al., 2019; Gupta et al., 2020; Tatischeff et al., 2021).

###  4.3 Intermittency

Figure 7: Time evolution of the particle distribution function. The parameters
of the SB are given in Table 1. Top panel: 100 massive stars. Second panel:
500 massive stars. Third panel: 1000 massive stars. Bottom panel: 100 massive
stars, with $\eta_{T}=1\%$ instead of 30%. Colours correspond to different
momenta (see the colour scale in units of GeV/c at the top). Dark blue: 10
MeV/c. Blue: 0.13 GeV/c. Light blue: 1.6 GeV/c. Turquoise blue: 20 GeV/c.
Green: 240 GeV/c. Yellow: 2.9 TeV/c. Orange: 35 TeV/c. Red: 0.43 PeV/c. Dark
red: 5.3 PeV/c.

The time evolution of the particle distribution function is shown in Figure 7
for various initial numbers of massive stars. The low energy bands ($p\lesssim
0.1$ GeV/c) are always intermittent because low energy particles are
stochastically reaccelerated in the turbulence as soon as they are injected
into the SB. Intermediate energies are less intermittent because the
stochastic acceleration is less efficient and the escape is slow. In
particular, for standard clusters containing hundreds of massive stars, the
GeV band is nearly stationary, except if the turbulence level is low (bottom
panel of Figure 7), in which case even these rather low energies are not
confined and escape between supernova explosions. High energies are very
intermittent because the escape is faster than the interval between two
supernova explosions. Three regimes of CR production should therefore be
distinguished. Whenever the stochastic acceleration time is lower than the
average time interval between two supernovae, i.e. $p^{2}/2D_{p}(p)<\Delta
t_{\rm SN}$, the CR production is expected to be intermittent. This regime
extends up to the momentum $p_{1}$ defined by $p_{1}^{2}/2D_{p}(p_{1})=\Delta
t_{\rm SN}$. In the test-particle approximation, we obtain:

$\beta(p_{1})^{3}p_{1}=4{~{}\rm
MeV/c}\left(\frac{N_{*}}{100}\right)^{-3}\left(\frac{\eta_{T}L_{*}}{10^{51}{~{}\rm
erg/Myr}}\right)^{2}\left(\frac{\mathcal{T}_{\rm SB}}{35{~{}\rm
Myr}}\right)^{3}\\\ \times\left(\frac{n}{0.01~{}{\rm
cm}^{-3}}\right)^{2}\left(\frac{B}{10{~{}\text{\textmu
G}}}\right)\left(\frac{R_{b}}{50~{}{\rm pc}}\right)^{-6}\,.$ (33)

Because the damping of the turbulence by non-thermal particles is generally
moderate, the test-particle approximation is expected to give reliable
estimates on average, even though the nonlinear timescale of diffusion may
decrease considerably right after a supernova explosion.

Figure 8: Instantaneous SB spectra. Each curve corresponds to a different time
represented by the colour scale in units of Myr at the top. The SB parameters
are given in Table 1, with varying number of stars. Top left panel:
$N_{*}=100$. Top right panel: $N_{*}=100$ with $\eta_{T}=1\%$. Bottom left
panel: $N_{*}=500$. Bottom right panel: $N_{*}=1000$.

Then, if both the stochastic acceleration time as well as the escape time are
larger than the average time interval between two supernovae, i.e.
$p^{2}/2D_{p}(p),R_{s}^{2}/6D_{x}(p)<\Delta t_{\rm SN}$, the CR production is
expected to be nearly stationary. This regime extends up to the momentum
$p_{2}$ defined by $R_{s}^{2}/6D_{x}(p_{2})=\Delta t_{\rm SN}$, which
translates into:

$p_{2}=0.09{~{}\rm
GeV/c}\left(\frac{N_{*}}{100}\right)^{3}\left(\frac{\eta_{T}L_{*}}{10^{51}{~{}\rm
erg/Myr}}\right)^{2}\left(\frac{\mathcal{T}_{\rm SB}}{35{~{}\rm
Myr}}\right)^{-3}\\\ \times\left(\frac{n}{0.01~{}{\rm
cm}^{-3}}\right)\left(\frac{B}{10{~{}\text{\textmu G}}}\right)^{-5}\,.$ (34)

Finally, for $p>p_{2}$, the escape dominates and the cosmic ray production is
expected to be intermittent again. The three regimes can be identified in
Figure 5, comparing the Fermi II acceleration time (red curve), the escape
time (blue curve) and the time interval between supernovae (black line). If
the black line is below the other timescales, supernovae can sustain the
reacceleration and escape. In other cases, intermittent reacceleration or
escape is expected.

It may happen that $p_{2}<p_{1}$. In this case, there are not enough stars to
sustain a stationary CR production against the escape, even at intermediate
energies. This translates into an intermittent CR production at all energies,
as it is the case for the cluster of 100 stars with inefficient turbulence
generation ($\eta_{T}=1\%$, see bottom panel of Figure 7).

For our fiducial clusters containing respectively 100, 500 and 1000 massive
stars, Equations 33 and 34 provide respectively $p_{1}\approx 0.1$, 0.04, 0.01
GeV and $p_{2}\approx 1$ GeV, 3 TeV, 1 PeV, which is consistent with the
curves shown in Figure 7. In these cases, the test-particle estimates are
proved to be reliable indicators within one order of magnitude.

Eventually, one notices in the top and bottom panels of Figure 7 that the
winds provide a stationary background below the intermittent emission due to
supernova explosions. Although the wind contribution is subdominant below
about 10 GeV; it is non negligible at higher energies, where the escape is
efficient.

###  4.4 Spectra

Examples of spectra resulting from our nonlinear computation are shown in
Figure 8. The overall shape agrees qualitatively with the conclusions drawn by
Ferrand & Marcowith (2010), that is, a hard component at low momenta competes
with a steep component at high momenta. We denote $p_{*}$ the momentum of the
transition. As we expressed the turbulence spectrum as function of the
mechanical energy of the stars rather than as function of the background
magnetic field, our expression for $p_{*}$ differs from that given in Ferrand
& Marcowith (2010) even in the test-particle regime:

$\beta(p_{*})^{3}p_{*}\approx 0.1m_{p}c\left(\frac{B}{10~{}\text{\textmu
G}}\right)^{-2}\left(\frac{\eta_{T}L_{*}}{10^{51}~{}\text{erg/Myr}}\right)^{2}\\\
\times\left(\frac{n}{0.01~{}\text{cm}^{-3}}\right)^{-1/2}\left(\frac{R_{b}}{100~{}\text{pc}}\right)^{-3}\,.$
(35)

One notices that $p_{*}$ decreases as the SB expands, which is a consequence
of the dilution of the mechanical energy of the stars in an increasing volume.
In the fiducial cluster of 100 massive stars considered in this work, $p_{*}$
is found around a few to tens of GeV. In typical clusters of several hundreds
of stars, $p_{*}$ can hardly be higher than a few tens of GeV, even
considering unrealistic turbulence generation efficiencies, because in this
case nonlinearities are expected to regulate the stochastic acceleration
process by damping the waves.

The shapes of the SB spectra are overall well understood by the timescales of
the competing reacceleration and escape processes. One may wonder why the high
energy “universal” asymptotic solution of successive shock reacceleration
computed in Vieu et al. (2021) is never retrieved, in particular for the
cluster of 1000 massive stars, where the average time interval between
supernovae is about 30 kyr, that is, much smaller than the escape time at GeV
energies. This is because we implicitly assumed, by considering the average CR
spectrum $n(p)$ instead of the distribution function $f(x,p)$, that the
(re)accelerated CRs instantaneously homogenise inside the SB. Under these
circumstances, preaccelerated particles have small chances to be reaccelerated
by the next SNR shock, which only spans a few percent of the SB volume. Such a
situation corresponds to a loose cluster, with large distances between the
stars. The model will be refined in Section 5.2 in order to include the case
of a compact cluster.

Finally, although the contribution from the winds is globally subdominant, it
yet provides a non-negligible nearly stationary steep component (see Figure
8).

When integrated over the lifetime of the cluster, the spectra display an
enlarged “Fermi II bump” and then a steep high energy tail above 10-100 GeV,
as shown in Figure 9, where the spectra have been corrected as $\Phi(E)=4\pi
p^{2}n(p)/(v(p)\tau(p))$. This corresponds to the CR flux escaping from the
cluster in the ISM, except above 1 PeV where we disregarded the very high
energy flux escaping upstream of SNRs. Interestingly, the high energy flux is
a power law scaling as $\Phi(E)\propto E^{-2.2}$, which would reproduce the
$E^{-2.7}$ scaling measured around the Earth if one assumes that the CR
diffusion coefficient scales as $D(E)\propto E^{0.5}$ in the ISM, as expected
in the Iroshnikov-Kraichnan regime of turbulence. Eventually, the transition
between the maximum energy achievable in WTS (assumed to be about 100 TeV) and
that achieved at SNRs (about 1 PeV), produces a spectral break between 0.1 and
1 PeV which has a certain resemblance to the knee observed in the galactic CR
spectrum.

Figure 9: Average CR flux escaping the SB.

##  5 Two-zone model

In this section we show that it is possible to refine the modelling in order
to compute supernova reacceleration in compact clusters in a more realistic
way as well as to consider the effect of the SB supershell.

###  5.1 Diffusion in a two-zone model

Let us split the interior of the spherical SB into two zones, as schemed in
Figure 10. The top panel presents the general setup while the two other panels
display two limit cases of physical interest, which will be discussed below.
In this section we stick to a general formalism, assuming that the inner
region has a radius $R_{1}$ and the outer region has a radius $R_{b}$, which
is the SB radius. The boundaries of both areas are therefore separated by a
length equal to $R_{b}-R_{1}$. The two regions are characterised by diffusion
coefficients $D_{1}$ and $D_{2}$, respectively. We denote $f_{1}$ and $f_{2}$
the average distributions of particles in both regions such that the
corresponding energy densities are $n_{1}=V_{1}f_{1}$ and $n_{2}=V_{2}f_{2}$,
where $V_{1}=4/3\pi R_{1}^{3}$ and $V_{2}=V_{\rm SB}-V_{1}$ with $V_{\rm SB}$
the SB volume.

Figure 10: Top: Sketch of the two-zone model. Bottom left: Limit case of a
compact cluster, where zone 1 is the region of supernova reacceleration.
Bottom right: Limit case of a supershell. Figure 11: Average spectra around a
compact cluster modelled as a two-zone model (inner/outer). Each curve
corresponds to a different time represented by the colour scale in units of
Myr at the top. The number of massive stars is assumed to be 100 (left), 500
(middle) and 1000 (right). The top panels show the spectrum in the inner
region spanned by SNRs while the bottom panels show the spectrum in the outer
region, which fills most of the SB volume.

The CR transport in the SB is then described as a two-zone diffusion-losses
model. The transport equation reads:

$\partial_{t}f_{i}=\nabla\cdot\left(D_{i}\cdot\nabla
f_{i}\right)+\mathcal{D}_{i}[f_{i}]+\mathcal{L}_{i}[f_{i}]+q_{i}\,,$ (36)

where $i=1$ (resp. $i=2$) in the inner (resp. outer) region,
$\mathcal{D}_{i}[f_{i}]=(1/p^{2})\partial_{p}\left(p^{2}D_{p,i}\partial_{p}f_{i}\right)$
is the stochastic reacceleration operator, with $D_{p,i}$ the momentum
diffusion coefficient and $\mathcal{L}_{i}$ an operator encoding the various
loss processes in both zones, and $q_{i}$ encodes the background wind
injection in zone 1 ($q_{1}=q_{w}$, with $q_{w}$ defined in Eq. 13) and is
zero in zone 2. Instead of solving directly Eq. 36, we average the
distribution functions $f_{i}$. Denoting $\tau_{12}$ the typical time it takes
for particles to diffuse from zone 1 to zone 2 (and conversely), and
$\tau_{2}$ the escape time from zone 2 to the ISM, we obtain the following
system of differential equations:

$\left\\{\begin{array}[]{lll}\partial_{t}n_{1}&=&\frac{n_{2}-n_{1}}{\tau_{12}}+\mathcal{D}_{1}[n_{1}]+\mathcal{L}_{1}[n_{1}]+q_{w}\,,\\\
\partial_{t}n_{2}&=&\frac{n_{1}-n_{2}}{\tau_{12}}-\frac{n_{2}}{\tau_{2}}+\mathcal{D}_{2}[n_{2}]+\mathcal{L}_{2}[n_{2}]+\phi_{\rm
ISM}\,,\\\ \end{array}\right.$ (37)

where $\phi_{\rm ISM}$ is the flux of interstellar particles entering the
outer region.

To estimate the typical residence time of particles in between zone 1 and zone
2, we compute the diffusion time between the centre of the bubble ($r=0$) and
the middle of the shell ($r=(R_{1}+R_{b})/2$). To do so, we follow the
approach presented in e.g. Drury (1983). The solution of the stationary
diffusion equation in spherical symmetry,
$\partial_{r}(r^{2}D\partial_{r}f)=r^{2}Q\delta(r)$, with a constant injection
of particles $Q$ at $r=0$, reads:

$\displaystyle\left\\{\begin{array}[]{l}r<R_{1}~{}~{}\text{(zone 1)}\implies
f=\frac{Q}{4\pi D_{1}}\left(\frac{1}{R_{1}}-\frac{1}{r}\right)+\frac{Q}{4\pi
D_{2}}\left(\frac{2}{R_{1}+R_{b}}-\frac{1}{R_{1}}\right)\,,\\\
r>R_{1}~{}~{}\text{(zone 2)}\implies f=\frac{Q}{4\pi
D_{2}}\left(\frac{2}{R_{1}+R_{b}}-\frac{1}{r}\right)\,.\\\ \end{array}\right.$
(40)

Integrating this solution from $0$ to $R_{1}$ and then from $R_{1}$ to
$(R_{1}+R_{b})/2$, we obtain $N$, the total number of particles in between
both regions:

$N=\frac{QR_{1}^{2}}{6D_{1}}\left(1+\frac{\rho_{b}}{\delta}+\frac{\rho_{b}^{2}}{4\delta}\right)\,,$
(41)

where we set $\rho_{b}\equiv(R_{b}-R_{1})/R_{1}$ and $\delta\equiv
D_{2}/D_{1}$. The typical residence time is therefore identified as:

$\tau_{12}=\frac{R_{1}^{2}}{6D_{1}}\left(1+\frac{\rho_{b}}{\delta}+\frac{\rho_{b}^{2}}{4\delta}\right)\,,$
(42)

which can also be interpreted as the time it takes for a particle to be
transferred from one zone to the other.

Eventually, we compute the escape time from zone 2 to the interstellar medium
by integrating Eq. 36 over the volume of zone 2, assuming that the spatial
diffusion coefficient in the ISM is much larger than that in zone 2. This
provides, by virtue of the divergence theorem:

$\partial_{t}n_{2}=\oiint\left(D_{2}\cdot\nabla
f\right)\cdot\mathrm{d}\boldsymbol{S}+\mathcal{D}_{2}\left[n_{2}\right]+\mathcal{L}_{2}\left[n_{2}\right]\,.$
(43)

Then we estimate the gradient of particles across the surface $S_{b}=4\pi
R_{b}^{2}$ by means of a linear extrapolation as $\nabla f_{2}\sim
2\alpha(f_{\rm ISM}-n_{2}/V_{2})/(R_{b}-R_{1})$, with $\alpha$ a geometrical
factor. The surface integral becomes trivial and we eventually get:

$\partial_{t}n_{2}=\frac{6\alpha D_{2}R_{b}^{2}(V_{2}f_{\rm
ISM}-n_{2})}{(R_{b}^{2}+R_{b}R_{1}+R_{1}^{2})(R_{b}-R_{1})^{2}}+\mathcal{D}_{2}\left[n_{2}\right]+\mathcal{L}_{2}\left[n_{2}\right]\\\
-\iint_{S_{1}}\left(D_{2}\cdot\nabla f\right)\cdot\mathrm{d}\boldsymbol{S}\,,$
(44)

where $f_{\rm ISM}$ is the diffuse CR distribution function in the ISM and
$S_{1}$ the surface of the interface between the two zones. The surface
integral over the inner surface $S_{1}$ encodes the transfer of particles from
zone 1 to zone 2 which has been already computed. From Eq. 44 we identify the
escape time from zone 2 to the ISM as well as the flux of interstellar
particles entering the SB:

$\displaystyle\tau_{2}$
$\displaystyle=\frac{(R_{b}^{2}+R_{b}R_{1}+R_{1}^{2})(R_{b}-R_{1})^{2}}{6\alpha
D_{2}R_{b}^{2}}\,,$ (45) $\displaystyle\phi_{\rm ISM}$
$\displaystyle=\frac{4\pi}{3\tau_{2}}\left(R_{b}^{3}-R_{1}^{3}\right)f_{\rm
ISM}\,.$ (46)

In the limit case $R_{1}\ll R_{b}$, one should set $\alpha=1$ to get
$\tau_{2}=R_{b}^{2}/(6D_{2})$ which is the diffusion time in spherical
symmetry. On the other hand, if $R_{1}\approx R_{b}$, one should set
$\alpha=4$ in order to get $\tau_{2}=((R_{b}-R_{1})/2)^{2}/(2D_{2})$, which is
the diffusion time in one dimension.

Now that all timescales have been determined, Eq. 37 can be written in matrix
form and solved iteratively. Such two-zone model can be used to compute for
instance the acceleration and transport of CRs in a SB where turbulence is
only efficiently generated around the stellar cluster, or to account for the
residence time of particles in the inner region spanned by successive SNRs, or
to model the shielding of the interstellar CRs by the dense shell, or last but
no least, to probe the modulation of the CR spectra induced by the supershell
with associate non-thermal emission of photons.

###  5.2 Cosmic ray reacceleration in compact clusters

As mentioned in the end of Section 4, the one-zone model does not properly
account for CR reacceleration in compact clusters, for it is implicitly
assumed that the accelerated particles homogenise instantaneously in the SB
after the passage of a remnant. The two-zone model can be used to describe the
reacceleration process more accurately. Indeed, the inner region can be
identified with the volume spanned by SNRs $\mathcal{V}_{\rm SNR}$444Because
this volume depends on the density of the SB interior, it is time-dependent
and the radius of the inner region is defined dynamically.. The SNR radius is
typically 10% the SB radius such that this situation is close to the limit
case $R_{1}\ll R_{b}$ sketched in the bottom left panel of Figure 10. For
simplicity we assume that the turbulence is homogeneous in both zones,
although in reality it should probably be stronger in the inner zone close to
the energy deposition region. In the case of a compact cluster, a collective
wind termination shock (WTS) is expected to form at a distance of a few tens
of pc around the cluster. Morlino et al. (2021) showed that a standard
$p^{-4}$ spectrum is expected up to a maximum energy of hundreds of TeV, at
most PeV. Therefore, Eq. 13 satisfactorily reproduces the injection of
particles at a collective WTS. In this section, we stick to a maximum energy
of 100 TeV achieved at the WTS, as a conservative estimate. Eventually, we
neglect the flux of interstellar particles leaking inside the SB ($\phi_{\rm
ISM}=0$). This assumption is justified for clusters hosting several hundreds
of massive stars as the energy density of the confined particles is expected
to be much larger than the interstellar CR energy density. Besides, low energy
ISM particles will be depleted by the ionisation losses in the SB shell, while
the steep high energy tail will always be negligible. For small clusters, the
ISM may only provide a noticeable contribution to the SB content around GeV
energies. Although it is straightforward to include this contribution in our
model, we chose to discard it in order to make the interpretation of the
spectra easier. In particular, this allows to directly compare the output of
the two-zone model with the results of Section 4.

Figure 11 displays the spectra obtained in this refined model. The first thing
to notice is that the overall normalisation is about one to two orders of
magnitude below that obtained using a one-zone model (Figure 8). The energy
transfer between the SNR and the CRs is much less efficient because of the
backreaction of the CRs onto the shocks in the inner region. Indeed,
preaccelerated CRs do not have time to entirely dilute inside the SB between
two supernova explosions. The particle energy density in the inner region is
high, which slows down the shock precursor and makes the acceleration of low-
energy particles less efficient (see Vieu et al., 2021). On the other hand, if
the massive stars are numerous enough, the time interval between supernovae
can be smaller than the escape time at intermediate or even high energies and
the particles are efficiently reaccelerated without leaving the region. As
expected from the “benchmark” asymptotic solution discussed in Vieu et al.
(2021), the spectrum resulting from the reacceleration by successive shocks is
concave, with a high energy slope close to 3.5 in the case where the escape of
the particles is completely neglected. As shown in the middle panel of Figure
11, concave spectra are indeed intermittently retrieved in a cluster
containing initially 500 massive stars. For a cluster of 1000 stars, as shown
in the right panel of Figure 11, the reacceleration is very efficient and the
benchmark asymptotic spectrum may be achieved up to the highest energies. Then
the concave spectrum realised in the inner region is transferred in the outer
region, where low energy bands experience a stochastic reacceleration. The
particles will eventually interact with the shell and concave gamma-ray
spectra should be expected, which could be a typical signature of massive
compact clusters.

Interestingly, hard gamma-ray spectra with slight concavities were recently
observed in the G25 region (Katsuta et al., 2017). Although the nature of the
source is yet unclear, it was suspected that a star forming region could be
responsible for this extended emission. Such spectra are indeed possible
signatures of shock reacceleration in a SB surrounding the hypothetical
G25.18+0.26 young stellar association.

Figure 12: Average spectra in a two-zone (top panels: interiors, bottom
panels: shells) model with two different values of the shell confinement
coefficient $\rho_{b}/\delta$ = 10 (left panels), 1000 (right panels). Each
curve corresponds to a different time represented by the colour scale in units
of Myr at the top.

###  5.3 Modelling the supershell

The two-zone model can also be used to account for the effect of the dense
supershell in which the interstellar matter and magnetic field swept-up by the
SB forward shock accumulate. The interstellar magnetic field is compressed and
strengthened in the direction tangential to the shell. In order to escape in
the interstellar medium, particles must diffuse mainly perpendicularly to this
strong background field. Deriving the expression for a perpendicular diffusion
coefficient from first principles is far from being trivial. On the other
hand, heuristic arguments as well as numerical simulations suggest that
$D_{\perp}\approx(\delta B/B)^{\alpha_{B}}D_{\parallel}$, with
$\alpha_{B}\approx 4$ (e.g. Casse et al., 2001; Mertsch, 2020). Therefore, the
diffusion coefficient in the shell could well be orders of magnitude below
that in the superbubble interior, which would enhance the confinement of the
particles. This can be accounted for in the two-zone model as the parameter
$\delta$ in Eq. 42 is precisely the ratio of the diffusion coefficients.
Eventually we need to specify the shell density and thickness. For simplicity
we assume that the shell thickness is 10% of the SB radius at any time, that
is, $\rho_{b}=0.1$. This parameter is not well constrained. Observations and
numerical simulations provide estimates ranging from a few percent (e.g.
Joubaud et al., 2019) to tens of percent (e.g. Krause et al., 2013) and even
up to 50% (e.g. Yadav et al., 2017). The hydrogen number density in the shell
$n_{H}$ follows by assuming that all the mass swept-up by the SB forward shock
has accumulated in the shell, which provides $n_{H}=n_{\rm
ISM}(1+\rho_{b})^{3}/((1+\rho_{b})^{3}-1)$. In the following we assume $n_{\rm
ISM}=10$ cm-3.

Eq. 42 predicts that the confinement of the particles will be enhanced when
$\rho_{b}/\delta>1$. Indeed the left panels of Figure 12 show that for
$\delta=0.01$ ($\rho_{b}/\delta=10$), the emission in a SB hosting 100 massive
stars is much less intermittent, which results in overall smoother and steeper
spectra because the flat wind component is always subdominant. For a high
“shell confinement coefficient” $\rho_{b}/\delta=10^{3}$ (right panels of
Figure 12), the particles are strongly confined. In this case, the CR energy
density in the shell is so high that the turbulence cascade terminates close
to the injection scale: the shell is not turbulent. This demonstrates that low
values of the parameter $\delta$ (that is, strong confinements) are not
unrealistic, although a self-consistent model computing $\delta$ from first
principles and accurately accounting for the backreaction of the particles on
the perpendicular transport is still missing.

Whenever the shell efficiently enhances the confinement, the particles
reaccelerated at SNRs stay inside the SB in such a way that the interior is
characterised by a homogeneous concave spectrum which converges toward the
asymptotic solution discussed in Vieu et al. (2021). Eventually, one notices
in the bottom panels of Figure 12 that the spectra are cut off below 1 GeV,
which is due to the enhanced energy losses suffered by low-energy particles
within the dense shell.

##  6 Superbubble contribution to galactic cosmic rays

Let us now provide a rough estimate of the contribution of SBs enclosed by
supershells to galactic CRs. The differential flux of particles reaching the
Earth from SBs is estimated as:

$\frac{\mathrm{d}N}{\mathrm{d}E}=\frac{f_{\rm SB}}{\langle V_{\rm
SB}\rangle}\frac{T_{\rm disk}(p)c}{\tau_{2}}p^{2}n_{\rm shell}(p)\,,$ (47)

where $f_{\rm SB}\approx 20\%$ is the volume filling factor of the hot gas in
the solar vicinity (Ferrière, 1998), $\langle V_{\rm SB}\rangle\approx 0.004$
kpc3 is the volume of a typical superbubble, $T_{\rm disk}(p)\approx 3(p/{\rm
GeV})^{-\delta_{\rm gal}}~{}{\rm Myr}$ is the residence time in the galactic
disk, with $\delta_{\rm gal}\approx 0.3-0.5$ (e.g. Génolini et al., 2019) and
$\tau_{2}$ is the diffusion time in the shell given by Eq. 45. Eq. 47 provides
an order of magnitude estimate of the flux normalisation without introducing
any arbitrary parameter. Figure 13 shows the differential flux computed from
the shell energy density $n_{\rm shell}$ obtained in the previous subsection
for $N_{*}=100$ and a confinement parameter $\rho_{b}/\delta=10$ (bottom left
panel of Figure 12).

Figure 13: Superbubble contribution to the galactic cosmic ray proton
spectrum, assuming a galactic diffusion coefficient with index
$\delta_{gal}=0.5$ (blue), $\delta_{gal}=0.4$ (red), $\delta_{gal}=0.3$
(green).

As seen in Figure 13, the flux normalisation is close to that of the available
data. The main parameter driving the CR density in SBs is the injection
efficiency at SNRs, that is the fraction of thermal particles injected in the
accelerator. Increasing the injection efficiency is however not expected to
change the result dramatically, for shock acceleration becomes less efficient
in the nonlinear regime. The agreement between this rough prediction and the
data is overall acceptable and this preliminary estimate of the contribution
of SBs to the galactic CRs show that they can account for the observed CR flux
with a realistic residence time and thus a realistic galactic diffusion
coefficient.

Let us eventually recall that the SB contribution to the galactic CR spectrum
computed beyond 1 PeV is not reliable, for we neglected the flux of particles
escaping upstream of SNRs beyond the maximum energy. An accurate description
of this component requires to solve the time-dependent reacceleration problem
in the nonlinear regime, at least phenomenologically, which could provide an
additional steep power law component at energies close to that of the knee.
Such promising analysis is left for future work.

##  7 Gamma-ray spectra

The gamma-ray emission of supershells due to neutral pion decay can eventually
be computed from the cosmic ray energy spectra as:

$\Phi_{\gamma}(E_{\gamma})=\frac{n_{H}c}{4\pi D_{S}^{2}}\times\\\
\int_{T_{p,min}}\mathrm{d}T_{p}4\pi p^{2}n_{\rm
shell}(p(T_{p}))\epsilon_{nh}(T_{p})\frac{\mathrm{d}\sigma_{pp}}{\mathrm{d}E_{\gamma}}(T_{p},E_{\gamma})\,,$
(48)

where $n_{H}$ is the hydrogen number density of the shell, $D_{S}$ is the
distance of the source, $T_{p,min}$ the threshold energy for neutral pion
production in p-p interactions, $\mathrm{d}\sigma/\mathrm{d}E_{\gamma}$ is the
differential cross-section of gamma-ray production and $\epsilon_{nh}$ is the
nuclear enhancement factor accounting for nucleus-nucleus interactions. The
computation follows the prescription of Kafexhiu et al. (2014).

Figure 14: Gamma-ray spectral energy distributions from p-p interactions in a
standard superbubble shell of hydrogen number density $n_{H}=40$ cm-3 ($n_{\rm
ISM}=10$ cm-3), with a shell confinement coefficient $\rho_{b}/\delta=10$
(left panel), and 1000 (right panel). The source is a loose cluster of 100
massive stars located at 1.5 kpc from the Earth. Spectra are plotted at
various times from 1 to 20 Myr corresponding to the various colours from blue
to red (see the colour scale at the top). The data points show the gamma-ray
flux measured in the Cygnus region: Fermi-LAT data in red (Abdollahi et al.,
2020) and HAWC data in blue (Abeysekara et al., 2021).

Figure 14 shows the resulting gamma-ray spectral energy distributions for a
cluster of 100 massive stars and two values of the parameter $\delta$. It is
striking to see how the modification of a single parameter can produce many
different spectral shapes. In standard SBs surrounding loose clusters, the CRs
are not very efficiently reaccelerated (left panel). The fluxes escaping
through the shell produce rather steep power law gamma-ray spectra, with a low
normalisation and characterised by a strong intermittency. The flux is
generally below the detection threshold of the Fermi-LAT telescope, which
could explain why several massive star clusters are not associated with gamma-
ray emissions (Maurin et al., 2016), and in particular the well-known Orion-
Eridani superbubble (Joubaud et al., 2020). On the other hand, with an
enhanced shell confinement coefficient $\rho_{b}/\delta$, the CRs are trapped
in the shell, which results in a more efficient production of gamma-rays
(right panel). The spectra at early times, driven by stellar winds, are
slightly steeper than $E^{-2}$ (e.g. $\Phi_{\gamma}\propto
E^{-2.1}-E^{-2.2}$), which is consistent with what is observed in the Cygnus-X
star forming region. In particular, the spectrum at 3 Myr shows a very good
agreement with the normalisation and shape of the spectral energy distribution
measured in the Cygnus region, although we did not adjust any parameter to
obtain the best fit. A thorough comparison with Fermi and HAWC data is left
for future work. Eventually, hard power law or concave spectra are expected if
the CRs are efficiently confined and reaccelerated, as it starts to be the
case for $\rho_{b}/\delta=1000$ (right panel of Figure 14) after about $5-10$
Myr. This is expected if the perpendicular diffusion is suppressed in the
shell or if the particles are efficiently reaccelerated in compact clusters.

Several observations of extended hard gamma-ray spectra possibly associated
with young massive stellar clusters or SB cavities have been reported in the
last decade, even though it is not established that these gamma-ray emissions
are produced by collective plasma effects, in particular because the
contribution of point-like sources, e.g. pulsar wind nebulae, is difficult to
isolate (Tibaldo et al., 2021). Slightly concave gamma-ray spectral energy
distributions were observed in the G25 region (Katsuta et al., 2017), flat
gamma-ray spectra were observed in the very compact and very massive
Westerlund 2 young stellar cluster (Aharonian et al., 2007; Yang et al.,
2018), and flat gamma-ray spectra extending up to TeV energies were observed
in the Westerlund 1 region, in the vicinity of the most massive star cluster
of our Galaxy (Abramowski et al., 2012; Ohm et al., 2013). It seems that all
these observations, despite differing quantitatively or even qualitatively,
could be explained by SB proton spectra. Our SB modelling would allow to
constrain the properties of the supershells, turbulence, and massive star
clusters in order to reproduce gamma-ray data. In particular, the model being
self-consistent, it predicts not only spectral shapes but realistic
normalisations as well. Eventually, protons of energies above 1 PeV could
produce ultra high energy photons above the sensitivity of the LHAASO
observatory, which recently detected a photon of 1.4 PeV in the direction of
the Cygnus OB2 massive cluster (Cao et al., 2021). We however recall that the
in order to accurately describe the CR spectrum produced in SBs near the
maximum energy, a time-dependent nonlinear model of CR reacceleration at SNR
shocks is needed.

##  8 Conclusions

Using a self-consistent model accounting for all relevant ingredients
rederived from first principles, we computed the acceleration of protons in
superbubble environments. We extended the model of Ferrand & Marcowith (2010)
in order to account for the dynamical evolution of the environment as well as
the stellar winds, the losses and the effect of the shell. We also refined the
model to include the backreaction of the particles onto the turbulence and
onto the shocks, the latter being properly computed using an up-to-date semi-
analytical model of nonlinear diffusive shock reacceleration. Our model
complements the work of Bykov (2001) who focused on the early phase of
stochastic acceleration in strong supersonic turbulence.

We found that the stars efficiently transfer their energy into non-thermal
particles by means of shock acceleration as well as stochastic acceleration in
turbulence, especially around compact clusters where particles are efficiently
reaccelerated in the inner region and in the case where a dense magnetised
supershell prevents the particles to escape in the ISM. High CR energy
densities are generally achieved in the SB, which calls for nonlinear models,
and in some cases the turbulence can be completely suppressed by the non-
thermal particles confined in the shell.

When the confinement is less efficient, e.g. in the limit of a thin supershell
or in small clusters, the typical spectra are rather intermittent, displaying
a typical “Fermi II bump” from the injection energy to the GeV band, and then
transitions toward a steep power law produced by the nearly stationary wind
contribution as well as the intermittent supernovae modulated by the escape.
Providing the level of turbulence is about a few percent, which requires the
stars to transfer a few tens of percent of their mechanical power into
hydromagnetic waves, the low energy particles are efficiently reaccelerated,
which gives rise to steep escape fluxes typically scaling as $E^{-2.2}$, which
is close to what is needed in order to account for the diffuse CR spectrum
observed near Earth. When the confinement of the particles is enhanced, the
spectra harden with typically flat gamma-ray signatures, or even small
concavities with somewhat hard components (about $E^{-1.8}$ at high energies)
due to the successive reaccelerations of the confined particles.

The variability of the spectra from cluster to cluster, and also during the
lifetime of a given cluster, provides a simple answer to the puzzling
discrepancies between the recent gamma-ray observations of SBs and star
clusters. Indeed, some SBs are not detected in gamma-rays (e.g. the Orion-
Eridani SB or the Rosette nebula), others display rather steep power law
spectra (e.g. the Cygnus region), or flat energy distributions (e.g. the
Westerlund 1 and 2 regions), and some data even suggest slight concavities
(e.g. in the G25 region). The specific gamma-ray signature of a SB can
therefore be used to constrain the properties of a given massive star cluster
(its number of massive stars, if it is compact or loose…), the surrounding
environment (the magnetic field, the turbulence level…), as well as the
properties of the supershell (the thickness, the density…). Although the model
depends on several parameters, some of them such as the cluster and shell
properties can hopefully be constrained by multiwavelength observations. The
main difficulty is probably to infer the properties of the turbulence
(intensity and spectrum), which at the moment are unknown. Using our model,
gamma-ray spectra could be used indirectly to constrain the diffusion
coefficients. Indeed, the momentum diffusion coefficient drives the CR
spectrum slope at low energies by means of the stochastic reacceleration,
while the spatial diffusion coefficient drives the CR spectrum slope at high
energies as well as the intensity of the gamma-ray flux, because it determines
the efficiency of the confinement. Such analysis could be applied for example
to the Cygnus region.

Eventually, the overall SB contribution to the galactic CR population can be
estimated using Monte-Carlo samplings. Thorough statistical computations
confronted to the observed CR spectrum could be used to probe the most likely
SB parameters, in particular the magnetic fields, turbulence levels and shell
properties, as well as their variance. The modulation of the galactic spectrum
due to the local bubble could also be included in such computation using our
two-zone modelling. Indeed, this model allows to compute the screening of the
interstellar CR flux when it penetrates a dense supershell, an idea which was
already raised in the 80’s (Streimatter et al. 1985, see also Bouyahiaoui et
al. 2019; Silsbee & Ivlev 2019; Phan 2020).

## Acknowledgements

TV acknowledges A. Marcowith, E. Parizot, V.H.M. Phan, L.M. Bourguinat for
helpful discussions and suggestions. SG and VT acknowledge support from Agence
Nationale de la Recherche (grant ANR-17-CE31-0014).

## Data Availability

No new data were generated or analysed in support of this research.

## References

  * Abdollahi et al. (2020) Abdollahi S., et al., 2020, ApJS, 247, 33
  * Abeysekara et al. (2021) Abeysekara A. U., et al., 2021, Nature Astronomy,
  * Abramowski et al. (2012) Abramowski A., et al., 2012, A&A, 537, A114
  * Abramowski et al. (2015) Abramowski A., et al., 2015, Science, 347, 406
  * Ackermann et al. (2011) Ackermann M., et al., 2011, Science, 334, 1103
  * Aharonian (2004) Aharonian F. A., 2004, Very high energy cosmic gamma radiation : a crucial window on the extreme Universe, doi:10.1142/4657.
  * Aharonian et al. (2007) Aharonian F., et al., 2007, A&A, 467, 1075
  * Aharonian et al. (2019) Aharonian F., Yang R., de Oña Wilhelmi E., 2019, Nature Astronomy, 3, 561
  * Berezinskii et al. (1990) Berezinskii V. S., Bulanov S. V., Dogiel V. A., Ptuskin V. S., 1990, Astrophysics of cosmic rays
  * Blasi et al. (2005) Blasi P., Gabici S., Vannoni G., 2005, MNRAS, 361, 907–918
  * Bouyahiaoui et al. (2019) Bouyahiaoui M., Kachelriess M., Semikoz D. V., 2019, J. Cosmology Astropart. Phys., 2019, 046
  * Brunetti et al. (2004) Brunetti G., Blasi P., Cassano R., Gabici S., 2004, MNRAS, 350, 1174
  * Butt & Bykov (2008) Butt Y. M., Bykov A. M., 2008, ApJ, 677, L21
  * Bykov (2001) Bykov A. M., 2001, Space Sci. Rev., 99, 317
  * Bykov (2014) Bykov A. M., 2014, A&ARv, 22, 77
  * Bykov & Fleishman (1992a) Bykov A. M., Fleishman G. D., 1992a, Soviet Astronomy Letters, 18, 95
  * Bykov & Fleishman (1992b) Bykov A. M., Fleishman G. D., 1992b, MNRAS, 255, 269
  * Bykov & Toptygin (1990) Bykov A. M., Toptygin I. N., 1990, Zhurnal Eksperimentalnoi i Teoreticheskoi Fiziki, 98, 1255
  * Bykov & Toptygin (1993) Bykov A. M., Toptygin I., 1993, Physics Uspekhi, 36, 1020
  * Bykov et al. (1995) Bykov A. M., Ptuskin V. S., Toptygin I. N., 1995, in International Cosmic Ray Conference. p. 337
  * Bykov et al. (2013) Bykov A. M., Gladilin P. E., Osipov S. M., 2013, MNRAS, 429, 2755
  * Bykov et al. (2015) Bykov A. M., Ellison D. C., Gladilin P. E., Osipov S. M., 2015, MNRAS, 453, 113
  * Bykov et al. (2018) Bykov A. M., Ellison D. C., Gladilin P. E., Osipov S. M., 2018, Advances in Space Research, 62, 2764
  * Bykov et al. (2020) Bykov A. M., Marcowith A., Amato E., Kalyashova M. E., Kruijssen J. M. D., Waxman E., 2020, Space Sci. Rev., 216, 42
  * Cao et al. (2021) Cao Z., et al., 2021, Nature
  * Casse et al. (2001) Casse F., Lemoine M., Pelletier G., 2001, Phys. Rev. D, 65, 023002
  * Castor et al. (1975) Castor J., McCray R., Weaver R., 1975, ApJ, 200, L107
  * Cesarsky & Montmerle (1983) Cesarsky C. J., Montmerle T., 1983, Space Sci. Rev., 36, 173
  * Cooper et al. (2004) Cooper R. L., Guerrero M. A., Chu Y.-H., Chen C. H. R., Dunne B. C., 2004, ApJ, 605, 751
  * Drury (1983) Drury L. O., 1983, Reports on Progress in Physics, 46, 973
  * Dunne et al. (2001) Dunne B. C., Points S. D., Chu Y.-H., 2001, ApJS, 136, 119
  * Dunne et al. (2003) Dunne B. C., Chu Y.-H., Chen C. H. R., Lowry J. D., Townsley L., Gruendl R. A., Guerrero M. A., Rosado M., 2003, ApJ, 590, 306
  * Ebinger et al. (2020) Ebinger K., Curtis S., Ghosh S., Fröhlich C., Hempel M., Perego A., Liebendörfer M., Thielemann F.-K., 2020, ApJ, 888, 91
  * Eilek (1979) Eilek J. A., 1979, ApJ, 230, 373
  * Fang et al. (2019) Fang K., Bi X.-J., Yin P.-F., 2019, MNRAS, 488, 4074
  * Ferrand & Marcowith (2010) Ferrand G., Marcowith A., 2010, A&A, 510, A101
  * Ferrière (1998) Ferrière K., 1998, ApJ, 503, 700
  * Finke & Dermer (2012) Finke J. D., Dermer C. D., 2012, ApJ, 751, 65
  * Fryer (1999) Fryer C. L., 1999, ApJ, 522, 413
  * Fryer (2003) Fryer C. L., 2003, Classical and Quantum Gravity, 20, S73
  * 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
  * Gaggero et al. (2018) Gaggero D., Zandanel F., Cristofari P., Gabici S., 2018, MNRAS, 475, 5237
  * Gallegos-Garcia et al. (2020) Gallegos-Garcia M., Burkhart B., Rosen A. L., Naiman J. P., Ramirez-Ruiz E., 2020, ApJ, 899, L30
  * Génolini et al. (2019) Génolini Y., et al., 2019, Phys. Rev. D, 99, 123028
  * Gupta et al. (2020) Gupta S., Nath B. B., Sharma P., Eichler D., 2020, MNRAS, 493, 3159
  * Higdon & Lingenfelter (2005) Higdon J. C., Lingenfelter R. E., 2005, ApJ, 628, 738
  * Higdon & Lingenfelter (2006) Higdon J. C., Lingenfelter R. E., 2006, Advances in Space Research, 37, 1913
  * Higdon et al. (1998) Higdon J. C., Lingenfelter R. E., Ramaty R., 1998, ApJ, 509, L33
  * Jaskot et al. (2011) Jaskot A. E., Strickland D. K., Oey M. S., Chu Y. H., Garc\́mathrm{i}a-Segura G., 2011, ApJ, 729, 28
  * Joubaud et al. (2019) Joubaud T., Grenier I. A., Ballet J., Soler J. D., 2019, A&A, 631, A52
  * Joubaud et al. (2020) Joubaud T., Grenier I. A., Casandjian J. M., Tolksdorf T., Schlickeiser R., 2020, A&A, 635, A96
  * Kafexhiu et al. (2014) Kafexhiu E., Aharonian F., Taylor A. M., Vila G. S., 2014, Phys. Rev. D, 90, 123014
  * Kalyashova et al. (2019) Kalyashova M. E., Bykov A. M., Osipov S. M., Ellison D. C., Badmaev D. V., 2019, Journal of Physics: Conference Series, 1400, 022011
  * Katsuta et al. (2017) Katsuta J., Uchiyama Y., Funk S., 2017, ApJ, 839, 129
  * Klepach et al. (2000) Klepach E. G., Ptuskin V. S., Zirakashvili V. N., 2000, Astroparticle Physics, 13, 161
  * Krause et al. (2013) Krause M., Fierlinger K., Diehl R., Burkert A., Voss R., Ziegler U., 2013, A&A, 550, A49
  * Krumholz et al. (2019) Krumholz M. R., McKee C. F., Bland-Hawthorn J., 2019, ARA&A, 57, 227
  * Limongi & Chieffi (2006) Limongi M., Chieffi A., 2006, ApJ, 647, 483
  * Lingenfelter (2018) Lingenfelter R. E., 2018, Advances in Space Research, 62, 2750
  * Lopez et al. (2020) Lopez L. A., Grefenstette B. W., Auchettl K., Madsen K. K., Castro D., 2020, ApJ, 893, 144
  * Mac Low & McCray (1988) Mac Low M.-M., McCray R., 1988, ApJ, 324, 776
  * Maddox et al. (2009) Maddox L. A., Williams R. M., Dunne B. C., Chu Y. H., 2009, ApJ, 699, 911
  * Mannheim & Schlickeiser (1994) Mannheim K., Schlickeiser R., 1994, A&A, 286, 983
  * Massey et al. (1995) Massey P., Johnson K. E., Degioia-Eastwood K., 1995, ApJ, 454, 151
  * Maurin et al. (2016) Maurin G., Marcowith A., Komin N., Krayzel F., Lamanna G., 2016, A&A, 591, A71
  * McCray & Kafatos (1987) McCray R., Kafatos M., 1987, ApJ, 317, 190
  * Melrose & Pope (1993) Melrose D. B., Pope M. H., 1993, Proceedings of the Astronomical Society of Australia, 10, 222
  * Mertsch (2020) Mertsch P., 2020, Ap&SS, 365, 135
  * Miller & Roberts (1995) Miller J. A., Roberts D. A., 1995, ApJ, 452, 912
  * Montmerle (1979) Montmerle T., 1979, ApJ, 231, 95
  * Morlino et al. (2021) Morlino G., Blasi P., Peretti E., Cristofari P., 2021, MNRAS,
  * Norman & Ferrara (1996) Norman C. A., Ferrara A., 1996, ApJ, 467, 280
  * Ohm et al. (2013) Ohm S., Hinton J. A., White R., 2013, MNRAS, 434, 2289
  * Parizot & Drury (1999) Parizot E., Drury L., 1999, A&A, 349, 673
  * Parizot et al. (2004) Parizot E., Marcowith A., van der Swaluw E., Bykov A. M., Tatischeff V., 2004, A&A, 424, 747
  * Phan (2020) Phan V. H. M., 2020, PhD thesis, http://www.theses.fr/2020UNIP7070
  * Ptuskin et al. (2003) Ptuskin V. S., Jones F. C., Moskalenko I. V., Zirakashvili V. N., 2003, in International Cosmic Ray Conference. p. 1929 (arXiv:astro-ph/0306366)
  * Ptuskin et al. (2005) Ptuskin V., Moskalenko I., Jones F., Strong A., Mashnik S., 2005, Advances in Space Research, 35, 162
  * Ptuskin et al. (2006) Ptuskin V. S., Moskalenko I. V., Jones F. C., Strong A. W., Zirakashvili V. N., 2006, ApJ, 642, 902
  * Ptuskin et al. (2017) Ptuskin V. S., Zirakashvili V. N., Seo E. C., 2017, Bulletin of the Russian Academy of Sciences: Physics, 81, 437
  * Rauw et al. (2007) Rauw G., Manfroid J., Gosset E., Nazé Y., Sana H., De Becker M., Foellmi C., Moffat A. F. J., 2007, A&A, 463, 981
  * Salpeter (1955) Salpeter E. E., 1955, ApJ, 121, 161
  * Seo et al. (2018) Seo J., Kang H., Ryu D., 2018, Journal of Korean Astronomical Society, 51, 37
  * Silsbee & Ivlev (2019) Silsbee K., Ivlev A. V., 2019, ApJ, 879, 14
  * Smith (2006) Smith N., 2006, MNRAS, 367, 763
  * Smith & Brooks (2007) Smith N., Brooks K. J., 2007, MNRAS, 379, 1279
  * Streimatter et al. (1985) Streimatter R. E., Balasubrahmanyan V. K., Protheroe R. J., Ormes J. F., 1985, A&A, 143, 249
  * Sukhbold et al. (2016) Sukhbold T., Ertl T., Woosley S. E., Brown J. M., Janka H. T., 2016, ApJ, 821, 38
  * Tatischeff & Gabici (2018) Tatischeff V., Gabici S., 2018, Annual Review of Nuclear and Particle Science, 68, 377
  * Tatischeff et al. (2021) Tatischeff V., Raymond J. C., Duprat J., Gabici S., Recchia S., 2021, MNRAS, 508, 1321
  * Thornbury & Drury (2014) Thornbury A., Drury L. O., 2014, MNRAS, 442, 3010
  * Tibaldo et al. (2021) Tibaldo L., Gaggero D., Martin P., 2021, Universe, 7, 141
  * Tolksdorf et al. (2019) Tolksdorf T., Grenier I. A., Joubaud T., Schlickeiser R., 2019, ApJ, 879, 66
  * Verma et al. (1996) Verma M. K., Roberts D. A., Goldstein M. L., Ghosh S., Stribling W. T., 1996, J. Geophys. Res., 101, 21619
  * Vieu et al. (2020) Vieu T., Gabici S., Tatischeff V., 2020, MNRAS, 494, 3166
  * Vieu et al. (2021) Vieu T., Gabici S., Tatischeff V., 2021, preprint (arxiv:2109.10639),
  * Voelk & Biermann (1988) Voelk H. J., Biermann P. L., 1988, ApJ, 333, L65
  * Voss et al. (2010) Voss R., Diehl R., Vink J. S., Hartmann D. H., 2010, A&A, 520, A51
  * Weaver et al. (1977) Weaver R., McCray R., Castor J., Shapiro P., Moore R., 1977, ApJ, 218, 377
  * Yadav et al. (2017) Yadav N., Mukherjee D., Sharma P., Nath B. B., 2017, MNRAS, 465, 1720
  * Yang et al. (2018) Yang R.-z., de Oña Wilhelmi E., Aharonian F., 2018, A&A, 611, A77
  * Zhou & Matthaeus (1990) Zhou Y., Matthaeus W. H., 1990, J. Geophys. Res., 95, 14881

##  Appendix A Notations

Table 2 provides an exhaustive list of the notations used throughout the
present paper.

CLUSTER AND SUPERBUBBLE |  |
---|---|---
Initial number of massive ($M>8M_{\odot}$) stars | $N_{*}$ | $100-1000$
Lifetime of the superbubble | $\mathcal{T}_{\rm SB}$ | 35 Myr
Time-dependent mechanical power | $\mathcal{P}_{\rm tot}$ |  $10^{38}-10^{39}$ erg/s
Average mechanical power (winds and supernovae) | $L_{*}$ |  $10^{38}-10^{39}$ erg/s
Fraction of the mechanical power working on the shell | $\xi_{b}$ | 22%
Radius of the superbubble | $R_{b}$ |  $50-150$ pc
Superbubble radius normalised per star at 1 Myr | $R_{0}$ | 10 pc
Volume of the superbubble | $V_{\rm SB}$ |
Density of the interstellar medium around the superbubble | $n_{\rm ISM}$ |  $10-100$ cm-3
Interior number density | $n$ |  $0.01-0.1$ cm-3
Interior density | $\rho$ |  $0.2-2\times 10^{25}$ g cm-3
Interior temperature | $T$ |  $5\times 10^{6}$ K
MASSIVE STAR WINDS |  |
Wind mechanical power | $\mathcal{P}_{w}$ |  $10^{37}-10^{38}$ erg/s
Radius of the wind termination shock | $R_{s}$ |  $1-10$ pc
Velocity of the wind outflows | $V_{w}$ | 2000 km/s
Injection momentum at winds | $p_{0}$ | 5 MeV/c
Maximum momentum at winds | $p_{w}$ |  105 GeV/c
Wind acceleration efficiency | $\eta_{w}$ | 10%
SUPERNOVA REMNANTS (SNR) |  |
Supernova energy | $E_{\rm SN}$ |  $10^{51}$ erg
Mass of supernova ejecta | $M_{e}$ | $10~{}M_{\odot}$
Initial velocity of supernova shocks | $V_{\rm SN}$ | 3000 km/s
Velocity of SNR starting the Sedov-Taylor phase | $u_{0}$ | 1500 km/s
“Equivalent” volume spanned by “stationary” SNR | $\mathcal{V}_{\rm SNR}$ |  1% $V_{\rm SB}$
Fraction of cosmic rays reaccelerated at SNR | $\chi$ | 1%
Injection parameter of SNR shocks | $\xi$ | 3.5
Maximum momentum in SNR shocks | $p_{\rm max}$ | 1 PeV/c
TURBULENCE |  |
Source of the turbulence at the largest turbulent scale | $S$ |
Energy spectrum of the turbulence | $W$ |
Efficiency of turbulence generation | $\eta_{T}$ |  $1-30$%
Relaxation time of the turbulence generated by SNR | $\tau_{T}$ | 1 Myr
Largest turbulent scale | $\lambda$ | 10 pc
Smallest turbulent wavenumber | $k_{0}$ |
Non-thermal turbulence damping rate | $\Gamma$ |
Large scale magnetic field | $B_{0}$ | 10 µG
Random component of the magnetic field | $\delta B$ |  $1-10$ µG
Total magnetic field | $B$ |  $10-20$ µG
Turbulence level | $\eta$ |  $1-10$%
Alfvén velocity | $v_{A}$ |  $10-100$ km/s
Root-mean-square velocity | $\delta u$ |  $10-100$ km/s
PARTICLE TRANSPORT |  |
Particle distribution function | $f(x,p,t)$ |
Particle spectrum | $n(p,t)$ |
Particle velocity | $v$ |  $0.3-3\times 10^{8}$ m/s
Particle kinetic energy | $\epsilon(p)$ |  10 MeV $-$ 10 PeV
Spatial diffusion coefficient | $D_{x}$ |  $10^{27-28}$ cm2/s at 1 GeV/c
Momentum diffusion coefficient | $D_{p}$ |  1 (GeV/c)2/Myr at 1 GeV/c
Escape time from the superbubble | $\tau$ | 1 Myr at 1 GeV/c
Flux of particles coming from the interstellar medium | $\phi_{\rm ISM}$ |
Flux of particles escaping in the interstellar medium | $\Phi(E)$ |
Table 2: Notations used in the present work. The third colum provides typical
values. The input parameters of the model are highlighted with a star.