# Massive star cluster origin for the galactic cosmic ray population at very-
high energies

T. Vieu,1 B. Reville,1
1Max-Planck-Institut für Kernphysik, Saupfercheckweg 1, D-69117 Heidelberg,
Germany E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We demonstrate that supernova remnant (SNR) shocks embedded within massive
star clusters can reproduce both the cosmic-ray proton and all-particle
spectra measured in the vicinity of the Earth up to hundreds of peta-
electronvolts (PeV). We model two classes of massive star clusters. The first
population are “loose clusters” which do not power a collective wind
termination shock. SNR shocks then expand in a low-density and weakly
magnetised medium, and this population mainly contributes up to the “knee” of
the CR spectrum around 1 PeV. The second population are young compact
clusters, which are powerful and compact enough to sustain a collective wind
outflow. SNR shocks then expand from the cluster into the strongly magnetised
wind and accelerate nuclei up to ultra-high energies. This population,
representing only about 15% of all Galactic massive star clusters,
nevertheless dominates the spectrum between $\sim$ 1 and 100 PeV. While these
two components alone can reproduce the shape of the CR spectrum up to hundreds
of PeV, adding a light sub-ankle extragalactic component motivated by
composition and anisotropy measurements, allows to reproduce the spectrum up
to the highest energies. Fitting parameters are systematically linked to
physical variables whose values are in line with theoretical expectations.

###### keywords:

cosmic rays – acceleration of particles – open clusters and associations:
general – shock waves

††pubyear: 2022††pagerange: Massive star cluster origin for the galactic
cosmic ray population at very-high energies–Massive star cluster origin for
the galactic cosmic ray population at very-high energies

##  1 Introduction

Several arguments suggest that the transition between galactic and
extragalactic cosmic rays (CR) begins around several hundreds of PeV, with a
subdominant extragalactic component below 1 EeV, which becomes dominant at the
“ankle” at about 3 EeV. For example, measurements of the CR composition by
several experiments, despite uncertainties in the modelling of hadronic air
showers, show that the mean logarithmic mass of the detected particles
increases with energy beyond about 1 PeV, peaks at several hundreds of PeV,
and then decreases up to the ankle, where it starts to increase again (see
Figure 9 in Supanitsky, 2022). This is consistent with a scenario in which CRs
are accelerated within galactic sources with successive charge dependent cut-
off, resulting in an excess of heavy elements around the maximum energy of
these galactic sources, before a light sub-ankle component starts to take
over. This picture is supported by measurements on the dipole anisotropy,
which show that the phase of the anisotropy starts to smoothly increase
towards the Galactic Anti-centre beyond 100 PeV (see Figure 9 in Kachelrieß &
Semikoz, 2019). The modest amplitude of the anisotropy in these energy bands
disfavours a galactic origin for the light CR component in this energy range
(Giacinti et al., 2012). Galactic sources thus struggle to explain
simultaneously composition and anisotropy measurements beyond 100 PeV, while
the measurements are consistent with an extragalactic scenario. Moreover, at
hundreds of PeV the Larmor radius of CRs in galactic magnetic fields becomes
comparable to the scale-height of the Galactic Disk, which further disfavours
a galactic origin at higher energies due to confinement issues. Well below the
ankle, it has been pointed out (e.g. Parizot, 2014) that a concave feature
separating two extended power-laws, such as the “knee” observed in the CR
spectrum at about 3 PeV, requires fine-tuning to be accounted for in a two-
component model, with the second component being subdominant yet non-
negligible. If the two components are completely unrelated, as is the case if
one assumes that the CRs below/above the knee are of galactic/extragalactic
origin, such tuning is unlikely and a smooth spectrum is not expected in
principle. On the other hand, the second knee at about 0.5 EeV, which is a
less pronounced concave feature, suffers less from this issue as the steep
power-law component extends beyond it over at most a decade, and therefore
this feature could mark the turn-over of the heavy galactic CRs, a few
hundreds of PeV before the light extragalactic component takes over (e.g.
Thoudam et al., 2016), especially if the latter is not a perfect power-law, as
suggested by a number of models (e.g. Unger et al., 2015; Globus et al., 2015;
Kachelrieß et al., 2017; Berezinsky et al., 2006).

Setting the transition between galactic and extragalactic CRs at hundreds of
PeV implies that at least a subset of the galactic sources – on the order of
20% (Axford, 1994; Parizot, 2014) – must be able to accelerate heavy nuclei up
to hundreds of PeV, i.e. protons beyond 10 PeV. Identifying these sources is
not a straightforward task and there is no consensus on the matter at the
present time. Indeed, it has long been suspected that supernova remnant (SNR)
shocks expanding in the warm interstellar medium (ISM) struggle to reach PeV
energies (Lagage & Cesarsky, 1983), even when considering CR-induced
instabilities in the dense wind blown by the progenitor star (Bell et al.,
2013), unless the SNR shock is exceptionally fast or exhibit extreme transport
properties (e.g. Voelk & Biermann, 1988; Bell & Lucek, 2001; Ptuskin et al.,
2010; Takamoto & Kirk, 2015). However, even in the most optimistic scenarios,
the maximum energy is an order of magnitude lower than what is required to
account for the GCR spectrum up to hundreds of PeV.

A number of alternative hypotheses accounting for the GCR spectrum in the
range $1-100$ PeV have been proposed, including, e.g., the re-acceleration of
particles by the galactic wind termination shock (Jokipii & Morfill, 1987;
Zirakashvili & Völk, 2006; Merten et al., 2018), the re-acceleration by
multiple scatterings in the Fermi bubbles (Cheng et al., 2012), the collective
effects in superbubbles (Bykov & Toptygin, 2001; Bykov et al., 2018). It is as
well possible that a different class of sources takes over beyond the knee,
e.g. the Galactic Centre (HESS Collaboration et al., 2016), microquasars
(Bosch-Ramon et al., 2005), hypernovae (Sveshnikova, 2003), or young massive
stars (Biermann & Cassinelli, 1993) and young stellar clusters blowing a
collective wind (Aharonian et al., 2019; Vieu et al., 2022c). The interest
towards the latter has been rapidly raising in the recent years since the
discovery that these sources are efficient gamma-ray emitters up to hundreds
of TeV (e.g. Abeysekara et al., 2021; Cao et al., 2021; Aharonian et al.,
2022).

It was realised early on (e.g. Higdon & Lingenfelter, 2005) that most massive
stars (and hence SN progenitors) form in clusters. The powerful stellar
feedback of the young early-type stars subsequently carve hot “superbubbles”
within the parent molecular cloud (Weaver et al., 1977; McCray & Kafatos,
1987). Efficient particle acceleration is expected to take place within these
environments, by the stellar winds (Montmerle, 1979; Cesarsky & Montmerle,
1983), by the turbulence (Bykov, 2001; Klepach et al., 2000), or by supernova
explosions (Parizot et al., 2004; Ferrand & Marcowith, 2010; Vieu et al.,
2022b). Several arguments support a superbubble origin of CRs; in particular,
the enrichment of the accelerated medium by the winds of Wolf-Rayet stars has
been shown to resolve the ${}^{22}\rm{Ne}/^{20}\rm{Ne}$ anomaly of the GCRs
(Gupta et al., 2020).

In a recent paper (Vieu et al., 2022c), we assessed the maximum proton
energies achievable under various particle acceleration scenarios around
massive star clusters. In particular, we have shown that, similarly to the
case of isolated SNR shocks, it is not straightforward to reach PeV energies
when realistic environmental parameters are considered. Nevertheless, we
highlighted a promising possibility, that is to accelerate particles around a
SNR shock launched from a compact cluster into the collective large-scale
($\sim 10$ pc) wind blown by the massive stars. Large magnetic fields are
amplified in the core of the cluster, potentially up to hundreds of µG, and
although the level of magnetic turbulence decreases in the powerful collective
wind, the average magnetic field in the region swept up by the SNR shock is
still expected to be about one order of magnitude larger than the average
interstellar value, which can facilitate effective particle acceleration.
Under reasonable assumptions on the efficiency of turbulence generation inside
the cluster, we computed the maximum energy of the particles as
$E_{max}\lesssim 5Zu_{5}$ PeV for standard cluster parameters, where $u_{5}$
is the initial SN shock velocity in units of 5000 km/s. This suggests that, in
favourable conditions, a fast ($>10$ 000 km/s) SNR shock can accelerate
protons beyond 10 PeV, and the fastest ($>30$ 000 km/s) shocks could almost
accelerate iron up to 1 EeV, therefore meeting the required bounds to explain
the GCR spectrum up to hundreds of PeV.

In the present work, we aim to compute the GCR spectrum including this
component due to particle acceleration around this special population of SNR
shocks. We start in Section 2 by a census on the galactic star-cluster
population, which is used to estimate the frequency of these SN explosions.
Section 3 is devoted to a description of a two-component model, including the
contribution of SNR shocks expanding in loose associations of massive stars,
where no collective wind termination can form, as well as the contribution of
clustered SNR shocks. We eventually compute the overall contribution of
massive star clusters to the GCRs in Section 4, where we discuss the agreement
between our model and measurements made in the vicinity of the Earth between 1
TeV and 1 EeV. We conclude in Section 5.

##  2 The massive star cluster population of the Galaxy

In this Section we review and categorise the galactic cluster census presented
in the recent analysis of the Gaia Data Release 2 (DR2) (Cantat-Gaudin et al.,
2020).

###  2.1 Definitions

Because star clusters are at the crossroads of several sub-fields of
astrophysics, e.g. cosmic-ray physics and astrometry, the nomenclature is not
well-defined. In the present work, we use the following designations, which
categorise massive star clusters in a way which will prove useful in what
follows.

Star cluster: a star cluster is defined as an association of stars over a
scale under a few tens of pc, which is not necessarily gravitationally
bounded. These clusters are referred to as open clusters in the astrometry
community, as opposed to globular clusters which are much older and extend to
much larger scales. We do not presently discuss globular clusters as they are
not expected to represent a sizeable fraction of the GCR sources.

Massive star cluster (MSC): a massive star cluster is a star cluster which
contains SN progenitors, i.e. stars with zero-age main-sequence (ZAMS) mass
larger than about 8$M_{\odot}$ (e.g. Poelarends et al., 2008; Smartt et al.,
2009)111In the following we disregard the fact that some stars with ZAMS mass
larger than $8M_{\odot}$ might not explode (e.g. Sukhbold et al., 2016)..
Importantly, this excludes all clusters which are more than about 40 Myr old,
as all SN progenitors are expected to explode within 40 Myr, assuming coeval
star formation (Limongi & Chieffi, 2006).

Young massive star cluster (YMSC): a young MSC is a MSC which still contains
stars expected to undergo a transition to the WR phase in the end of their
life. Assuming roughly that stars with ZAMS mass larger than 20$M_{\odot}$ are
expected to become WR stars (Crowther, 2007; Sander et al., 2012), young
clusters are clusters which are less than about 10 Myr old, assuming coeval
star formation.

Evolved massive star cluster: in contrast to YMSCs, evolved MSCs are defined
as MSCs which do not contain WR stars anymore (i.e. which are more than $10$
Myr old). Their wind mechanical power is then much smaller than that of YMSCs,
such that their energy budget is largely dominated by SN explosions (Vieu et
al., 2022b).

Compact massive star cluster: in the present work, we consider that a MSC is
compact if its half-mass radius is lower than $3(t_{\rm
c}/\mathrm{Myr})^{2/5}$ pc, where $t_{\rm c}$ is the cluster age. The physical
motivation lying behind this somewhat arbitrary threshold is clarified below.

Extended massive star cluster: in contrast to compact MSCs, extended MSCs are
MSCs whose stellar distribution extends beyond $3(t_{\rm
c}/\mathrm{Myr})^{2/5}$ pc.

Wind-blowing cluster (WBC): if a YMSC is compact enough, the powerful winds of
the young embedded early-type and Wolf-Rayet (WR) stars are expected to merge
into a collective outflow. This collective wind terminates at the wind
termination shock (WTS), whose radius can be estimated theoretically as
(Weaver et al., 1977):

$R_{\rm WTS}\approx N_{*}^{3/10}\left(\frac{n_{\rm ext}}{100~{}{\rm
cm}^{-3}}\right)^{-3/10}\left(\frac{t_{\rm c}}{1~{}{\rm
Myr}}\right)^{2/5}~{}{\rm pc}\,,$ (1)

where $N_{*}\approx 10-1000$ is the total number of massive ($M\gtrsim
8M_{\odot}$) stars in the cluster, $n_{\rm ext}\approx 1-100$ cm-3 is the
density of the medium beyond the superbubble shell, $t_{\rm c}=1-10$ Myr is
the age of the (young) cluster. In order to express the WTS radius as a
function of the number of massive stars $N_{*}$, we assumed an average
mechanical power per star of $10^{36}$ erg/s, which, as will be discussed in
details in Section 4.1, is only expected in the case of a young cluster. If
the clustered massive stars are distributed over a region smaller than $R_{\rm
WTS}$, then the collective wind can form. Otherwise, the theoretical model
breaks down and the individual winds are not expected to merge in a collective
outflow. We conveniently define wind-blowing clusters as the subset of YMSCs
which are compact enough to power a collective wind beyond the boundary of the
cluster. Eq. 1 provides an estimate of $R_{\rm WTS}$ ranging from about
$2(t_{\rm c}/\mathrm{Myr})^{2/5}$ pc for small clusters to about $20(t_{\rm
c}/\mathrm{Myr})^{2/5}$ pc for the most powerful star clusters evolving in a
rather dilute environment. The strong dependency of the WTS size on the
cluster and environmental parameters, as well as the difficulty to link the
theoretical “cluster boundary” to the observable “half-mass radius”, or the
difficulty to estimate the total number of massive members from observations,
make it challenging to set a simple compactness threshold for the formation of
a WTS. In the following, in order to constrain the fraction of WBCs, we will
assume a theoretical WTS shock radius of $3(t_{\rm c}/\mathrm{Myr})^{2/5}$ pc,
as the average cluster contains a few tens of massive stars and evolves in a
dense external medium. We consider this to be a conservative value and can be
understood as an upper bound, which only has a meaning when averaged over a
large number of MSCs.

Loose cluster (LC): in contrast to the population of WBCs, we define loose
massive star clusters as the MSCs which are unable to power a collective
outflow, which includes both extended and evolved clusters.

###  2.2 Massive cluster population as inferred from Gaia DR2

Figure 1: Left: Radial distribution (centred at the Sun) of the MSCs extracted
from the open cluster catalogue of Cantat-Gaudin et al. (2020). Right: MSCs
located within 3 kpc, classified in two categories: wind-blowing clusters (red
disks) and loose clusters (grey disks). The size of the disk is proportional
to the number of members detected in the cluster. The histogram in the
background shows the age distribution of the clusters.

Cantat-Gaudin et al. (2020) have performed a systematic analysis of the
clusters identified in the Gaia DR2. They used artificial neural network
techniques in order to characterise 1867 star clusters distributed within a
disk of radius about 6 kpc around the Sun. They estimated, among other
properties, the age of the clusters, their half-mass radii and their stellar
content with membership probability larger than 70%. As particle acceleration
is expected to be predominantly powered by massive stars, we restrict the
sample to the massive star clusters, which are less than 40 Myr old. This
represents a total of 430 MSCs. As shown in the left panel of Figure 1, the
distribution of these clusters as a function of distance from the Sun drops
beyond 3 kpc, due to the difficulty of inferring reliable cluster parameters
from reddened and blurred colour-magnitude diagrams observed in regions of
high extinction. In order to extract a representative sample of the cluster
distribution in the Galaxy, we disregard in the following the tail of the
distribution beyond 3 kpc. We are then left with 320 “nearby” MSCs located
within and in-between the Local, Perseus and Sagittarius arms.

Knowing the half-mass radii and the ages, we can classify the clusters in our
sample according to the definitions previously set. The result is shown on the
right panel of Figure 1. Although most of the clusters are compact, with a
half-mass radius $\lesssim 5$ pc as seen in Figure 1, only 50 are young enough
to be wind-blowing clusters. These numbers are summarised in Table 1.

The histogram showing the distribution of MSCs as a function of their age,
which is plotted in the background of the right panel of Figure 1, does not
display any correlation. This is expected since the lifetime of the massive
members in a star cluster is much smaller than the cluster disruption time
(Krumholz et al., 2019).

| Young | Evolved
---|---|---
Compact | 16% | 80%
Extended | 1% | 3%
Table 1: Representation of massive star clusters in the catalogue of Cantat-
Gaudin et al. (2020) located within 3 kpc from the Sun, classified according
to the definitions established in Section 2.1.

###  2.3 On the difficulty of estimating the number of massive star clusters
in our Galaxy

One can try to estimate the number of MSCs expected in the Galaxy on
theoretical grounds, considering an average star formation rate of
$2M_{\odot}$/yr (Molinari et al., 2014) and a clustering fraction of 80%
(Higdon & Lingenfelter, 2005). This results in a total stellar mass of about
$M_{\rm tot}=3\times 10^{7}M_{\odot}$ over 20 Myr – the average MSC age. The
sample of observed MSCs are distributed according to an IMF scaling as $M_{\rm
cl}^{-2}$ (Krumholz et al., 2019), such that the total number of MSCs in the
Galaxy reads:

$N_{\rm cl}=\frac{M_{\rm tot}/M_{\rm min}}{\ln\left(M_{\rm max}/M_{\rm
min}\right)}\,,$ (2)

where $M_{\rm min}$ ($M_{\rm max}$) are the lower (higher) ends of the MSC
distribution in the Milky Way. The maximum mass of a star cluster (typically
$10^{5}M_{\odot}$) only enters the logarithm in Eq. 2 and is therefore not a
critical parameter, in contrast to the minimal mass, which must be precisely
specified, although it is not well-defined physically. Indeed, the $M_{\rm
cl}^{-2}$ distribution extends down to very low masses, e.g. $10^{2}M_{\odot}$
(Lamb et al., 2010), but it is not relevant to talk about a “massive” star
cluster if it only hosts one or two massive members. If $M_{\rm min}$ varies
in the range $10^{2}M_{\odot}-10^{3}M_{\odot}$, one obtains $N_{\rm cl}\approx
6000-45000$, i.e., assuming that these are distributed over a galactic disk of
radius 15 kpc and that the local star formation rate is half that of the
Galaxy average (Miller & Scalo, 1979), we should expect 100 to 1000 MSCs in
the 3 kpc vicinity of the Sun. The number of 320 clusters catalogued in
Cantat-Gaudin et al. (2020) lies within this range and we could extrapolate it
to a total of about 16k clusters in the Galaxy, which would be obtained from
Eq. 2 assuming a minimum mass of $450M_{\odot}$ (6 massive members in the
cluster). This appears to be a reasonable estimate a posteriori, though the
completeness of the local cluster catalogue is far from certain. Finally,
given that the star formation rate is inhomogeneous on the scale of the entire
Galaxy, one might wonder if these extrapolations make any sense at all. In the
following we will alleviate this issue by considering, in the computation of
the overall normalisation of the GCR flux, the supernova rate instead of the
total number of star clusters.

##  3 Spectrum of particles accelerated by clustered supernovae

###  3.1 First population: supernovae in loose clusters

By definition, “loose” clusters do not sustain a collective outflow, as they
are either too old (in which case there is not enough stellar power to form a
collective WTS) or too extended (in which case the winds do not interact to
produce a collective outflow). Individual stellar winds, as well as SNRs,
deposit energy in the environment, which heats an expanding cavity filled by
turbulent MHD outflows. Because the turbulent magnetic fields are diluted in
this large cavity (sometimes referred to as a “superbubble”), the average
magnetic field probed by SNR shocks is not expected to exceed a few µG and the
maximum energy of the particles is estimated as (Lagage & Cesarsky, 1983;
Hillas, 2005):

$\displaystyle E_{\rm max}\approx 0.1\frac{ZB}{2~{}\text{\textmu
G}}u_{5}^{1/3}\left(\frac{E_{\rm
SN}}{10^{51}~{}\mathrm{erg/s}}\right)^{1/3}\left(\frac{n}{0.01{{\rm
cm}^{-3}}}\right)^{-\frac{1}{3}}\text{~{}PeV}\,,$ (3)

where $B$ is the average magnetic field, $n$ is the density within the
superbubble cavity and $E_{\rm SN}$ is the explosion energy. This computation
assumes a constant shock velocity in the free expansion phase. A similar
estimate is obtained when considering a decelerating shock, which is expected
to be the case after the collision with the (individual) WTS of the progenitor
star. Only the dependency on velocity changes, e.g. $E_{\rm max}\propto
u_{5}^{2/3}$ for $R\propto t^{4/7}$, but this is not expected to affect the
conclusions of this work. We set, for a loose cluster SN (LCSN) remnant shock
expanding into a typical superbubble environment:

$E_{\rm max,LCSN}\approx 0.1Zu_{5}^{1/3}~{}{\rm PeV}\,,$ (4)

where the numerical factor can be taken as a parameter of the model.

Eventually we model the escaping differential energy spectrum for a given
nucleus of charge $Z$ accelerated from a LCSN remnant shock as follows:

$S_{\rm LCSN}=\mathcal{A}_{\rm LC}\frac{E_{\rm
SN}}{(p_{0}c)^{2}}f_{Z}\left(\frac{p}{p_{0}}\right)^{-2-\Delta
s}\exp\left(\frac{-pc}{E_{\rm max,LCSN}}\right)\,,$ (5)

where $\mathcal{A}_{\rm LC}$ is a normalisation constant computed numerically
such that $\int_{p_{0}}\mathrm{d}pc\,ES_{\rm LCSN}\equiv f_{Z}E_{\rm
SN}$222$S_{\rm LCSN}$ is related to the particle distribution function $f$ as
$S_{\rm LCSN}=4\pi p^{2}f$., with $p_{0}$ the injection momentum. $f_{Z}$ is
therefore an acceleration efficiency which depends on the considered nucleus,
in principle as a function of the composition of the medium and injection
efficiency. In practice, we set the proton acceleration efficiency
$f_{1}\approx 0.1$ as an input parameter. The parametrization of the
acceleration efficiencies $f_{Z}$ of heavier nuclei is then set a posteriori
to match the composition measured at 1 TeV in the local GCR spectrum, from the
fluxes tabulated in Hörandel (2003).

The parameter $\Delta s$ allows for spectral slope deviation from the
canonical $p^{-4}$ scaling. This is motivated partly by gamma-ray observations
displaying spectral indices typically around $2.1-2.5$ (see e.g. Table 1 in
Caprioli, 2011). Although there is presently no convincing explanation of the
discrepancy between these observations and the prediction of the standard
diffusive shock acceleration theory, a number of arguments have been proposed,
for example: the effects of nonlinear processes around the SNR shock (e.g.
Kirk et al., 1996; Caprioli, 2012; Vieu et al., 2022a; Bell et al., 2019),
magnetic-field obliquity (Bell et al., 2011; Xu & Lazarian, 2022), the
weakening of the SNR shock after its collision with the wind termination shock
of the progenitor star (e.g. Das et al., 2022), or the effect of the time-
dependent particle leakage from the system (e.g. Celli et al., 2019). Within
superbubbles, such a steepening might be even aggravated during the
propagation of the particles in the turbulent environment, possibly up to very
high energies (Bykov, 2001; Ferrand & Marcowith, 2010).

###  3.2 Second population: supernovae in wind-blowing clusters

By definition, “wind-blowing clusters” (WBC) are young and compact enough to
power a collective wind. In this case, the energy deposited by stellar winds
and SNRs is expected to excite strong turbulence within the cluster. The
turbulent energy decreases in the free-wind region which extends beyond the
boundary of the cluster (e.g. as $1/r^{2}$). Therefore, the magnetic field
induced within the cluster stays confined within the cluster and can in
principle reach large values, up to hundreds of µG (Vieu et al., 2022c). This
phenomenological argument is supported by recent 3D MHD simulations of the
inner core of a WBC (Badmaev et al., 2022), which showed that even without
resolving the small-scale turbulence, the kinetic stellar energy is
efficiently converted to magnetic energy via the Cranfill effect, thereby
creating magnetic fields up to several hundreds of µG in filamentary
structures. Although these simulations do not probe the effect of SN
explosions nor the properties of the collective outflow beyond the boundaries
of the cluster, they do showcase that efficient magnetic field amplification
is to be expected in and around WBCs. A SNR shock expanding into these large
magnetic fields could accelerate protons up to a maximum energy which is
limited by the available time. The acceleration rate, assuming diffusion at
the Bohm limit, is $\mathrm{d}p/\mathrm{d}t\approx 3qB/(\alpha c)\dot{R}^{2}$,
where $\alpha\approx 9$ if one assumes that the upstream magnetic field $B$ is
compressed by a factor $\sqrt{11}$ across the shock. Because the magnetic
field decreases in the collective wind beyond the boundary of the cluster
($B\propto 1/r$), the acceleration rate is higher around a SNR shock expanding
inside the cluster than around a SNR shock launched close to the edge and
expanding in the cluster wind. However the fate of SNRs expanding deep inside
the cluster is unclear. Collisions with the strong winds of the neighbouring
massive stars are expected to hinder, at least partly, the propagation of the
shock, eventually converting part of the SNR shock pressure into MHD
turbulence. Although a system of converging shocks might be in itself an
efficient accelerator of very-high energy CRs (Bykov et al., 2015), it is yet
unclear if the increase of the maximum energy compared to a single shock is a
generic outcome of this class of systems (Vieu et al., 2020). In the
following, we therefore assume that such collisions between shocks in compact
clusters do not leave noticeable spectral features at ultra-high energies.

If a SN explodes deep within the cluster, its remnant shock typically reaches
the edge of the cluster at the time of the transition to the Sedov-Taylor
phase, which implies that particles at the highest energy are not further
accelerated in the cluster wind, but on the contrary suffer adiabatic losses
when advected from the cluster edge ($\sim 1$ pc) to the cluster WTS ($\sim
10-20$ pc). This shifts the spectrum back to lower energies, such that, in the
end, the maximum CR energy of the particles escaping the spectrum is typically
higher if the SN shock is launched close to the cluster edge and expands
mostly in the collective wind. As a conservative approach, we only consider
the latter in the following. The fraction of SN which produce shocks expanding
far in the wind can be estimated as follows. Assuming that $N_{*}$ massive
stars are homogeneously distributed in the cluster, the distance between them
reads $L=(2/3)^{1/3}R_{\rm c}N_{*}^{-1/3}$. About $2\pi R_{\rm c}/L$ stars lie
at the edge of the cluster, which represents a fraction $\eta_{\rm SN}\approx
7N_{*}^{-2/3}\approx 0.3-0.5$ of all SN progenitors. The parameter $\eta_{\rm
SN}$ can be taken as an efficiency which globally accounts for the
aforementioned uncertainties about the fate of SN shocks launched within
compact star clusters.

Let us now estimate the maximum energy of particles accelerated around a SNR
shock which expands in the collective wind. Assuming a $1/r$ magnetic field
profile in this region, we get, from the acceleration rate:

$E_{\rm max}\approx 3qB_{\rm c}R_{\rm c}/(\alpha c)u_{\rm SN}\ln\left(R_{\rm
max}/R_{\rm c}\right)\,,$ (6)

where $B_{c}$ is the value of the magnetic field at the edge of the cluster,
$R_{\rm c}$ is the radius of the cluster, $u_{\rm SN}$ is the initial SN shock
velocity, $R_{\rm max}$ is either the distance reached at the end of the free
expansion phase, or the radius of the WTS. For fast SN shocks (30 000 km/s)
launched close to the edge of the cluster, one can compute333We assume that
the energy released in each SNR is $\approx 10^{51}$ ergs $R_{\rm max}/R_{\rm
c}\approx 5$, while for slow SN shocks (5000 km/s), the shock collides with
the WTS while being still in the free expansion phase and $R_{\rm max}/R_{\rm
c}\sim 10$. Neglecting this small logarithmic correction, we obtain:

$E_{\rm max}\approx 3Z\left(\frac{B_{\rm c}}{100~{}\text{\textmu
G}}\right)\left(\frac{R_{\rm c}}{5~{}\text{pc}}\right)u_{5}~{}\mathrm{PeV}\,,$
(7)

where $u_{5}$ is again the initial shock velocity in units of 5000 km/s.
Taking into account the deceleration of the shock in the wind profile
($u\propto t^{-1/8}$, Finke & Dermer, 2012; Gaggero et al., 2018) leads to a
similar estimate, as detailed in Vieu et al. (2022c). We can generically
write:

$E_{\rm max,WBCSN}\approx Zu_{5}f_{\rm c}(R_{\rm c},N_{*},n_{\rm c},\eta_{\rm
T})~{}\text{PeV}\,,$ (8)

where $f_{\rm c}\sim O(1)$ is a function of the cluster properties (cluster
radius $R_{\rm c}$, number of massive stars $N_{*}$, density within the
cluster $n_{\rm c}$, and the fraction of the stellar mechanical power injected
into the MHD turbulence cascade $\eta_{\rm T}$). As we aim to obtain the
average contribution from all clusters in the Galaxy, and because the SN shock
velocity $u_{5}$ is not expected to correlate with the properties of the
cluster, the function $f_{\rm c}$ can be estimated on average. It is
furthermore not meaningful to compute it in detail as most of the cluster
parameters are uncertain, both from the observational and theoretical
perspectives. In particular, it is difficult to constrain the turbulence
generation efficiency $\eta_{\rm T}$ or the density outside of the superbubble
cavity. In the following we take the following estimate for an average
galactic cluster (assuming $\eta_{\rm T}=10\%$):

$E_{\rm max,WBCSN}\approx 4Zu_{5}~{}\text{PeV}\,,$ (9)

which can be taken as an input parameter of the model, justified by reasonable
values of the average cluster properties in Eq. 8 (details can be found in
Vieu et al., 2022c).

In the scenario considered here, the magnetic field is not amplified via
streaming instabilities but rather from a deposition of energy at the stellar
scale within the boundaries of the cluster. It is assumed that, if necessary,
CRs still sustain the scattering fields required to ensure Bohm diffusion.
These SNR shocks are therefore likely to produce hard spectra of high-energy
particles which we assume to match the standard $p^{-4}$ power-law achieved
around a strong shock, such that we can write:

$S_{\rm WBCSN}=\mathcal{A}_{\rm WBC}\frac{E_{\rm
SN}}{2}f_{Z}(pc)^{-2}\exp\left(\frac{-pc}{E_{\rm max,WBCSN}}\right)\,,$ (10)

where $\mathcal{A}_{\rm WBC}$ is a normalisation constant set such that
$\int_{p_{0}}\mathrm{d}pc\,ES_{\rm WBCSN}\equiv f_{Z}E_{\rm SN}/2$. Only about
half of the shock is launched in the collective wind region if the supernova
explodes close to the cluster edge, hence the additional factor of 1/2 in the
normalisation. The acceleration efficiency $f_{Z}$ is assumed to be equal to
that of the LCSN remnant shocks, even though the composition of the wind blown
by a WBC differs from that of a LC, due to the presence of WR stars. Figure 3
in Tatischeff et al. (2021) suggests that the contribution from the latter
does not make a dramatic difference and thus will not affect our results
substantially. A detailed study of the WBC plasma composition relative to that
of LCs, taking into account the material launched by the stellar population,
the material deposited by expanding SNR, and the mass-loading from the ISM, is
beyond the scope of this paper.

##  4 Galactic CR spectrum from massive star clusters

In this section we consider the entire population of galactic MSCs as analysed
in Section 2, in order to compute their contribution to the flux of particles
detected near the Earth.

###  4.1 Spectrum of particles accelerated by stellar winds

In the previous section we computed the contribution of the clustered SNR
shocks to the GCR population. However, massive stars blow powerful winds
during their lives, which create wind termination shocks (WTS), and it was
early realised that these could be favourable sites of particle acceleration
(e.g. Casse & Paul, 1980; Cesarsky & Montmerle, 1983). In this section we
consider the contribution of stellar winds in the model.

We first estimate the average mechanical power per massive star in young and
evolved clusters, which is critical to infer the normalisation of the wind
components. Stellar evolution models can be used to follow the evolution of a
population of clustered massive stars from the formation of the cluster
(usually one assumes a coeval population of stars) to the explosion of the
last SN, which by our definition marks the end of the MSC life. In the young
phase, the mechanical power of MSCs is dominated by the outflows of early O
and WR stars, while in the evolved phase, where by definition there is no
remaining WR stars, the power drops significantly. Assuming a Salpeter IMF and
using the fits of stellar evolution models provided by Seo et al. (2018) for
main sequence O and early B stars, we compute an average mechanical power per
massive star of $5\times 10^{32}$ erg/s for evolved clusters. By comparison,
each clustered supernova provides $10^{51}$ erg, which, averaged over a
cluster lifetime of 40 Myr, represents a power of $8\times 10^{35}$ erg/s. The
mechanical power input of stellar winds embedded in evolved clusters can
therefore be neglected.

Figure 2: Left: Mechanical power of WR stars as a function of their mass and
type. The blue, red and yellow dashed lines show the linear fit to the
individual types whilst the solid black line shows the average contribution
assuming an equal number of WC, WN and WNh stars. These fits take into account
the updated data from Sander et al. (2019); Hamann et al. (2019) (points), to
be compared with the transparent lines which fit the oldest data from Nugis &
Lamers (2000) (triangles). Right: Evolution of the mechanical power per
massive star as function of the age of an average massive star cluster,
assuming a coeval stellar population initially distributed according to a
Salpeter IMF.

To compute the mechanical power of massive stars embedded in young clusters,
it is critical to know the contribution from the Wolf-Rayet stars. Seo et al.
(2018) estimated this quantity by fitting observational data from Nugis &
Lamers (2000) (shown by the triangles in the left panel of Figure 2). However
the lack of observations of WR stars more massive than $25M_{\odot}$ and the
dispersion of the data make it difficult to extrapolate this analysis to high
stellar masses. In 2019, the population of galactic WN and WC stars was
revisited and complemented in light of the Gaia DR2 (Sander et al., 2019;
Hamann et al., 2019). The updated list of single WR stars contains WC stars up
to $30M_{\odot}$, WN stars up to $45M_{\odot}$ and WNh stars up to
$70M_{\odot}$, which makes the extrapolation of the data more meaningful. As
shown in the left panel of Figure 2, the updated relation between mechanical
power and star mass (solid black line) deviates significantly from that
obtained by fitting the data of Nugis & Lamers (2000) (transparent black
line), with a much shallower slope, which implies that high mass WR stars are
less powerful than previously thought, by up to one order of magnitude.
Including these revised relations in a time-dependent cluster model, we
compute, for young clusters following a Salpeter IMF, an average mechanical
power per massive star of $4\times 10^{35}$ erg/s if the fit is performed
ignoring the WR spectral type. The average increases to $7\times 10^{35}$
erg/s if we fit the three WR spectral types (WC, WN, WNh) independently and
then assume that WR stars are equally distributed between these three
types444This is likely providing an upper bound as it does not account for the
fact that the WC type is actually a stellar evolution phase subsequent to the
WN phase, hence the average mass of WC stars is expected to be significantly
smaller than that of WN stars, as suggested by the observational data (Figure
2).. In the following, we will assume an average mechanical power per massive
star embedded in a young cluster of $5\times 10^{35}$ erg/s555In comparison,
the relations found by fitting the outdated data from Nugis & Lamers (2000)
provide a significantly higher average power of $1.7\times 10^{36}$ erg/s..
Integrated over several Myr, this power provides an energy input comparable to
that of a supernova, hence the contribution of the young cluster stellar winds
to the galactic CRs cannot be disregarded.

In order to compute the normalisation of the wind component, we also need to
estimate the total number of massive stars in the Galaxy. This can be done in
two independent ways, either considering the birth or death rate. First, the
average star formation rate of $2M_{\odot}$/yr in the Milky Way is expected to
generate a population of about 450k massive stars, assuming that massive stars
are distributed according to a Salpeter initial mass function and live on
average for 17 Myr. Second, to sustain a rate of core-collapse supernova of 2
per century, there must be at least 340k massive stars of average lifetime 17
Myr. As a compromise we assume the total number of massive stars to be 400k.
Then, knowing the fraction of WBC as determined in Section 2, and assuming
that about $80-90$% of massive stars are clustered (Higdon & Lingenfelter,
2005; Lingenfelter, 2018), we can infer the total number of stars embedded in
WBCs. Finally, one estimates the total wind power of galactic WBCs as:

$P_{\rm w,WBC}=f_{\rm WBC}f_{\rm c}N_{\rm*,Gal}\bar{P}_{\rm w}\,,$ (11)

where $f_{\rm WBC}\approx 15\%$ is the fraction of WBCs in the galactic MSC
population (Section 2), $f_{\rm c}\approx 80-90\%$ is the probability that a
massive star is embedded in a MSC, $N_{\rm*,Gal}\approx 400k$ is the total
number of massive stars in the Galaxy, and $\bar{P}_{\rm w}\approx 5\times
10^{35}$ erg/s is the average mechanical power per massive star in a young
cluster as discussed above.

WBCs by definition blow a collective wind which terminates at a distance of
order 10 pc from the centre of the cluster. Particle acceleration is expected
to take place at the wind termination shock up to about 1 PV, for favourable
conditions (Gupta et al., 2020; Morlino et al., 2021; Vieu et al., 2022c). We
therefore model the spectrum of accelerated particles as follows:

$\phi_{\rm w}=\mathcal{A}_{\rm w}f_{Z}(pc)^{-2}\exp\left(\frac{-pc}{E_{\rm
max,WTS}}\right)\,,$ (12)

where the normalisation constant $\mathcal{A}_{\rm w}$ is set such that
$\int_{p_{0}}\mathrm{d}pc\,E\phi_{\rm w}\equiv f_{Z}P_{\rm w,WBC}$, and the
acceleration efficiency $f_{Z}$ is assumed to be the same as for the case of
strong SNR shocks. Gamma-ray observations of star clusters typically display
spectra which are slightly steeper than $E^{-2}$ (e.g. Yang et al., 2018;
Abeysekara et al., 2021; Aharonian et al., 2022), however the hard wind
component calculated from Eq. 12 will be already found in the following to
provide a subdominant contribution to the GCR flux – which can be anticipated
as the integrated power of stellar winds is about half that of SN explosions
–, such that the shape of the spectrum given by Eq. 12 is not expected to
impact our results.

###  4.2 Escaping spectra from a distribution of sources

The maximum energy achieved in SNR shocks depends on the initial velocity of
the shock (Eqs. 4 and 9). We therefore need to introduce a distribution of
velocities denoted $\xi$, and defined as:

$\mathrm{d}f_{\rm SN}=\xi(u)\mathrm{d}u\,,$ (13)

such that $f_{\rm SN}=\int\xi(u)\mathrm{d}u$ is the total rate of core-
collapse SN (CCSN) in our Galaxy. This distribution is difficult to constrain.
As discussed in Sveshnikova (2003), it is challenging to interpret the
distribution of SN absolute-magnitudes derived from observations, and it is
further difficult to constrain the initial velocity of the SNR shock. In the
following we assume that SN shock velocities range between 3000 km/s and 30
000 km/s (e.g. SN1987A), with a power-law distribution, e.g. $\xi(u)\propto
1/u$. The sensitivity of our results on the shape of the distribution $\xi$ is
discussed in the end of Section 4.4.

Finally, we gather the two SN components together with the wind component to
obtain the total flux released by MSCs in our Galaxy:

$\phi_{\rm c}(E,Z)=(1-f_{\rm WBC})f_{\rm c}\int\mathrm{d}u\,\xi(u)S_{\rm
LCSN}\\\ +f_{\rm WBC}f_{\rm c}\eta_{\rm SN}\int\mathrm{d}u\,\xi(u)S_{\rm
WCSN}+\phi_{\rm w}\,,$ (14)

where $f_{\rm WBC}\approx 15\%$ is the fraction of WBC in the galactic MSC
population obtained in Section 2, $f_{\rm c}\approx 80\%$ is the clustering
fraction of massive stars and, as discussed in Section 3.2, $\eta_{\rm SN}$
accounts for the fact that not all WBCSN remnant shocks are able to accelerate
particles up to the maximum energy given by Eq. 9, especially if their
progenitor explodes deep inside the cluster.

###  4.3 Propagation

Eq. 14 provides the spectrum of a given element of atomic number $Z$ released
in the Galaxy. Diffusive propagation in the galactic magnetic fields is then
taken into account using a one-zone model, which allows one to write the
spectral energy distribution of relativistic particles (per unit surface, per
unit solid angle, per unit time) as:

$\left.\frac{\mathrm{d}N}{\mathrm{d}E}\right|_{Z}=\frac{H_{\rm
disk}^{2}}{V_{\rm disk}}\frac{c}{4\pi D_{\rm Gal}(E,Z)}\phi_{\rm c}(E,Z)\,,$
(15)

where $H_{\rm disk}$ and $V_{\rm disk}$ are respectively the height and volume
of the Galactic disk, and $D_{\rm Gal}$ is the galactic diffusion coefficient
parametrized as follows:

$D_{\rm Gal}(E,Z)=D_{0}\left(\frac{E}{Z}\right)^{\delta}\,.$ (16)

The value of the galactic diffusion coefficient at 1 GV only affects the
global normalisation of the spectrum and can be adjusted together with the
acceleration efficiency. In the following we assume $D_{0}=10^{28}$ cm2/s. The
parametrization given by Eq. 16 is only valid up to the rigidity for which the
particle mean free path $D/c$ equals the scale height of the Galactic disk.
Assuming the latter to be 500 pc gives a break rigidity of $(5\times
10^{3})^{1/\delta}$ GV, which, for consistency, must be larger than the
maximum rigidity achieved around the fastest WBCSN remnant shocks. This
implies that $\delta$ must be smaller than 0.5 (smaller than 0.4 if we assume
a disk height of 100 pc). In the following, the best fit value will be found
to be $\delta=0.35$, therefore meeting this requirement.

Parameter | Notation | Value
---|---|---
Minimum initial velocity of a SN shock | $u_{\rm min}$ | 3000 km/s
Maximum initial velocity of a SN shock | $u_{\rm max}$ | 30000 km/s
Scaling of SN shocks velocity distribution | $\xi(u)$ | $\propto u^{-1}$
Fraction of massive stars which are clustered | $f_{\rm c}$ | 80%
Rate of CCSN | $f_{\rm SN}$ |  $2/$century
Average energy of a SN explosion | $E_{\rm SN}$ |  $10^{51}$ erg
∗Fraction of WBC in the MSC population | $f_{\rm WBC}$ | 16%
∗Probability to launch a fast shock in the collective WBC wind | $\eta_{\rm SN}$ | 0.5
Wind mechanical power per young massive star | $\bar{P}_{\rm w}$ |  $5\times 10^{35}$ erg/s
Number of massive stars in the Galaxy | $N_{*,\rm Gal}$ | 400k
Maximum energy achieved at WTS | $E_{\rm max,\rm WTS}$ | 1 PeV
∗Proton injection efficiency | $f_{1}$ | 0.18
Composition at 1 TeV | $f_{Z}$ | Hörandel (2003)
∗Spectral index deviation between LCSN and WBCSN spectra | $\Delta s$ | 0.17
CR injection momentum | $p_{0}$ | 10 MeV/c
Maximum energy of protons from a 5000 km/s LCSN shock | - | 0.1 PeV
Maximum energy of protons from a 5000 km/s WBCSN shock | - | 4 PeV
∗Galactic diffusion coefficient scaling | $\delta$ | 0.35
Galactic diffusion coefficient at 1 GV | $D_{0}$ |  $10^{28}$ cm2/s
Height of the Galactic disk | $H_{\rm disk}$ | 0.5 kpc
Volume of the Galactic disk | $V_{\rm disk}$ | 400 kpc3
Table 2: Exhaustive list of input parameters. Parameters marked with a star
are “fitting” parameters, which are allowed to vary within reasonable bounds
(see text).

###  4.4 All-particle spectrum

We compute the all particle spectrum by summing Eq. 15 up to $Z=31$. While the
normalisation of the proton spectrum is set by the acceleration efficiency
parameter, the contribution from each nucleus is then normalised with respect
to that of protons to ensure that we retrieve the elemental abundances at 1
TeV given in Hörandel (2003). The resulting spectrum depends on a number of
parameters, which are listed in Table 2. Highlighted in bold are parameters
which we allow to vary in order to fit the data. The other parameters are
fixed to the values indicated in the table.

Among the fitting parameters, the proton injection efficiency sets the overall
normalisation of the spectrum. The slope is affected by the galactic diffusion
coefficient scaling $\delta$, while the relative steepness between the LCSN
and WBCSN components (i.e. the shape of the “knee”) can be adjusted by varying
the spectral index deviation $\Delta s$. The fraction of WBC in the MSC
population sets the relative normalisation of the two dominant LCSN and WBCSN
components, while the parameter $\eta_{SN}$, which encodes the uncertainties
about the fate of SNRs expanding deep inside WBCs, impacts the relative
normalisation of the wind component. We vary these five parameters in order to
fit the all-particle spectrum of CRs. It should be noted that these parameters
are tightly constrained theoretically. For instance, the galactic diffusion
coefficient index is not expected to be smaller than $1/3$ (which corresponds
to a Kolmogorov turbulence regime). However, a higher diffusion coefficient
index leads to a steeper spectrum and therefore requires an increase of the
normalisation in order to match the measured flux. This can only be achieved
by raising the proton injection efficiency above 20%, which becomes higher
than what is suggested by theory. In the end, the five parameter values
displayed in Table 2 are the result of a compromise between the normalisation
of the components and the shape of the spectrum.

Figure 3: All-particle spectrum obtained from our computation, compared with
recent data, partly extracted from the compilation by C. Evoli
(https://github.com/carmeloevoli/The_CR_Spectrum): Auger (Verzi, 2019; Abreu
et al., 2021), HAWC (Alfaro et al., 2017), IceCube/IceTop (Aartsen et al.,
2019), KASCADE (Finger, 2011), NUCLEON (Grebenyuk et al., 2019), Telescope
Array (Matthews & Telescope Array Collaboration, 2017), Tibet (Amenomori et
al., 2008), Tunka (Budnev et al., 2020).

The all-particle spectrum obtained subsequently to this fitting procedure is
shown in Figure 3. It displays a good agreement with the measurements, given
the spread of the data. In particular, our calculation results in an excellent
match with the HAWC, Nucleon, Kascade, Tunka and Auger data, from 10 TeV to a
few hundreds of PeV. The galactic contribution at UHE, dominated by WBCSN (red
dashed line), succeeds in reproducing the data up to 200 PeV, and starts to
undershoot the data beyond 300 PeV. This demonstrates that supernovae
exploding within MSCs can account for the UHE GCR spectrum with reasonable
choice of the parameters. In particular, we show that a few fast SN shocks
(2/millenium) can alone account for the particles seen beyond 10 PeV, up to
hundreds of PeV.

As mentioned in the introduction, composition measurements point to a
transition from GCRs to extragalactic CRs taking place at a few hundreds of
PeV, with data suggesting the mean atomic mass peaks at $100$ PeV and then
decreases up to the ankle (e.g. Supanitsky, 2022). In order to probe how a
sub-ankle extragalactic component complements our scenario for the origin of
GCRs, we added the contribution modelled in Unger et al. (2015) (considering
the photo-disintegration of nuclei close to the extragalactic sources). To
this aim, we reproduce the “UFA” component calculated in Thoudam et al.
(2016), which includes propagation effects. It is shown by the grey dashed
line in Figure 3. This showcases that a sub-ankle extragalactic component can
indeed compensate the fall-off of the galactic contribution, leading to a very
good fit of the Auger data up to the ankle. Similar conclusions are obtained
using the calculations of Globus et al. (2015) or the model by Kotera &
Lemoine (2008).

Figure 4: All-particle spectrum obtained by varying the fraction of wind-
blowing clusters $f_{\rm WBC}$ with values indicated in the legend. The other
parameters are kept as listed in Table 2.

Any model trying to reproduce, with two components, a concave transition
between two featureless power-laws must require some level of tuning in order
to avoid unexpected bumps or drops (Parizot, 2014). We show in Figure 4 that
the fraction of WBCs, which sets the normalisation of the second, subdominant,
component in our model, must indeed be within 5% to 25% in order to reproduce
a knee shape. A fraction smaller than 5% produces a drop around 100 PeV, while
a fraction larger than 25% produces a bump around 1 PeV. Without prior
knowledge on the nature of this second component, it might seem to be a
coincidence to fall in this narrow range. However, the knowledge that this
component comes from a subset of the galactic MSCs already sets a strong
constraint on its normalisation factor, as it implies that it must be smaller
than 1 and, in fact, most likely smaller than 50% since we know that the
majority of clusters are not blowing a collective wind. More detailed
investigation on the star cluster population in the Milky Way, as done in
Section 2, show that the fraction of WBCs falls exactly in the required range
(about 16%, keeping in mind the uncertainties discussed in Section 2) and no
posterior adjustment on this parameter is needed.

Our model relies on the assumption that the distribution of SN shock
velocities can be described by an extended tail. With a distribution of SN
shock velocities scaling as $1/u$, one expects two SN per millennium launching
a shock with velocity 30 000 km/s. This is obviously required since the
maximum proton energy for a standard 5000 km/s shock expanding close to the
edge of a WBC is 4 PeV, leading to a maximum iron energy of 100 PeV, which
would not be enough to reproduce the all-particle spectrum up to the ankle,
even considering a sub-ankle extragalactic component. We show in Figure 5 the
results which would have been obtained assuming different scalings of the
distribution of SN shock initial velocities $\xi(u)$. Only a distribution
$\xi(u)=1/u$ is able to reproduce the Auger data up to 1 EeV. A $1/u^{2}$
scaling fails to match the systematics errors and, as expected, an
exponentially decreasing distribution falls short of the data beyond 100 PeV.
This shows that the $1/u$ distribution of SNR shock velocities is a critical
ingredient of our model, together with the maximum proton energy achieved
around standard SNR shocks expanding close to WBCs.

In contrast, the analysis done by Sveshnikova (2003) suggests a steeper
distribution tail, scaling as $\sim u^{-3}$. This is however inferred from an
energy distribution tail extending up to about $10^{53}$ erg (“hypernovae”)
under the hypothesis of constant ejected mass. The latter is questionable, and
one can find a hint in recent numerical simulations that the ejected mass is
anti-correlated with the explosion energy (Ebinger et al., 2020). For
instance, SN1987A had an estimated explosion energy around $10^{51}$ erg but
launched a shock of velocity beyond 30 000 km/s (Blinnikov et al., 2000).
Nevertheless, we can incorporate the “Sveshnikova” distribution in our model,
assuming $\xi(u)\propto u^{-3}$, and, for simplicity, $E_{\rm
SN}=u_{5}^{2}\times 10^{51}$ erg. The result is shown by the light green curve
in Figure 5. Due to the assumption of constant ejected mass, the fastest SN
shocks are also the most energetic. This makes their contribution to the CR
spectrum larger, which compensates for their scarcity. Incidentally, to obtain
the result shown in Figure 5 in this specific case, it was necessary to reduce
the proton injection efficiency to 14% in order to avoid overshooting the
data. This is the only parameter which differs compared to the other curves.

In the end, it is unclear which model of SN distribution is the closest to
reality, which is why, in the present paper, we agnostically chose to stick to
a constant explosion energy of $10^{51}$ erg and to set fiducial velocity
distributions.

Figure 5: All-particle spectrum obtained from various assumptions on the
distribution of SN shock velocities $\xi$: three different power-laws, an
exponential cut-off at $u=3000$ km/s, and a model inspired by the work of
Sveshnikova (2003) where the explosion energy follows as well a power-law
distribution and the normalisation of the CR spectra is set accordingly (see
text). All parameters are kept as listed in Table 2, except for the latter
case where $f_{1}=0.14$.

###  4.5 Proton spectrum

Having computed the all-particle spectrum, we can extract the spectrum of the
protons and compare it with available data. The result is shown in Figure 6.
Interpreting the match between our model and the data is less straightforward
in the case of the protons. Indeed, the discrepancy between experiments,
mostly due to the high uncertainties in the modelling of hadronic air showers,
does not set tight constraints beyond 1 PeV. We nevertheless find good
agreement below 1 PeV, where our model meets most of the data within error
bars. It does seem that the WBCSN contribution tends to overshoot the Kascade
and Tibet measurements above 1 PeV, however this is actually expected, as
WBCSN accelerate mostly material blown by WR stars, whose composition contains
less hydrogen than the average superbubble composition around loose clusters
(about 80%, as suggested by Figure 3 in Tatischeff et al., 2021).

Figure 6: Proton spectrum obtained from our computation, compared with recent
data (Panov et al., 2009; Adriani et al., 2019; Yoon et al., 2017; An et al.,
2019; Aartsen et al., 2019; Antoni et al., 2005; Apel et al., 2013; Grebenyuk
et al., 2019; Amenomori et al., 2006; Epimakhov et al., 2013).

###  4.6 Composition

Finally, we compute the mean logarithmic mass as inferred from our model and
compare it with available data. The result is shown in Figure 7. Once again,
the variance of the data from experiment to experiment is large, due to the
uncertainty in the modelling of the air shower, in particular when computing
the average depth of the shower maximum of a given species. Most recent data
(e.g. from the TALE or Tunka observatories) suggest that the composition
around 100 PeV is lighter than previously thought (e.g. compared to older
measurements by Kascade). Nevertheless, our model matches reasonably the trend
of the data, and in particular the peak around 100 PeV which is due, in our
scenario, to the transition between the heaviest GCRs which are the only ones
which WBCSN shocks can accelerate beyond 100 PeV, and the light extragalactic
component due to, e.g. the photo-disintegration of nuclei around sources
(Unger et al., 2015; Globus et al., 2015).

Figure 7: CR composition (mean logarithmic atomic mass) obtained from our
computation, compared with recent data taken from Sveshnikova et al. (2013);
Supanitsky (2022), and references therein.

##  5 Conclusions

Only fast SNR shocks expanding in the collective outflow around a wind-blowing
cluster might, under favourable conditions, accelerate protons beyond 10 PeV,
which is required to account for the galactic component of the local CR
spectrum. In this work, we have estimated the rate of these SNR shocks by
analysing the massive star cluster population of the Milky Way. We found that
about 15% of the MSCs blow a collective wind. Taking into account a clustering
fraction of massive stars of 80% and a probability of 0.5 that a SN explodes
close enough to the cluster edge in order to launch a fast shock in the
collective wind, the rate of these events is less than 6% that of galactic
CCSN, i.e. about 1 per millennium for 3000 km/s shocks, and less than 1 per 10
kyr for 30 000 km/s shocks. We nevertheless demonstrated that this is enough
to account for the spectrum of galactic CRs from TeV bands up to hundreds of
PeV, where a light sub-ankle extragalactic component is expected to take over.
The knee of the CR spectrum is then found to be the signature of the
transition between the contribution of standard SNR shocks expanding in loose
clusters (too evolved or too extended to blow a collective wind) and that of
scarce SNR shocks expanding in the vicinity of wind-blowing clusters (young
and compact enough to sustain a collective wind). The continuity at the knee
is found to be a consequence of the second population being naturally a
fraction of order $10-20\%$ the galactic CCSN population. A good agreement is
obtained as well with the spectrum of the protons and data on composition.

Our model thus provides a possible explanation for GCRs, as produced in star
clusters without ad-hoc hypotheses, and with a sensible set of free-
parameters. Crucially, all parameter values are physically motivated and, as
summarised in Table 2, within expected bounds. The main uncertainty of the
model is the distribution of initial velocities of SN shocks. As fast SN
shocks are needed to accelerate heavy elements up to hundreds of PeV, it is
necessary that such a distribution displays a rather extended tail. Such a
hypothesis is not excluded by available data, although the interpretation of
the latter is not straightforward. Different assumptions can be set for
simplicity, e.g. constant explosion energy, or constant ejected mass, which
all succeed to account for the CR spectrum. Finally, it yet remains to be
proven that SNR shocks expanding in cluster winds are able to accelerate
protons up to several PeV, which in particular requires a high level of
turbulence in the cluster, which is at the moment mostly supported by
theoretical arguments.

## Acknowledgements

TV acknowledges L. Harer, S. Gabici, L.-M. Bourguinat and M. Tsirou for
helpful discussions.

## Data Availability

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

## References

  * Aartsen et al. (2019) Aartsen M. G., et al., 2019, Phys. Rev. D, 100, 082002
  * Abeysekara et al. (2021) Abeysekara A. U., et al., 2021, Nature Astronomy, 5, 465
  * Abreu et al. (2021) Abreu P., et al., 2021, European Physical Journal C, 81, 966
  * Adriani et al. (2019) Adriani O., et al., 2019, Phys. Rev. Lett., 122, 181102
  * Aharonian et al. (2019) Aharonian F., Yang R., de Oña Wilhelmi E., 2019, Nature Astronomy, 3, 561
  * Aharonian et al. (2022) Aharonian F., et al., 2022, arXiv e-prints, p. arXiv:2207.10921
  * Alfaro et al. (2017) Alfaro R., et al., 2017, Phys. Rev. D, 96, 122001
  * Amenomori et al. (2006) Amenomori M., et al., 2006, Advances in Space Research, 37, 1938
  * Amenomori et al. (2008) Amenomori M., et al., 2008, ApJ, 678, 1165
  * An et al. (2019) An Q., et al., 2019, Science Advances, 5, eaax3793
  * Antoni et al. (2005) Antoni T., et al., 2005, Astroparticle Physics, 24, 1
  * Apel et al. (2013) Apel W. D., et al., 2013, Astroparticle Physics, 47, 54
  * Axford (1994) Axford W. I., 1994, ApJS, 90, 937
  * Badmaev et al. (2022) Badmaev D. V., Bykov A. M., Kalyashova M. E., 2022, MNRAS,
  * Bell & Lucek (2001) Bell A. R., Lucek S. G., 2001, MNRAS, 321, 433
  * Bell et al. (2011) Bell A. R., Schure K. M., Reville B., 2011, MNRAS, 418, 1208
  * Bell et al. (2013) Bell A. R., Schure K. M., Reville B., Giacinti G., 2013, MNRAS, 431, 415
  * Bell et al. (2019) Bell A. R., Matthews J. H., Blundell K. M., 2019, MNRAS, 488, 2466
  * Berezinsky et al. (2006) Berezinsky V., Gazizov A., Grigorieva S., 2006, Phys. Rev. D, 74, 043005
  * Biermann & Cassinelli (1993) Biermann P. L., Cassinelli J. P., 1993, A&A, 277, 691
  * Blinnikov et al. (2000) Blinnikov S., Lundqvist P., Bartunov O., Nomoto K., Iwamoto K., 2000, ApJ, 532, 1132
  * Bosch-Ramon et al. (2005) Bosch-Ramon V., Aharonian F. A., Paredes J. M., 2005, A&A, 432, 609
  * Budnev et al. (2020) Budnev N. M., et al., 2020, Astroparticle Physics, 117, 102406
  * Bykov (2001) Bykov A. M., 2001, Space Sci. Rev., 99, 317
  * Bykov & Toptygin (2001) Bykov A. M., Toptygin I. N., 2001, Astronomy Letters, 27, 625
  * 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
  * Cantat-Gaudin et al. (2020) Cantat-Gaudin T., et al., 2020, A&A, 640, A1
  * Cao et al. (2021) Cao Z., et al., 2021, Nature, 594, 33
  * Caprioli (2011) Caprioli D., 2011, J. Cosmology Astropart. Phys., 2011, 026
  * Caprioli (2012) Caprioli D., 2012, J. Cosmology Astropart. Phys., 2012, 038
  * Casse & Paul (1980) Casse M., Paul J. A., 1980, ApJ, 237, 236
  * Celli et al. (2019) Celli S., Morlino G., Gabici S., Aharonian F. A., 2019, MNRAS, 490, 4317
  * Cesarsky & Montmerle (1983) Cesarsky C. J., Montmerle T., 1983, Space Sci. Rev., 36, 173
  * Cheng et al. (2012) Cheng K. S., Chernyshov D. O., Dogiel V. A., Ko C. M., Ip W. H., Wang Y., 2012, ApJ, 746, 116
  * Crowther (2007) Crowther P. A., 2007, ARA&A, 45, 177
  * Das et al. (2022) Das S., Brose R., Meyer D. M. A., Pohl M., Sushch I., Plotko P., 2022, A&A, 661, A128
  * 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
  * Epimakhov et al. (2013) Epimakhov S., et al., 2013, in International Cosmic Ray Conference. p. 818
  * Ferrand & Marcowith (2010) Ferrand G., Marcowith A., 2010, A&A, 510, A101
  * Finger (2011) Finger M. R., 2011, PhD thesis, KIT, Karlsruhe
  * Finke & Dermer (2012) Finke J. D., Dermer C. D., 2012, ApJ, 751, 65
  * Gaggero et al. (2018) Gaggero D., Zandanel F., Cristofari P., Gabici S., 2018, MNRAS, 475, 5237
  * Giacinti et al. (2012) Giacinti G., Kachelrieß M., Semikoz D. V., Sigl G., 2012, J. Cosmology Astropart. Phys., 2012, 031
  * Globus et al. (2015) Globus N., Allard D., Parizot E., 2015, Phys. Rev. D, 92, 021302
  * Grebenyuk et al. (2019) Grebenyuk V., et al., 2019, Advances in Space Research, 64, 2546
  * Gupta et al. (2020) Gupta S., Nath B. B., Sharma P., Eichler D., 2020, MNRAS, 493, 3159
  * HESS Collaboration et al. (2016) HESS Collaboration et al., 2016, Nature, 531, 476
  * Hamann et al. (2019) Hamann W. R., et al., 2019, A&A, 625, A57
  * Higdon & Lingenfelter (2005) Higdon J. C., Lingenfelter R. E., 2005, ApJ, 628, 738
  * Hillas (2005) Hillas A. M., 2005, Journal of Physics G Nuclear Physics, 31, R95
  * Hörandel (2003) Hörandel J. R., 2003, Astroparticle Physics, 19, 193
  * Jokipii & Morfill (1987) Jokipii J. R., Morfill G., 1987, ApJ, 312, 170
  * Kachelrieß & Semikoz (2019) Kachelrieß M., Semikoz D. V., 2019, Progress in Particle and Nuclear Physics, 109, 103710
  * Kachelrieß et al. (2017) Kachelrieß M., Kalashev O., Ostapchenko S., Semikoz D. V., 2017, Phys. Rev. D, 96, 083006
  * Kirk et al. (1996) Kirk J. G., Duffy P., Gallant Y. A., 1996, A&A, 314, 1010
  * Klepach et al. (2000) Klepach E. G., Ptuskin V. S., Zirakashvili V. N., 2000, Astroparticle Physics, 13, 161
  * Kotera & Lemoine (2008) Kotera K., Lemoine M., 2008, Phys. Rev. D, 77, 023005
  * Krumholz et al. (2019) Krumholz M. R., McKee C. F., Bland-Hawthorn J., 2019, ARA&A, 57, 227
  * Lagage & Cesarsky (1983) Lagage P. O., Cesarsky C. J., 1983, A&A, 125, 249
  * Lamb et al. (2010) Lamb J. B., Oey M. S., Werk J. K., Ingleby L. D., 2010, ApJ, 725, 1886
  * Limongi & Chieffi (2006) Limongi M., Chieffi A., 2006, ApJ, 647, 483
  * Lingenfelter (2018) Lingenfelter R. E., 2018, Advances in Space Research, 62, 2750
  * Matthews & Telescope Array Collaboration (2017) Matthews J., Telescope Array Collaboration 2017, in 35th International Cosmic Ray Conference (ICRC2017). p. 1096
  * McCray & Kafatos (1987) McCray R., Kafatos M., 1987, ApJ, 317, 190
  * Merten et al. (2018) Merten L., Bustard C., Zweibel E. G., Becker Tjus J., 2018, ApJ, 859, 63
  * Miller & Scalo (1979) Miller G. E., Scalo J. M., 1979, ApJS, 41, 513
  * Molinari et al. (2014) Molinari S., et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 125
  * Montmerle (1979) Montmerle T., 1979, ApJ, 231, 95
  * Morlino et al. (2021) Morlino G., Blasi P., Peretti E., Cristofari P., 2021, MNRAS, 504, 6096
  * Nugis & Lamers (2000) Nugis T., Lamers H. J. G. L. M., 2000, A&A, 360, 227
  * Panov et al. (2009) Panov A. D., et al., 2009, Bulletin of the Russian Academy of Sciences, Physics, 73, 564
  * Parizot (2014) Parizot E., 2014, Nuclear Physics B - Proceedings Supplements, 256-257, 197–212
  * Parizot et al. (2004) Parizot E., Marcowith A., van der Swaluw E., Bykov A. M., Tatischeff V., 2004, A&A, 424, 747
  * Poelarends et al. (2008) Poelarends A. J. T., Herwig F., Langer N., Heger A., 2008, ApJ, 675, 614
  * Ptuskin et al. (2010) Ptuskin V., Zirakashvili V., Seo E.-S., 2010, ApJ, 718, 31
  * Sander et al. (2012) Sander A., Hamann W. R., Todt H., 2012, A&A, 540, A144
  * Sander et al. (2019) Sander A. A. C., Hamann W. R., Todt H., Hainich R., Shenar T., Ramachandran V., Oskinova L. M., 2019, A&A, 621, A92
  * Seo et al. (2018) Seo J., Kang H., Ryu D., 2018, Journal of Korean Astronomical Society, 51, 37
  * Smartt et al. (2009) Smartt S. J., Eldridge J. J., Crockett R. M., Maund J. R., 2009, MNRAS, 395, 1409
  * Sukhbold et al. (2016) Sukhbold T., Ertl T., Woosley S. E., Brown J. M., Janka H. T., 2016, ApJ, 821, 38
  * Supanitsky (2022) Supanitsky A. D., 2022, Galaxies, 10, 75
  * Sveshnikova (2003) Sveshnikova L. G., 2003, A&A, 409, 799
  * Sveshnikova et al. (2013) Sveshnikova L. G., et al., 2013, in International Cosmic Ray Conference. p. 1538
  * Takamoto & Kirk (2015) Takamoto M., Kirk J. G., 2015, ApJ, 809, 29
  * Tatischeff et al. (2021) Tatischeff V., Raymond J. C., Duprat J., Gabici S., Recchia S., 2021, MNRAS, 508, 1321
  * Thoudam et al. (2016) Thoudam S., Rachen J. P., van Vliet A., Achterberg A., Buitink S., Falcke H., Hörandel J. R., 2016, A&A, 595, A33
  * Unger et al. (2015) Unger M., Farrar G. R., Anchordoqui L. A., 2015, Phys. Rev. D, 92, 123001
  * Verzi (2019) Verzi V., 2019, in 36th International Cosmic Ray Conference (ICRC2019). p. 450
  * Vieu et al. (2020) Vieu T., Gabici S., Tatischeff V., 2020, MNRAS, 494, 3166
  * Vieu et al. (2022a) Vieu T., Gabici S., Tatischeff V., 2022a, MNRAS, 510, 2529
  * Vieu et al. (2022b) Vieu T., Gabici S., Tatischeff V., Ravikularaman S., 2022b, MNRAS, 512, 1275
  * Vieu et al. (2022c) Vieu T., Reville B., Aharonian F., 2022c, MNRAS, 515, 2256
  * Voelk & Biermann (1988) Voelk H. J., Biermann P. L., 1988, ApJ, 333, L65
  * Weaver et al. (1977) Weaver R., McCray R., Castor J., Shapiro P., Moore R., 1977, ApJ, 218, 377
  * Xu & Lazarian (2022) Xu S., Lazarian A., 2022, ApJ, 925, 48
  * Yang et al. (2018) Yang R.-z., de Oña Wilhelmi E., Aharonian F., 2018, A&A, 611, A77
  * Yoon et al. (2017) Yoon Y. S., et al., 2017, ApJ, 839, 5
  * Zirakashvili & Völk (2006) Zirakashvili V. N., Völk H. J., 2006, Advances in Space Research, 37, 1923