# Between the cosmic-ray ‘knee’ and the ‘ankle’: Contribution from star
clusters

Sourav Bhadra Raman Research Institute, Sadashiva Nagar, Bangalore 560080,
India  Joint Astronomy Programme, Department of Physics, Indian Institute of
Science, Bangalore 560012, India  Satyendra Thoudam  Department of Physics,
Khalifa University, PO Box 127788, Abu Dhabi, United Arab Emirates.  Biman B
Nath  Raman Research Institute, Sadashiva Nagar, Bangalore 560080, India
Prateek Sharma  Joint Astronomy Programme, Department of Physics, Indian
Institute of Science, Bangalore 560012, India

###### Abstract

We show that massive young star clusters may be possible candidates that can
accelerate Galactic cosmic rays (CRs) in the range of $10^{7}\hbox{--}10^{9}$
GeV (between the ‘knee’ and ‘ankle’). Various plausible scenarios such as
acceleration at the wind termination shock (WTS), supernova shocks inside
these young star clusters, etc. have been proposed, since it is difficult to
accelerate particles up to the $10^{7}\hbox{--}10^{9}$ GeV range in the
standard paradigm of CR acceleration in supernova remnants. We consider a
model for the production of different nuclei in CRs from massive stellar winds
using the observed distribution of young star clusters in the Galactic plane.
We present a detailed calculation of CR transport in the Galaxy, taking into
account the effect of diffusion, interaction losses during propagation, and
particle re-acceleration by old supernova remnants to determine the all-
particle CR spectrum. Using the maximum energy estimate from the Hillas
criterion, we argue that a young massive star cluster can accelerate protons
up to a few tens of PeV. Upon comparison with the observed data, our model
requires a CR source spectrum with an exponential cutoff of $5\times 10^{7}Z$
GeV ($50\,Z$ PeV) from these clusters together with a cosmic-ray injection
fraction of $\sim 5\%$ of the wind kinetic energy. We discuss the possibility
of achieving these requirements in star clusters, and the associated
uncertainties, in the context of considering star clusters as the natural
accelerator of the ‘second component’ of Galactic cosmic rays.

(ISM:) cosmic rays, galaxies: star clusters: general, ISM: bubbles, stars:
winds, outflows, acceleration of particles, shock waves

##  1 Introduction

Cosmic rays (CRs hereafter) are high-energy particles that span an extensive
range of energy from $1$ GeV to $\sim 10^{11}$ GeV. Lower energy CRs up to
$\sim 10^{5-6}$ GeV are believed to be accelerated by supernova shocks (Lagage
& Cesarsky 1983; Axford 1994). This dominant acceleration mechanism, revealed
by both theoretical (Fermi 1949; Axford et al. 1977; Bell 1978; Blandford &
Ostriker 1978; Blasi 2013; Caprioli 2015) and observational (Drury et al.,
1994; Ackermann et al., 2013; H. E. S. S. Collaboration et al., 2018) studies,
is diffusive shock acceleration (DSA), a first-order Fermi acceleration
process in which $\sim 10$% of the shock energy is expected to be converted to
cosmic rays. Although the acceleration mechanism continues to work throughout
the active stage of a supernova remnant (SNR) until it becomes
indistinguishable from the ambient interstellar medium after $\sim
10^{5}-10^{6}$ years, most of the particle acceleration occurs during the un-
decelerated blast wave phase, which lasts for $\leq 10^{3}$ years (Lagage &
Cesarsky, 1983). This limits the maximum CR energy that can be accelerated in
SNRs because the acceleration time of CRs cannot be longer than the age of the
SNR (Morlino, 2017). Considering nonlinear effects such as the scattering of
the cosmic rays by the waves they generate themselves and assuming the
magnetic flux density of the interstellar magnetic field to be $\sim\mu$G,
Lagage & Cesarsky (1983) estimated the upper limit of CR energy in supernova
remnants to be $\sim 10^{5}$ GeV per nucleon.

Preliminary observations seem to align with these theoretical concepts. Suzuki
et al. (2022) reported cutoff energy of around TeV from $\gamma$-ray
observations of $15$ supernova remnants. Recently, LHAASO (Large High Altitude
Air Shower Observatory), and Tibet air shower observations have identified a
number of PeVatron candidates (Cao et al. 2021; Tibet AS$\gamma$ Collaboration
et al. 2021), which may include a few SNRs. These theoretical and
observational developments suggest cutoff energy in the range $10^{5}-10^{6}$
GeV for SNRs. At the high-energy end, cosmic rays above $\sim 10^{9}$ GeV are
considered to have an extragalactic origin, possibly originating from galaxy
clusters (Kang et al., 1996), radio galaxies (Rachen & Biermann, 1993), AGN
jets (Mannheim et al., 2000) or gamma-ray bursts (Waxman, 1995).

There is a gap between the contribution from SNRs and the extragalactic
component, which lies in the range of $\sim 10^{7}-10^{9}$ GeV, the region
between the so-called ‘knee’ and ‘ankle’ (also known as the ‘shin’ region). To
explain this gap in the all-particle CR spectrum, a few models have been
proposed in the literature. Biermann & Cassinelli (1993); Thoudam et al.
(2016) have discussed the explosion of Wolf-Rayet stars embedded in the wind
material from the same stars as a potential acceleration site of CRs in the
range of $\sim 10^{7}-10^{9}$ GeV. However, there may be some problems with
this scenario. A uniform distribution of Wolf-Rayet stars in the Galaxy was
assumed, which is unrealistic. Moreover, there are many uncertainties in the
crucial parameter of the magnetic field of the Wolf-Rayet stars. For a proton
cutoff energy of $1.1\times 10^{8}$ GeV, the surface magnetic field of a Wolf-
Rayet star is required to be $\approx 10^{4}$ G in this model (Thoudam et al.,
2016), while realistic predictions for the same are in the range of a few
hundred G (Neiner et al., 2015; Blazère et al., 2015). Although no direct
magnetic signature has been detected in any of the Wolf-Rayet stars, using
Bayesian statistics, Bagnulo et al. (2020) have estimated their surface
magnetic fields to be of the order of a few kiloGauss.

The idea of Galactic wind termination shock to accelerate high-energy CRs also
has problems. The effect of Galactic winds on the transport of cosmic rays in
the Galaxy has been discussed in detail (Lerche & Schlickeiser 1982; Bloemen
et al. 1993; Strong & Moskalenko 1998; Jones et al. 2001; Breitschwerdt et al.
2002). Following these developments, Jokipii & Morfill (1987); Zirakashvili &
Völk (2006); Thoudam et al. (2016) introduced these CRs originating from
Galactic wind termination shock as the possible candidate for the ‘second’
(between ‘knee’ and ‘ankle’) component of Galactic cosmic rays. The cosmic
rays originating from the Galactic wind (GW-CRs) are believed to mostly
contribute to the higher energy range. This is due to the increasing effect of
advection over diffusion at lower energy, preventing particles from reaching
the Galactic disk. Higher energy particles, which diffuse relatively faster,
can overcome the advection and reach the disk more effectively. Thoudam et al.
(2016) have used a distance of $\sim 100$ kpc for the Galactic wind
termination shock. Bustard et al. (2017) have argued that in order for the CRs
to reach $10^{8}$ GeV, either the outflow speed needs to be of order $\sim
1000$ km s-1 or the magnetic field needs to be amplified. However, realistic
simulations of outflows from Milky Way-type galaxies do not find signatures of
such strong outflows/shocks. Sarkar et al. (2015) showed that the outer shock
due to the Galactic wind weakens and continues to propagate as a sound wave
through the circum-galactic medium. The termination shock remains confined
within $\lesssim 10$s of kpc and disappears after the mechanical power is
stopped being injected. Also, the observed nuclear abundances suggest lighter
nuclei in contrast to the expectation from the Galactic wind model in the
$10^{7}-10^{9}$ GeV energy range. Thus, this model has been disfavoured. In
order to explain the observed all-particle spectrum in the range
$10^{7}-10^{9}$ GeV, an appropriate model of CRs is still required.

Coming back to the DSA mechanism of CR acceleration in supernova shocks, this
standard scenario is known to bear several ailing problems (e.g., Gabici et
al. 2019). The acceleration scenario cannot explain some of the observed
features of cosmic rays like excess of Ne22/Ne20 in CRs compared to standard
cosmic abundances in ISM (Binns et al. 2008; Wiedenbeck et al. 1999), proton
acceleration up to greater than PeV ($10^{6}$ GeV) energy range, and so on.
Various additional CR acceleration sites are reported in the literature to
solve these paradigms; young massive star clusters are one of those other
possible sources of Cosmic rays in our Galaxy (Bykov, 2014; Knödlseder, 2013;
Aharonian et al., 2019). Recently, the $\gamma$-ray observations by LHAASO,
HESS, Fermi-LAT, and HAWC have provided evidence of CR acceleration up to very
high energy in a few massive star clusters like Westerlund1 and Cygnus
(Aharonian et al., 2019; Abeysekara et al., 2021). These star-forming regions
have been discussed as potential candidates for CR accelerators (e.g. Bykov
2014 ); these $\gamma$-ray observations have strengthened the hypothesis of CR
acceleration in these environments. Recently, Gupta et al. (2020) has shown
that the excess (22Ne/20Ne) ratio can be explained by considering wind
termination shock (WTS) of massive star clusters as CR accelerators. Recently,
Tatischeff et al. (2021) showed that the refractory elements of Galactic
cosmic rays are produced in super-bubbles. This theoretical and observational
evidence prod us to study the total contribution of Cosmic rays originating
from the distribution of massive star clusters in our Galaxy.

Star clusters, which are the birthplace of massive stars (that ultimately
explode as SNe), give rise to continuous mass outflow in the form of stellar
wind. These are mainly located in dense molecular clouds and weigh of the
order of several thousand solar masses (Longmore et al., 2014). Star clusters
host massive stars as well as supernova explosions, which produce a low-
density bubble around them (Weaver et al., 1977; Gupta et al., 2018). Young
star clusters contain sufficient kinetic energy supplied by interacting
stellar winds, which can accelerate protons up to $\sim 10^{7}$ GeV.
Considering heavier nuclei, this cosmic ray component originating from star
clusters can, therefore, be considered as the second component of Galactic
cosmic rays, which can explain the observed all-particle spectrum in the
energy range of $10^{7}-10^{9}$ GeV. Bhadra et al. (2022), using hydrodynamic
simulation, showed that the observed distribution of $\gamma$-rays can be
explained by invoking cosmic ray acceleration in the Westerlund1 cluster.

Following the above discussion, it is clear that: (1) Galactic supernovae can
accelerate particles up to a few times $10^{6}$ GeV energy, and (2)
extragalactic components can explain the all-particle spectrum above $\sim
10^{9}$ GeV energy. The gap in the energy range cannot be explained using only
these two components, and we require another Galactic component to explain the
observed data in the range $10^{7}-10^{9}$ GeV. Our main focus in this paper
is the second component of Galactic cosmic rays. In this regard, we propose CR
contribution from the population of massive star clusters as a source of the
observed all-particle CR spectrum in the range $\sim 10^{7}-10^{9}$ GeV. This
may act as a bridge between the SNR component and the extragalactic component
and fill the gap in the desired energy range.

We begin with some basics in Section 2. The details of our proposed model are
described in Section 3. In Section 4, we present our results, and in Section 5
we compare our model with other models. In Section 6, we consider various
models for the extragalactic CR component. This is followed by further
discussion in Section 7 and a conclusion in Section 8.

##  2 Existing components of Cosmic rays

###  2.1 First Galactic component: SNR-CRs

As mentioned in the Introduction, supernova remnants are the most likely
candidate for cosmic-ray acceleration up to $\sim 10^{6}$ GeV energy (Lagage &
Cesarsky, 1983). The diffusive shock acceleration at strong shocks produces a
power-law spectrum with an index of $\sim-2$ (Krymskii 1977; Bell 1978;
Blandford & Ostriker 1978; Caprioli et al. 2011). We have adopted the model of
Thoudam et al. (2016) for the CR component from Galactic supernova remnants
(SNR-CR component). After the acceleration in the strong shock of supernova
remnants, CR particles escape the remnants and propagate through the
interstellar medium via diffusion. These CR particles can be re-accelerated
repetitively by expanding supernova remnant shock waves already existing in
the interstellar medium during their propagation. These shocks are mainly
produced by older remnants and are relatively weak.

We use the same contribution of the SNR-CR component as presented in Thoudam
et al. (2016). Their calculation assumes an exponential cutoff for the proton
source spectrum at $E_{c}=2.5\times 10^{6}$ GeV. This value has been chosen by
fitting the observed all-particle spectrum. The maximum energy of SNR-CRs
corresponds to the cutoff energy of iron nuclei, which is $26\times
E_{c}=6.5\times 10^{7}$ GeV. This result shows that SNR-CRs can contribute
only $\sim 30\%$ of the total observed intensity above $\sim 2\times 10^{7}$
GeV (Thoudam et al., 2016). Therefore, additional components are required to
explain the all-particle spectrum in the $\gtrsim 10^{7}$ GeV range.

###  2.2 Extragalactic component

Various previous works have already pointed out that the ‘ankle’-like feature
of the CR spectrum at $\gtrsim 10^{9}$ GeV can be explained if we consider the
propagation effects of the extragalactic component (mainly proton) in the
evolving microwave background (Hillas, 1967; Berezinskii & Grigor’eva, 1988;
Berezinsky et al., 2006; Aloisio et al., 2012, 2014). We consider two
different models for the extragalactic component: the ‘UFA model’ (Unger et
al., 2015) and a combination of minimal (di Matteo 2015), and PCS model
(Rachen 2015; Thoudam et al. 2016). We refer to this combined model as the
‘MPCS’ model.

Unger et al. (2015) considers acceleration of energetic nuclei at the shocks
associated with gamma-ray bursts or tidal disruption events, and photo-
disintegration of these particles in the photon background present inside the
source region. In this model, only the highest energy particles having an
escape time shorter than the photo-disintegration time can escape the source
region leading to a strong proton component in the energy region below the
ankle. We call this the ‘UFA’ model of the extragalactic component. In
addition to the all-particle CR spectrum, data of the primary composition in
the ultra-high energy range have become available in the last few years.

The ‘minimal model’ has been derived from CR composition measured by the
Pierre Auger Observatory (di Matteo 2015) and assumes uniformly distributed
sources in a comoving volume that produce a power-law cosmic ray spectrum with
some cutoff at a particular rigidity $R_{c}$ (rigidity is defined as $Apc/Ze$,
where $A/Z$ is nuclear mass/charge and $p$ is momentum, $e$ the charge of the
electron, and $c$ is the speed of light in vacuum). Above $\sim 3\times
10^{10}$ GeV, the spectrum exhibits a steep cut-off that is mainly due to the
intrinsic cut-off in the injection spectrum (di Matteo, 2015), and not due to
the GZK absorption (Greisen, 1966; Zatsepin & Kuz’min, 1966) during the
propagation.

The PCS (primordial cluster shock) model is based on the universal scaling
argument. It takes into account the acceleration of primordial proton and
helium mixture by primordial cluster shocks, which are mainly the accretion
shocks expected from clusters of galaxies during the structure formation. In
this scenario, the acceleration of CR particles is limited by losses due to
pair production in the CMB. This component is not expected to reach ultra-high
energies. Consequently, the minimal model plus the “primordial cluster
component” was introduced by Rachen (2015), where the acceleration of heavy
nuclei at shocks of gamma-ray bursts or in tidal disruption events are
considered.

Figure 1: Schematic diagram of a stellar wind bubble. The position of
termination shock is $R_{s}$; $R_{\rm cd}$ and $R_{\rm fs}$ are contact
discontinuity and forward shock positions, respectively.

##  3 Second Galactic component: cosmic rays from star clusters

The all-particle cosmic ray spectrum has two main features: a steepening of
the spectral index from $-2.7$ to $-3.1$ at about $3$ PeV, commonly known as
the ‘knee’, and a flattening back to $-2.7$ at about $4\times 10^{9}$ GeV,
generally known as the ‘ankle’. Therefore, we need to assume a cut-off in the
Galactic component immediately below the ‘ankle’ to explain the observed
spectrum. This is a ‘second knee’ feature in the CR spectrum. For a typical
magnetic field of $3$ $\mu$G in the Galaxy, cosmic rays with energy Z $\times
10^{8}$ GeV have a Larmor radius of $36/Z$ pc, which is much smaller than the
extent of the diffusion halo of the Galaxy. This implies that cosmic rays with
the energy around the second knee remain confined within the Galaxy. This also
suggests the observed cut-off at this energy is due to some CR accelerators
different from SNRs, as the latter accelerate particles only up to a few
$10^{6}$ GeV.

In the following, we discuss one potential scenario of another Galactic
component of CRs: the acceleration of cosmic rays by the young massive star
clusters, which we briefly mentioned in Section 1. It has especially been
speculated that the winds of massive stars may be a suitable location for the
acceleration of CRs (Cesarsky & Montmerle 1983; Webb et al. 1985; Gupta et al.
2018; Bykov et al. 2020). CRs can be accelerated in the fast stellar wind of
star clusters, and in particular, two scenarios can be important. Firstly, CR
acceleration in the wind termination shock (WTS) (Gupta et al., 2018), and
secondly, acceleration by supernova shocks embedded in the stellar winds (Vieu
et al., 2022). Those cosmic rays accelerated in young star clusters with age
$\leq 10$ Myr can contribute significantly to the observed total flux of CRs
(Gupta et al., 2020). Recently LHAASO has observed $\gamma$-rays in the PeV
energy range from young massive star clusters (Cao et al., 2021), which can be
associated with cosmic ray acceleration in those clusters.

Figure 1 shows a schematic diagram of a stellar wind bubble around a compact
star cluster. There are several distinct regions inside the bubble, such as
(a) the free wind region, where the stellar wind originating from the source
expands adiabatically, (b) the wind termination shock (WTS), (c) the shocked
wind region containing slightly denser gas, and (d) the outermost dense shell
containing the swept-up ambient gas. Cosmic rays can be accelerated in the
central region as well as in the shocks of the cluster. After getting
accelerated to very high energy, cosmic rays will diffuse outward from the
source into the ISM.

###  3.1 Distribution of star clusters in Galactic plane

Star clusters are distributed all over the Galactic plane, and each star
cluster creates a superbubble-like structure around itself (Weaver et al.,
1977; Gupta et al., 2018). Bronfman et al. (2000) observed $748$ OB
associations across the Galactic disk and found their distribution to peak at
$R_{p}=0.55\,R_{0}$, ($R_{0}=8.5$ kpc is the solar distance from the Galactic
center). We find that their inferred (differential) star cluster distribution
can be roughly fitted by

$dN_{c}(r)=\Sigma_{0}\,e^{\frac{-(r-R_{p})^{2}}{\sigma^{2}}}\,2\pi r\,dr,$ (1)

where $r$ is the Galactocentric distance and $\sigma=3$ kpc and $\Sigma_{0}$
is the normalization constant in the unit of kpc-2. This denotes the number of
star clusters in an annular ring of radius $r$ to $r+dr$. The surface density
of the clusters can be obtained by dividing this number by the surface area of
the annular ring ($2\pi r\,dr)$, as

$\displaystyle\nu(r)=\Sigma_{0}\,e^{\frac{-(r-R_{p})^{2}}{\sigma^{2}}}$ (2)

where, $\Sigma_{0}\sim 14$ kpc-2 (Nath & Eichler, 2020). We have used a
minimum number of 10 and a maximum number of $1000$ OB stars in a cluster
(these are somewhat arbitrary, and we later discuss the impact of these
choices).

The actual distribution of cosmic ray sources is expected to follow the
distribution of young stars and dense gas in the form of a spiral structure,
for instance, as traced by the FIR luminosity (Bronfman et al., 2000), or the
Ly-$\alpha$ radiation (Higdon & Lingenfelter, 2013), which show a bit of
asymmetry and trace the spiral arms to some extent. We emphasize that although
we assume an axisymmetric source distribution with smooth radial distribution,
this assumption yields a spectrum for cosmic-ray protons/nuclei that closely
resembles that derived from the spiral-arm feature of the source distribution
(e.g., Werner et al., 2015). Beyond $\sim 10$ GeV, Werner et al. (2015) show
that the spectrum’s shape remains largely unchanged and flux varies by less
than 2% when spiral-arm features are introduced. However, it is worth noting
that the presence of nearby star clusters associated with spiral arms can
introduce noticeable effects on cosmic-ray anisotropy at Earth. This is an
important topic for future investigation, but beyond the scope of the current
paper.

###  3.2 Transport of CRs originating from star clusters in the Galaxy

After getting accelerated in SNR and star cluster shocks, CRs propagate
through the Galaxy. This propagation is mainly dominated by diffusion and
energy loss due to interaction with ISM material. Some fraction of the
propagating CRs can be re-accelerated up to higher energy by the interaction
with existing weaker shocks that have been generated from older supernova
remnants in the ISM. This process has been discussed in detail in Thoudam &
Hörandel (2014). The transport equation for cosmic ray nuclei in a steady
state can be written as

$\displaystyle\nabla\cdot(D\nabla N)-[n{\rm v}\sigma+\zeta]\delta(z)N+$
$\displaystyle\Big{[}\zeta
sp^{-s}\,\int_{p_{0}}^{p}\,du\,N(u)\,u^{s-1}\Big{]}\delta(z)=-Q(r,p)\delta(z).$
(3)

Here we include spatial diffusion (first term on the left-hand side), re-
acceleration (terms with coefficient $\zeta$), and interaction losses
($\propto\sigma$, the loss cross-section) of the CR particles, as mentioned
above. The diffusion coefficient $D(p)$ depends on the momentum of CR
particles. Here $n$ represents the averaged surface density (number density
per unit area) of interstellar atoms in the Galaxy, $v(p)$ is the CR particle
velocity, $\sigma(p)$ is the cross-section of inelastic collision, $N$ is the
differential number density (number per unit volume per momentum) and $\zeta$
is the rate of re-acceleration. The third term involving the momentum integral
represents the generation of higher energy particles via the re-acceleration
of lower energy particles. It has been assumed that a CR population is
instantaneously re-accelerated to form a power-law distribution with an index
of $s\sim 4.5$ (Thoudam & Hörandel, 2014). We consider a cylindrical geometry
for the diffusion halo denoted by the radial coordinate $r$ and vertical
direction $z$. The diffusive halo has upper and lower boundaries at $z=\pm H$
and a radial boundary at $20$ kpc. A significant fraction of cosmic rays that
reach the earth is produced from those sources located within a distance $\sim
5$ kpc (Taillet & Maurin, 2003).

The term on the right side, $Q(r,p)\delta(z$), represents the injection rate
of cosmic rays per unit volume in the momentum bin $[p,p+dp]$ by the sources.
The $\delta(z)$ term denotes that all sources are confined to the Galactic
plane $z=0$. Similarly, re-acceleration and loss regions are confined within
the Galactic mid-plane.

The injection term $Q(r,p)$ can be written as a combination of a space-
dependent part and a momentum-dependent part, i.e.,

$\displaystyle Q(r,p)=\nu(r)\,H[R-r]\,H[p-p_{0}]\,Q(p),$

where $\nu(r)$ (see equation 2) represents the number of star clusters per
unit surface area on the Galactic disk (see section 3.1 for details),
$H[t]=1(0)$ for $t>0(<0)$ is the Heaviside step function, and $p_{0}$ (which
is the lower limit in the integral in Equation 3) is the low-momentum cutoff
introduced to approximate the ionization losses. Wandel et al. (1987) showed
that the ionization effects could be taken into account if we truncate the
particle distribution below $\sim 100$ MeV/nucleon. In our calculation, we
introduce a low-energy cutoff of $100$ MeV/nucleon. Our assumed distribution
of star clusters, motivated by observations, has a peak at $\sim 4.6$ kpc
($0.55R_{0}$, where $R_{0}$ is the distance of the Earth from the Galactic
center $\sim 8.5$ kpc) and then decreases rapidly at large distances (for
details see section 3.1).

The expression for surface density of star clusters $\nu$ has been calculated
in section 3.1, and the power-law source spectrum is described in section 3.3.
The energy-dependent diffusion coefficient as a function of particle rigidity
follows

$\displaystyle D(\rho)=D_{0}\,\beta(\rho/\rho_{0})^{\delta},$

where $D_{0}$ is the diffusion constant, $\rho=Apc/Ze$ is the particle
rigidity, $\beta=v(p)/c$ where $v(p)$ is the CR particle velocity and $c$ is
the speed of light, $\delta=0.33$ is the diffusion index, and $\rho_{0}=3$ GV
is a normalisation constant.

In this injection-diffusion-reacceleration111For typical parameters,
reacceleration only affects the CR spectrum below $10^{5}$ GeV (Thoudam &
Hörandel, 2014) and so is irrelevant for the energy range considered in this
work. model, the rate of reacceleration depends on the rate of supernova
explosions and the fractional volume occupied by SNRs in the Galaxy. The
reacceleration parameter $\zeta$ can be expressed as, $\zeta=\eta V\nu_{SN}$,
where $V=4\pi\mathfrak{R}^{3}/3$ is the volume occupied by each SNR of radius
$\mathfrak{R}$ re-accelerating the cosmic rays. Here, $\eta$ is a correction
factor that takes care of the actual unknown size of the remnants, and
$\nu_{SN}$ is the rate of supernova explosions per unit surface area in the
Galactic disk. The values of $\mathfrak{R}$ and $\nu_{SN}$ have been taken as
$100$ pc and $25$ SNe Myr-1kpc-2 respectively (Thoudam & Hörandel, 2014).

Figure 2: Top: Schematic distribution of star clusters in the Galactic plane
(face-on view), each star indicates a star cluster on the plane; bottom: the
surface density (number per area) of star clusters ($\Sigma$; see Eq. 2) as a
function of distance from the Galactic center.

The solution of equation 3 can be obtained by invoking the Green’s function
method and by considering the two separate transport equations for the regions
below and above the Galactic disk ($z<0$ and $z>0$ respectively), and by
connecting the two solutions at Galactic plane, i.e., $z=0$, via a jump
condition. Following this procedure one can get the Green’s function
$G(r,r^{\prime},z,p,p^{\prime})$ (equation A.20 of Thoudam & Hörandel 2014).
After convolving the obtained Green’s function with the assumed source
distribution and integrating it over the Galactic plane, one can get the final
solution (see equation 6 of the same paper) for the CR density $N(r,z,p)$.
Following the procedure described in Thoudam & Hörandel (2014), we get the
solution of the transport equation 3 as

$\displaystyle
N(r,z,p)=2\pi\,\int_{0}^{\infty}dp^{\prime}\int_{0}^{\infty}r^{\prime}dr^{\prime}\,G(r,r^{\prime},z,p,p^{\prime})\,Q(r^{\prime},p^{\prime}).$

Substituting the obtained $G(r,r^{\prime},z,p,p^{\prime})$ (Thoudam &
Hörandel, 2014) and the assumed source distribution in the above equation, the
cosmic ray density at the Earth ($r=8.5$ kpc) can be calculated by evaluating
the above solution at $z=0$ since our Solar system lies close to the Galactic
plane. More explicitly, the differential number density measured at the
location of Earth is

$\displaystyle N(r,p)$ $\displaystyle=$
$\displaystyle\int_{r^{\prime}=0}^{R}\,\int_{k=0}^{\infty}\Sigma_{0}\frac{J_{0}[k(r-r^{\prime})]}{L(p)}\,k\,dk\,e^{\frac{-(r^{\prime}-R_{p})^{2}}{\sigma^{2}}}\,r^{\prime}\,dr^{\prime}$
$\displaystyle\Big{[}Q(p)+\zeta
sp^{-s}\,\int_{p_{0}}^{p}\,{p^{\prime}}^{s}\,dp^{\prime}\,Q(p^{\prime})\,A(p^{\prime})\,$
$\displaystyle\times\exp\Big{(}\zeta
s\int_{p^{\prime}}^{p}\,A(u)du\Big{)}\Big{]},$

where $R=20$ kpc is the radial boundary of the Galaxy, $\Sigma_{0}$ is the
number density of star clusters, $J_{0}$ is the Bessel function of order zero
and the functions $A(p)$ and $L(p)$ are given by (see Thoudam & Hörandel 2014
for details),

$\displaystyle L(p)$ $\displaystyle=$ $\displaystyle
2D(\rho)k\,\coth(kH)+nv\sigma(p)+\zeta\,,$ (5) $\displaystyle A(p)$
$\displaystyle=$ $\displaystyle\frac{1}{pL(P)}.$ (6)

Equation LABEL:eq:number_density gives the differential number density, i.e.,
number per unit volume per unit momentum of cosmic ray particles measured at
earth. All the necessary terms needed to solve equation
LABEL:eq:number_density are discussed in sections 3.2–3.5.

###  3.3 Injection spectra of cosmic ray nuclei

The cosmic ray source spectrum $Q(p)$ from star clusters is assumed to follow
a power-law in total momentum $Ap$, with an exponential cut-off, where $A$ is
the mass number of the nucleus. We write the differential number of CR
particles with nucleon number $A$, having momentum per nucleon in the range
$(p\,,p+dp)$, as,

$Q(p)=Q_{0}(Ap)^{-q}\exp\left(-\frac{Ap}{Zp_{\rm max}}\right)\,.$ (7)

Here $Q_{0}$ is a normalization constant that is proportional to the fraction
of total wind kinetic energy $f$ channeled into cosmic rays by a single star
cluster. We call this ‘injection fraction’, which is a free parameter and can
be estimated by comparing the model result with observations. Also, $q$ is the
spectral index, $p_{\rm max}$ is the cutoff momentum (for a single nucleon),
and $Ap$ is the total momentum of a particle with the mass number $A$ and the
atomic number $Z$.

###  3.4 Maximum energy estimate of accelerated particles

For the estimation of the maximum accelerated energy of cosmic ray particles,
we consider two different acceleration scenarios inside a young star cluster:
acceleration at WTS and acceleration of particles around SNR shock inside a
star cluster.

####  3.4.1 Acceleration at wind termination shock (WTS):

In equation 7, $p_{\rm max}$, which represents the maximum momentum of
accelerated CRs, depends on the extension of the accelerating region for a
stellar wind bubble of the cluster. Typically this accelerating region can be
taken as the distance to the WTS ($R_{\rm WTS}$) from the center of the
cluster. The maximum energy is achieved when the diffusion length becomes
comparable to the size of the shock (in this case the WTS), for beyond this
limit, the particles escape out of the accelerating region. In the case of
Bohm diffusion, $\kappa=pc^{2}/(\zeta qB)$, the maximum energy is then (Vieu
et al., 2022):

$\displaystyle E_{\rm max}\sim\zeta\,q\,B_{\rm WTS}\,{R_{\rm
WTS}}\,\frac{V_{w}}{c}\,.$ (8)

Here $R_{\rm WTS}$ is the radius of WTS. In the above equation, $V_{w}$ is the
velocity of stellar wind, $B_{\rm WTS}$ is the value of the magnetic field at
the WTS position, $\zeta=3r_{g}/\lambda$, with $\lambda$ the mean free path
due to the magnetic field. The Bohm diffusion, which is the most optimistic
scenario, corresponds to the limit $\zeta=3$.

We follow the arguments advocated by Vieu et al. (2022) to estimate the
magnetic field in the cluster core,

$\displaystyle B_{c}\sim 150\left(\frac{n_{\rm c}}{10\,{\rm
cm^{-3}}}\right)^{1/6}\left(\frac{\eta_{T}}{0.1}\right)^{1/3}\left(\frac{N_{\rm
OB}}{100}\right)^{2/9}$ $\displaystyle\left(\frac{R_{c}}{1{\rm
pc}}\right)^{-2/3}\,\mu G\,.$ (9)

Here, $n_{\rm c}$ is the core density, $\eta_{T}$ is the efficiency of
generation of turbulence, $N_{\rm OB}$ is the number of OB stars in the
cluster, $R_{c}$ is the core radius of the cluster. The magnetic field
advected into the free wind region has a $1/r$ radial profile. Therefore, the
magnetic field at the position of the wind termination shock and cluster core
can be related using $B_{c}R_{c}=B_{\rm WTS}R_{\rm WTS}$. Therefore, equation
8 can be expressed as,

$\displaystyle E_{\rm max}\sim\zeta\,q\,B_{c}\,{R_{c}}\,\frac{V_{w}}{c}\,.$
(10)

This leads to a maximum estimate:

$\displaystyle E_{\rm max}\sim 6\left(\frac{\zeta}{3}\right)\left(\frac{n_{\rm
c}}{10\,{\rm
cm^{-3}}}\right)^{1/6}\left(\frac{\eta_{T}}{0.2}\right)^{1/3}\left(\frac{N_{\rm
OB}}{1000}\right)^{2/9}$ $\displaystyle\left(\frac{R_{c}}{1~{}{\rm
pc}}\right)^{1/3}\,\left(\frac{v_{w}}{2000\,{\rm km\,s^{-1}}}\right)\,{\rm
PeV}$ . (11)

Equation 11 gives a conservative estimate of ${\rm E_{max}=6\,PeV}$ ($6\times
10^{6}$ GeV) for the maximum attainable energy of protons. Note that this
value is a few times higher than the maximum accelerated energy for the SNR-CR
scenario.

However, in the realistic scenario, the magnetic field may be amplified in the
accelerating region due to the existence of turbulence, and due to
instabilities driven by cosmic ray streaming in the upstream region of WTS,
which can therefore increase the estimated value of maximum accelerated
energy. There are other uncertainties as well (e.g., in $\eta_{T}$, wind
velocity $v_{w}$) that can conceivably increase the maximum energy by a factor
of a few.

####  3.4.2 Acceleration at SNR shock inside star clusters:

Another potential scenario for CR acceleration inside the young star cluster
is the SNR shocks propagating in the free wind region of the cluster. CR
particles can be accelerated up to $10^{8}$ GeV if the SNR shocks advance
through fast and highly magnetised stellar winds (Voelk & Biermann, 1988;
Biermann & Cassinelli, 1993). Non-linear effects in the acceleration process
(Bell & Lucek, 2001) also contribute to this high-energy acceleration. Bell
(2013); Schure & Bell (2013) highlight that the outer shocks of SNR can
accelerate CR beyond the ‘knee’ if the shock propagates into a magnetic field
much larger than a typical interstellar field, that can be present inside a
star cluster. Particles will be accelerated during the expansion of the SNR
shock in upstream of the WTS. This idea has been studied extensively by Vieu
et al. (2022), and the maximum energy has been estimated by the authors as
follows:

$\displaystyle E_{\rm max}\sim 21\,\left(\frac{V_{c}}{5000\,{\rm
km\,s^{-1}}}\right)\,\left(\frac{\zeta}{3}\right)\,\left(\frac{R_{c}}{1\,{\rm
pc}}\right)\,\left(\frac{N_{\rm OB}}{1000}\right)^{2/9}\,$
$\displaystyle\left(\frac{n_{\rm c}}{10\,{\rm
cm^{-3}}}\right)^{1/6}\left(\frac{\eta_{T}}{0.2}\right)^{1/3}\,\left[1-\left(\frac{R_{c}}{R_{\rm
WTS}}\right)^{1/7}\right]\,\,\,\,{\rm PeV}\,\,.$

For a typical young cluster $R_{\rm WTS}/R_{c}\sim 5-30$, which gives
$\left(1-(R_{c}/R_{\rm WTS})^{1/7}\right)\sim 0.2-0.4$ (Vieu et al., 2022).
This estimate can give a maximum energy of a few PeV for protons. However, if
a supernova launches a very fast shock in the free wind region of a compact
cluster with velocity $\geq 2\times 10^{4}$ km s-1, it can accelerate protons
up to a few tens of PeV energy. Note that, this high velocity of SNR shock
inside a clumpy star cluster may efficiently drive MHD turbulence to generate
a high value of the magnetic field, which will likely result in a higher value
of maximum energy.

Figure 3: Model prediction for the star cluster model as a second galactic
component considering an injection fraction $\sim 5\%$. The thick solid maroon
line represents the total contribution from Galactic star clusters. Thin
dashed lines represent the flux of individual elements. For the CRs generated
from star clusters, an exponential energy cut-off for protons at
$E_{c}=5\times 10^{7}$ GeV ($50$ PeV) is assumed. High-energy data: IceTop
(Aartsen et al., 2013), Tibet III (Amenomori et al., 2008), the Pierre Auger
Observatory (The Pierre Auger Collaboration et al., 2013), and HiRes II (High
Resolution Fly’S Eye Collaboration et al., 2009). Low energy data have been
taken from CREAM (Ahn et al., 2009; Yoon et al., 2011), ATIC-2 (Panov et al.,
2007), AMS-02 (Aguilar et al., 2015), PAMELA (Adriani et al., 2011), CRN
(Mueller et al., 1991), HEAO (Engelmann et al., 1990), TRACER (Obermeier et
al., 2011), KASCADE (Antoni et al., 2005), DAMPE (An et al., 2019). We have
only shown the high-energy data points with different symbols in the figure.
Low data points: Proton (black square), Helium (grey square), Oxygen (purple
solid plus), Carbon (red circle), Iron (blue circle), Neon (green circle),
Silicon (magenta circle), Magnesium (black stars). The lower energy data from
various experiments are represented together by one symbol. The error bars for
proton and helium have been shown and the rest are not shown in the figure.

The recent detection of $\gamma$-rays above PeV by LHAASO from some sources
indeed indicates these sources can accelerate particles up to at least a few
tens of PeV because, at high energy, the $\gamma$-ray energy can be
approximated as $E_{\rm cr}\approx 10E_{\gamma}$. Some of those sources
possibly are young massive clusters (see extended table 2 of (Cao et al.,
2021). The $\gamma$-ray photons from the LHAASO J2032+4102 source have the
highest energy of 1.4 PeV, which corresponds to tens of PeV for cosmic ray
protons.

###  3.5 Elemental abundances in star cluster winds

We consider a simple stellar population formed at time $t=0$ with an initial
mass function (IMF) $\frac{dn}{dm}\propto m^{-2.35}$ (Salpeter, 1955). We can
calculate the elemental abundances in the wind material following the
procedure described in Roy et al. (2021). Now,

$M_{w}(X,m,t)=\int_{0}^{t}\,\dot{m}_{w}(X,m,t^{\prime})\,dt^{\prime}$ (13)

is the cumulative mass of element $X$ ejected in winds by a star of initial
mass $m$ up to age $t$ where,

$\displaystyle\dot{m}_{w}(X,m,t^{\prime})$ $\displaystyle=$ $\displaystyle{\rm
mass\,\,\,fraction}(X,m,t^{\prime})\times\dot{m}(m,t^{\prime})$ (14)
$\displaystyle=$ $\displaystyle f(X,m,t^{\prime})\times\frac{dm_{\rm
star}}{dt^{\prime}}.$

We use the mass loss rate for each nucleus $\dot{m}(X,m,t^{\prime})$ using
models for nucleosynthesis in massive stars and their return to the ISM via
winds (A. Roy, private communication). Hence the elemental abundance of a
particular element $X$ can be calculated using,

$\displaystyle f(X,m)=\frac{M_{w}(X,m,t)}{M_{w,\rm
tot}(m,t)}=\frac{\int_{0}^{t}\,\dot{m}_{w}(X,m,t^{\prime})\,dt^{\prime}}{\int_{0}^{t}\,\dot{m}(m,t^{\prime})\,dt^{\prime}}.$
(15)

We have taken into account evolution until the core carbon burning time, which
implies the maximum time of the evolution of a star with mass $m$ as the upper
limit of the integration. The mass-weighted elemental abundance of element X
can be calculated using the following expression invoking the Salpeter mass
function,

$\displaystyle\langle
f(X)\rangle=\frac{\int_{0}^{m}\,f(X,m)\,{Am^{-2.35}}dm}{\int_{0}^{m}\,{Am^{-2.35}}dm}$
(16)

Using this method, we have calculated the mass-weighted mean individual
elemental abundance in the ejected stellar wind material. We have used the
results of a state-of-the-art evolutionary model. The elemental abundances
have been calculated considering the rotation-driven instabilities inside the
star, the correct abundances of elements, and the mass loss rate from the
stellar surface.

Figure 4: Model prediction for the all-particle spectrum using the Galactic
star cluster CR model as the second galactic component. For the star cluster
component, the considered injection fraction is $\sim 5\%$, and the cutoff is
at $5\times 10^{7}Z$ GeV. The thick solid maroon line represents the total
SNR-CRs, the thick dashed maroon line represents star cluster CRs, and the
thick maroon dotted line represents the UFA model of extragalactic CR
component (EG-UFA) taken from Unger et al. (2015), and the thick solid blue
line represents the total all-particle spectrum. The thin lines represent the
total spectra for the individual elements i.e., a combination of both SNR-CR
and the CRs originating from star clusters. The figure shows the $E^{3}$ times
the cosmic ray flux $I(E)=(c/4\pi)N(E)$ at the position of the earth measured
by different experiments as a function of cosmic ray energy, where $N(E)$ is
the differential number density of cosmic ray particles. High energy and low
energy data are the same as figure 3.

###  3.6 Average kinetic luminosity of clusters:

Our assumption requires a certain fraction of the total wind kinetic energy to
go into CRs. Therefore, we need the value of the average kinetic luminosity of
a cluster using a distribution of OB associations over the luminosity range.
Oey & Clarke (1997) assume that the mechanical luminosity function of OB
association is given by $\phi(L)\propto L^{-2}$. We use this distribution to
calculate the average luminosity of clusters with kinetic luminosity in the
range $L_{\rm min}=10^{37}$ erg s-1 (corresponds to $N_{\rm OB}=10$) to
$L_{\rm max}=10^{39}$ erg s-1 (corresponds to $N_{\rm OB}=1000$) as following,

$\displaystyle\langle L_{w}\rangle$ $\displaystyle=$
$\displaystyle\frac{\int_{L_{\rm min}}^{L_{\rm
max}}\phi(L)\,L\,dL}{\int_{L_{\rm min}}^{L_{\rm max}}\phi(L)\,dL}\sim
4.5\times 10^{37}\,\,{\rm erg\,s^{-1}}.$ (17)

Note, we adopt a minimum $10$ number of OB stars for the production of CRs,
and the largest OB association in our Galaxy has $1000$ OB stars. The
dependence of $\langle L_{w}\rangle$ on $L_{\rm max}$ is weak, but there is
the sensitive dependence on $L_{\rm min}$, the implications of which we
discuss later.

Figure 5: Mean logarithmic mass $\langle\ln A\rangle$ of cosmic rays as a
function of energy, predicted using the combination of SNR-CR, CRs from star
clusters (these two are Galactic components), and EG-UFA model (extragalactic
component, Unger et al. 2015). Data: KASCADE (Antoni et al. 2005), TUNKA
(Berezhnev 2015), Yakutsk (Knurenko & Sabourov 2011), the Pierre Auger
observatory (The Pierre Auger Collaboration et al. 2015) and the different
optical measurement compiled in Kampert & Unger 2012. The two different
colored (black and grey) sets of data points correspond to two models EPOS-LHC
and QGSJET-II-04, respectively, which have been used to convert $X_{\rm max}$
values to ${\langle\ln A\rangle}$ (see equation 18).

##  4 Model prediction for the second component of galactic cosmic rays

The values of cosmic ray propagation parameters ($D_{0},\,\delta$; the
normalization of the diffusion coefficient and its power-law dependence on
momentum) and re-acceleration parameters ($\eta\,,s$; the SNR filling factor
and reacceleration power-law index) have been calculated by comparing the
observed Boron to Carbon abundance ratio with the value obtained by the
adopted model. The best fit values are $D_{0}=9\times 10^{28}$ cm2s-1,
$\eta=1.02$, $s=4.5$, and $\delta=0.33$ (Thoudam & Hörandel, 2014). We have
also used these values in our model. For the interstellar matter density
($n$), the averaged density in the Galactic disk within a radius equal to the
size of the diffusion halo $H$ was considered. We choose $H=5$ kpc (Thoudam et
al., 2016), which gives an averaged surface density of atomic hydrogen of
$n=7.24\times 10^{20}$ atoms cm-2 (Thoudam & Hörandel, 2014). To account for
the helium abundance in the interstellar medium, we add an extra $10$% to $n$.
The radial extent of the source distribution is taken as $R=20$ kpc. The
inelastic cross-section for proton ($\sigma(p)$) is taken from Kelner et al.
(2006).

Since we are interested in the acceleration of CRs in WTS, as well as around
SNR shocks inside the free wind region of the star cluster, the relevant
abundances correspond to that in the stellar wind for massive stars. For this
purpose, we use the stellar wind abundances for massive stars beginning with
the Zero Age Main Sequence (ZAMS) phase. We have used the surface abundance of
massive stars as a function of time, calculated after properly taking into
account the effect of stellar rotation. The spectral indices for different
elements are given in Table 1. Note that these values are slightly different
from the spectral indices assumed in Thoudam et al. (2016) for the SNR-CR
component. Also, the stellar wind elemental abundances are mentioned in Table
1, which are calculated using the method described in Roy et al. (2021)
(provided to us by A. Roy, private communication). We have then averaged the
abundances over time and mass distribution of stars in the cluster, as
described in section 3.5. Using these values of various parameters, we
calculate the particle spectra for different cosmic ray elements (proton,
helium, carbon, oxygen, neon, magnesium, silicon, and iron). The CR spectral
indices ($q$) of source spectra for the individual elements are very similar
to each other and are chosen to match the observed individual nuclear
abundances in CRs closely.

Elements | $q$ | Fractional abundances in winds
---|---|---
Proton | 2.25 | 0.86
Helium | 2.23 | 0.13
Carbon | 2.20 | $3.32\times 10^{-3}$
Oxygen | 2.24 | $8.51\times 10^{-4}$
Neon | 2.24 | $8.83\times 10^{-5}$
Magnesium | 2.28 | $3.62\times 10^{-5}$
Silicon | 2.24 | $3.42\times 10^{-5}$
Iron | 2.24 | $3.72\times 10^{-5}$
Table 1: Source spectral indices $q$ and fractional abundances of different
elements in the wind material. The elemental abundances are calculated
following Roy et al. (2021).

Figure 3 shows the star cluster contribution to CRs using different parameters
mentioned earlier. We have used the maximum energy for the proton as $5\times
10^{7}$ GeV ($50$ PeV) and the injection fraction of $\sim 5\%$. These values
of the parameters are chosen to match the observed spectra with our
theoretical model. It is important to mention that, in section 3.4 we have
estimated the maximum accelerated energy considering different scenarios in a
star cluster. The maximum energy can go up to a few tens of PeV (especially
for the SNR shock inside the star cluster scenario), although our used value
is admittedly on the higher side. Also, recently Vieu & Reville (2023) have
shown that the SNR shocks inside a star cluster scenario can explain the all-
particle cosmic ray spectrum in the region between ‘knee’ and ‘ankle’.
Therefore, star clusters are likely a possible candidate for cosmic ray
acceleration between a few times $10^{6}$ and $10^{9}$ GeV.

Also, if one uses a higher lower limit of $N_{\rm OB}=30$ instead of $10$,
then the injection fraction will need to be increased to match the observed
spectrum. The data points correspond to different measurements. For lower
energy ranges, the individual spectra are fitted to the observed elemental
spectra. We consider $8$ elements: proton, helium, carbon, oxygen, neon,
magnesium, silicon, and iron for our calculations, and the total contribution
(solid brown curve in the figure 3) is a combination of these $8$ elements.

##  5 All-particle spectrum of cosmic rays

Figure 4 combines all three CR components to get the total all-particle
spectrum of cosmic rays and compares it with various observations. The SNR-CR
component shown in this figure is calculated following the procedure mentioned
in Thoudam et al. (2016), assuming a uniform distribution of SNRs in the
Galactic plane and a proton spectrum cut-off of $\sim 2.5\times 10^{6}$ GeV.
For the extragalactic component, we have adopted the UFA model (Unger et al.,
2015), which considers a significant contribution of extragalactic CRs below
the ankle to reproduce the observed CR energy spectrum as well as $X_{\rm
max}$ (the depth of the shower maximum) and the variance of $X_{\rm max}$
above the ankle observed at the Pierre Auger Observatory (di Matteo, 2015).
With these two models (SNR & extragalactic), we have combined our proposed
star cluster model with a proton spectrum cut-off at $5\times 10^{7}$ GeV
($50$ PeV).

The total contributions from all these three components can explain the
observed features in the all-particle spectrum. Also, the spectra of the
individual elements can be explained well with the model. The flux of
different elements has been measured well in the lower energy region, but in
the higher energy range, i.e., above $10^{5-6}$ GeV, the observation data for
individual elements are not available. Observed data points for all-particle
CR spectra have been taken from various experiments like TIBET III (Amenomori
et al., 2008), IceTop (Aartsen et al., 2013), Auger (The Pierre Auger
Collaboration et al., 2013), HiRes II (High Resolution Fly’S Eye Collaboration
et al., 2009), etc.

Several ground-based experiments such as KASCADE, TUNKA, LOFAR, and the Pierre
Auger Observatory have provided measurements of the composition of CRs at
energies above $\sim 10^{6}$ GeV. Heavier nuclei interact at a higher altitude
in the atmosphere, which results in smaller values of $X_{\rm max}$ as
compared to lighter nuclei. For comparison with the theoretical predictions,
$\langle\ln A\rangle$, the mean logarithmic mass of the measured cosmic rays,
is of utmost importance. This can be obtained from the measured $X_{\rm max}$
values using the following relation mentioned in Hörandel (2003),

$\displaystyle{\ln\,A}_{i}\,=\,\Big{(}\frac{X_{\rm max}^{i}-X_{\rm max}^{\rm
p}}{X_{\rm max}^{\rm Fe}-X_{\rm max}^{\rm p}}\Big{)}\,\times{\ln\,A_{Fe}}.$
(18)

Here $X_{\rm max}^{\rm p}$ and $X_{\rm max}^{\rm Fe}$ represent the average
maximum depths of the shower for protons and iron nuclei, respectively, and
${\rm A_{Fe}}$ is the mass number of iron nuclei. In figure 5, we have also
shown the obtained mean logarithmic mass using our model and compared it with
the observational data.

We calculate the mean mass in the following way,

$\displaystyle\langle{\ln\,A}\rangle=\frac{{\sum_{i}\ln\,A_{i}}\times{\rm
Flux_{i}}}{\sum_{i}\,{\rm Flux_{i}}}$ (19)

where ${A_{i}}$ denotes the mass number of an element $i$ (we have considered
$8$ elements: proton, helium, carbon, oxygen, neon, magnesium, silicon, and
iron), and ${\rm Flux_{i}}$ is the obtained flux of element $i$ using our
model. Figure 5 shows that the results obtained using our star cluster model
(green curve) follow the observed trend for the mean logarithmic mass in the
total energy range from $10^{8}$ GeV to $10^{11}$ GeV when combined with the
UFA model for the extragalactic CRs. In the energy range of about $2\times
10^{7}$ and $10^{8}$ GeV, our prediction shows some deviation from the
observed trend but still lies within limits presented in Kampert & Unger
(2012).

To reiterate, the primary focus of this work has been to present a model
incorporating stellar wind shocks that can explain the observed all-particle
CR spectrum, especially in the energy range between $10^{7}$ and $10^{9}$ GeV.
The discussion so far shows that we can indeed explain the observed data in
this energy range using a CR component originating from massive star clusters.
The required CR injection fraction for this component of $\sim 5\%$ and an
energy cutoff of $5\times 10^{7}\,Z$ GeV ($50$ PeV), as suggested by the
fitting of the all-particle CR spectrum with our proposed stellar wind model,
are entirely reasonable, and therefore lend support to the idea that CR from
massive star clusters can fill the CR spectrum gap between the ‘knee’ and the
‘ankle’. We will discuss the implications of our result in Section 7. Before
that, we discuss the dependence of the required injection fraction and cutoff
energy on the chosen extragalactic component in section 6 below.

##  6 Varying the extragalactic component

Figure 6: Top panel: All-particle CR spectrum when combined with SNR-CRs and
EG-MPCS model (Rachen 2015) for the extragalactic CRs. Bottom panel: Mean
logarithmic mass when combined with the EG-MPCS (red curve) and the EG-UFA
(green curve, same as 5 ) models. Data are the same as in Figure 5.

As mentioned in 2.2, we consider two different models of extragalactic CRs:
UFA model (Unger et al., 2015) and a combination of PCS and Minimal model
(MPCS model) (Rachen, 2015; Thoudam et al., 2016). Depending on the chosen
extragalactic component, the value of injection fraction and maximum cutoff
energy can slightly change. The UFA and MPCS models predict a significant
contribution of extra-galactic cosmic rays below the ‘ankle’. All these
different extragalactic models can explain the observations when combined with
the SNR-CR component and the CR component from star clusters, although the UFA
model somehow shows a smooth transition (Figure 4) between the Galactic and
extragalactic components. The sharp increase near $10^{9}$ GeV in the MPCS
model (top panel, figure 6) is due to the dip in the proton spectrum. It
results from the intersection of the minimal model and the components from
galaxy clusters. Below $10^{9}$ GeV, both the UFA and MPCS models give similar
results and can explain the observed spectra. We have also shown the mean
logarithmic mass plot for a combination of each different extragalactic model
with the two different Galactic components (bottom panel of figure 6). It is
clear from the plot that all these different models for the extragalactic
component, in combination with the Galactic components, follow the observed
trend of mean logarithmic mass for the whole energy range.

##  7 Discussion

Our study demonstrates that the cosmic rays originating from spatially
distributed young massive star clusters in the Galactic plane fit well the
all-particle CR spectrum, particularly in the $10^{7}-10^{9}$ GeV energy
range, and therefore this can be a potential candidate for the ‘second
Galactic component’ of CRs. We also show that the observed all-particle
spectrum, as well as the cosmic ray composition at high energies, can be
explained with the following two types of Galactic sources: (i) SNR-CRs,
dominating the spectrum up to $\sim 10^{7}$ GeV, and (ii) star cluster CRs,
which dominate in the range $10^{7}-10^{9}$ GeV.

The SNR-CR component can only contribute up to maximum energy, corresponding
to a proton cut-off energy of $2.5\times 10^{6}$ GeV (Thoudam et al., 2016).
Such a high value of energy cannot be achieved if we consider the DSA
mechanism with typical values of the magnetic field in the ISM. However, some
numerical simulations have indicated that supernova shocks can amplify the
magnetic field near them several times larger than the value in the ISM (Bell
& Lucek, 2001; Reville & Bell, 2012). Such a strong magnetic field can
accelerate CR protons up to the cut-off energy used in this study. Also,
recently detected $\gamma$-rays from a few SNRs have also identified a few
SNRs as cosmic ray PeVatrons that can accelerate particles up to a few times
$\sim 10^{6}$ GeV energy.

Maximum CR energy in star clusters: According to our model, the component of
CRs that is plausibly generated in star clusters can contribute significantly
towards the total CR flux, especially in the $10^{7}-10^{9}$ GeV range, if one
considers that the protons can be accelerated up to $5\times 10^{7}$ GeV
energy. For other elements with atomic number $Z$, the maximum energy is
$5\times 10^{7}$ Z GeV in these young compact star clusters, and a cosmic ray
injection fraction of $\sim 5\%$.Note that this value of maximum energy
required for proton is slightly on the higher side, but can be justified under
the assumption of the very high initial shock velocity of SNRs inside compact
clusters, faster wind velocity, and possible amplification of magnetic field
inside the cluster. Our maximum CR energy from Hillas criterion may be an
overestimate, but we should point out that the requirement of $E_{\rm max}=50$
PeV ($5\times 10^{7}$ GeV) depends on the assumption of elemental abundance
ratios in the wind material, an aspect that remains uncertain at present. A
higher abundance of heavy elements would increase $E_{\rm max}$, as required
to fit the observed spectrum. Also note that the recently detected
$\gamma$-ray photons by LHAASO in PeV range from $12$ objects, some of which
are associated with massive star clusters, have indicated that these clusters
can accelerate particles at least up to a few tens of PeV (Cao et al., 2021),
consistent with our estimates.

CR anisotropy: The total CR anisotropy $\Delta$ can be calculated from our
model using the diffusion approximation given by Mao & Shen (1972),

$\displaystyle\Delta=\frac{3D}{c}\frac{|\nabla N|}{N}\,\,\,,$ (20)

where $N$ is the CR number density, $D$ is the diffusion coefficient and $c$
is the velocity of light. From our model, we get $\Delta\sim 0.04-0.2$ in the
range $1-100$ TeV. However, it is noteworthy that our calculated estimates
exhibit higher values compared to the measured anisotropy, which is
approximately in the range of $(0.5-1)\times 10^{-3}$ for the same energy
spectrum. Notably, our findings align with the earlier estimates proposed by
Blasi & Amato 2012; Thoudam & Hörandel 2012, although, like the case of
previous calculations, they are larger than the observed anisotropy in the
same energy range.

Inside an OB association, individual SNR shocks, as well as colliding shocks,
can accelerate particles on a time scale below $1000$ years. An OB association
may enter the evolutionary stage of multiple SN explosions on a time scale
larger than a few hundred thousand years. It can create large bubbles of $\sim
50$ pc size, and the injected mechanical power can reach $\sim 10^{38}$ erg
s-1 over $10$ Myr—the lifetime of a superbubble. This process is supplemented
by the formation of multiple shocks, large-scale flows, and broad spectra of
MHD fluctuations in a tenuous plasma with frozen-in magnetic fields. The
collective effect of multiple SNRs and strong winds of young massive stars in
a superbubble is likely to energize CR particles up to hundreds of PeV in
energy (see Montmerle 1979; Cesarsky & Montmerle 1983; Bykov & Fleishman 1992;
Axford 1994; Higdon et al. 1998; Bykov & Toptygin 2001; Marcowith et al. 2006;
Ferrand & Marcowith 2010) and even to extend the spectrum of accelerated
particles to energies well beyond the ‘knee’ (Bykov & Toptygin, 2001).

Gupta et al. (2020) pointed out the advantage of considering CRs accelerated
in massive star clusters in explaining several phenomena (e.g., Neon isotope
ratio) that are left unexplained by the paradigm of CR production in SNRs.
They also proposed that this component need not be considered entirely
independent and separate from the SNR component since massive stars (which are
the progenitors of SNRs) always form in clusters. Therefore, the two
components (SNR-CR and star cluster CRs) arise from similar sources, with some
differences. The SNR-CRs can be thought of as CRs produced in individual SNRs,
which arise from very small clusters with one or two massive stars, whereas
the second component can be thought of as arising from different shocks that
occur in the environment of massive star clusters. Hence, the two components
can be put on the same platform, and the combined scenario offers a fuller,
more complete picture of the phenomenon of CR acceleration in the Galaxy.

###  7.1 Caveats of our model

Finally, we discuss some caveats of our model.

The injection fraction, a free parameter in our analysis whose value is
obtained by fitting the all-particle CR spectrum, ultimately depends on many
other factors, such as diffusion coefficient, the assumed lower limit of the
number of OB stars in a cluster, and so on. With a higher value of diffusion
coefficient, the required injection fraction should be increased in order to
match our results with observations. A larger diffusion coefficient implies
that particles would diffuse out of the source more rapidly, which will
decrease the particle density in the vicinity of clusters. This is why one
needs a larger injection fraction to explain the observational data. On the
other hand, the total number of OB associations depends on the lower cut-off
in the distribution of cluster masses. E.g., the number of OB associations
which has a minimum of 30 OB stars is lower than the number of OB associations
with a minimum of 10 OB stars. For the second case, the required injection
fraction will be lower (since the number of OB associations is higher). Also,
the location of the peak of the cluster spatial distribution has a significant
effect on the observed flux and may introduce some uncertainty to the value of
the injection fraction parameter.

The efficiency is likely independent of the number of OB stars inside a single
star cluster. Considering the gamma-ray luminosity of massive star clusters,
observations indicate the gamma-ray luminosity is $\sim 10^{-3}\times L_{w}$
(where $L_{w}$ the wind kinetic energy), irrespective of the total Stellar
number (e.g., Ackermann et al. 2013), as we have assumed here (see also Gupta
et al. 2018). However, the required cosmic ray injection fraction to explain
the observed data depends on the total number of star clusters in the Galaxy.
For the current study, we have considered the spatial distribution of OB stars
but we did not classify them according to their age and mass. According to the
recent GAIA survey, $\sim 20\%$ of the total clusters are compact young, and
massive (Vieu et al., 2022). However, we have assumed all the OB associations
are young and compact. If we instead take a fraction (say $20\%$) of these
clusters to be young then we would need a higher (by a factor of 5, say)
fraction of cosmic ray injection efficiency to match the observed flux.

Elemental abundance: There are other uncertainties that arise from the assumed
abundances of the eight elements considered here. This elemental abundance
depends on the rotational velocity of the stars. For our calculations, we have
used the abundances of stars, which rotate with a velocity that is $60\%$ of
the critical velocity of the star. Varying the rotational velocity would
change the elemental abundances. This will give an uncertainty between $2-3\%$
in the mean logarithmic mass plot (Figure 5). Results may also change if
abundances from other previous works (Heger et al. 2000, 2005) were to be
used. Comparing the abundances from these works with the ones used here, we
find that it would introduce an uncertainty of $5-7\%$ in Figure 5. However,
these variations will not significantly change the shape of our predicted
$\langle\ln A\rangle$.

CR propagation: Another aspect that is important to mention is the mode of
cosmic ray propagation. Our calculation of the second component assumes
diffusion from the source. Nevertheless, it is crucial to acknowledge that the
diffusion approximation may cease to be valid beyond a specific energy
threshold, leading to a shift from diffusion to drift motion in the transport
process. If we consider a mean magnetic field of $3$ $\mu$G in the Galactic
plane then the transition from diffusion to drift will occur at $\sim Z\times
10^{17-18}$ eV (Kääpä et al., 2023). The maximum energy for proton in WTS is
$\sim 50$ PeV ($5\times 10^{16}$ eV) and is below the transition region.
Therefore, the diffusion approximation works well for the energy range
considered by us. However, we modified the diffusion coefficient above
$10^{17}$ eV in a manner that mimics ballistic propagation beyond this energy
threshold and found the spectra do not change significantly as the second
component in this energy range is mainly dominated by the exponential cutoff.

##  8 Conclusions

In this paper, we suggest that the ‘second Galactic component of CRs’, needed
to explain the observed flux of CRs in the range between the ‘knee’ and the
‘ankle’ ($10^{7}$ GeV to $10^{9}$ GeV ), may arise from a distribution of
massive star clusters.

This component can bridge the gap between the SNR-CR component, which
dominates below $\sim 10^{7}$ GeV, and the extragalactic component, which
dominates above $\sim 10^{9}$ GeV. It has been previously noted that SNR-CRs
and CRs from star clusters need not be considered two separate components, but
rather originating from similar sources, viz. massive star clusters, the less
massive ones leading to individual SNRs and SNR-CRs, while the more massive
ones can accelerate CRs in a variety of strong shocks appearing in the dense
cluster environment. We have argued that there is a possibility of
acceleration of protons up to a few tens of PeV by considering the particle
acceleration around the WTS, as well as acceleration by SNR shocks inside
massive star clusters. This value is larger than that possible in the standard
paradigm of CR acceleration inside supernova remnants present in the ISM. In
this paper, we have carried out a detailed calculation of the propagation of
cosmic rays in the Galaxy and demonstrated that this model can possibly
explain the all-particle CR spectrum measured at the Earth. Our calculation
considers a realistic distribution of star clusters in the Galaxy and also
includes all the important transport processes of CRs including re-
acceleration by the old SNRs in the Galaxy.

Our analysis requires a proton cut-off energy of $\sim 5\times 10^{7}$ GeV
($50$ PeV) for the CRs accelerated in star clusters. A comparison of our
analytical results with the observed all-particle CR spectrum yields an
injection fraction (the fraction of kinetic energy of shocks being deposited
in CRs) of $\sim 5\%$ (depending on the choice of the extragalactic
component). Furthermore, the variation of the mean logarithmic mass with CR
energy (especially in the energy range of around $10^{7}\hbox{--}10^{9}$ GeV)
supports the argument that the suggested CR component from star clusters can
be considered as the second Galactic component of CRs.

## Acknowledgements

We would like to thank Arpita Roy for the stellar wind elemental abundance
data and acknowledge the Australian supercomputer facilities GADI and AVATAR,
where those simulations have been run. We thank Manami Roy, and Alankar Dutta
for the valuable discussions. We also thank Thibault Vieu for the helpful
discussions. SB acknowledges the Prime Minister’s Research Fellowship (PMRF)
and Govt. of India for financial support. ST acknowledges funding from the Abu
Dhabi Award for Research Excellence (AARE19-224) and the Khalifa University
Emerging Science & Innovation Grant (ESIG-2023-008). PS acknowledges a
Swarnajayanti Fellowship (DST/SJF/PSA-03/2016-17) and a National
Supercomputing Mission (NSM) grant from the Department of Science and
Technology, India.

## Data Availability

The data not explicitly presented in the paper will be available upon
reasonable request from the first author.

## References

  * Aartsen et al. (2013) Aartsen, M. G., Abbasi, R., Abdou, Y., et al. 2013, Phys. Rev. D, 88, 042004, doi: 10.1103/PhysRevD.88.042004
  * Abeysekara et al. (2021) Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2021, Nature Astronomy, 5, 465, doi: 10.1038/s41550-021-01318-y
  * Ackermann et al. (2013) Ackermann, M., Ajello, M., Allafort, A., et al. 2013, Science, 339, 807, doi: 10.1126/science.1231160
  * Adriani et al. (2011) Adriani, O., Barbarino, G. C., Bazilevskaya, G. A., et al. 2011, Science, 332, 69, doi: 10.1126/science.1199172
  * Aguilar et al. (2015) Aguilar, M., Aisa, D., Alpat, B., et al. 2015, Phys. Rev. Lett., 115, 211101, doi: 10.1103/PhysRevLett.115.211101
  * Aharonian et al. (2019) Aharonian, F., Yang, R., & de Oña Wilhelmi, E. 2019, Nature Astronomy, 3, 561, doi: 10.1038/s41550-019-0724-0
  * Ahn et al. (2009) Ahn, H. S., Allison, P., Bagliesi, M. G., et al. 2009, ApJ, 707, 593, doi: 10.1088/0004-637X/707/1/593
  * Aloisio et al. (2014) Aloisio, R., Berezinsky, V., & Blasi, P. 2014, J. Cosmology Astropart. Phys, 2014, 020, doi: 10.1088/1475-7516/2014/10/020
  * Aloisio et al. (2012) Aloisio, R., Berezinsky, V., & Gazizov, A. 2012, Astroparticle Physics, 39, 129, doi: 10.1016/j.astropartphys.2012.09.007
  * Amenomori et al. (2008) Amenomori, M., Bi, X. J., Chen, D., et al. 2008, ApJ, 678, 1165, doi: 10.1086/529514
  * An et al. (2019) An, Q., Asfandiyarov, R., Azzarello, P., et al. 2019, Science Advances, 5, eaax3793, doi: 10.1126/sciadv.aax3793
  * Antoni et al. (2005) Antoni, T., Apel, W. D., Badea, A. F., et al. 2005, Astroparticle Physics, 24, 1, doi: 10.1016/j.astropartphys.2005.04.001
  * Axford (1994) Axford, W. I. 1994, ApJS, 90, 937, doi: 10.1086/191928
  * Axford et al. (1977) Axford, W. I., Leer, E., & Skadron, G. 1977, in International Cosmic Ray Conference, Vol. 11, International Cosmic Ray Conference, 132
  * Bagnulo et al. (2020) Bagnulo, S., Wade, G. A., Nazé, Y., et al. 2020, A&A, 635, A163, doi: 10.1051/0004-6361/201937098
  * Bell (1978) Bell, A. R. 1978, MNRAS, 182, 147, doi: 10.1093/mnras/182.2.147
  * Bell (2013) —. 2013, Astroparticle Physics, 43, 56, doi: 10.1016/j.astropartphys.2012.05.022
  * Bell & Lucek (2001) Bell, A. R., & Lucek, S. G. 2001, MNRAS, 321, 433, doi: 10.1046/j.1365-8711.2001.04063.x
  * Berezhnev (2015) Berezhnev, S. F. 2015, Proc. 33rd ICRC, Paper ID 326
  * Berezinskii & Grigor’eva (1988) Berezinskii, V. S., & Grigor’eva, S. I. 1988, A&A, 199, 1
  * Berezinsky et al. (2006) Berezinsky, V., Gazizov, A., & Grigorieva, S. 2006, Phys. Rev. D, 74, 043005, doi: 10.1103/PhysRevD.74.043005
  * Bhadra et al. (2022) Bhadra, S., Gupta, S., Nath, B. B., & Sharma, P. 2022, MNRAS, 510, 5579, doi: 10.1093/mnras/stac023
  * Biermann & Cassinelli (1993) Biermann, P. L., & Cassinelli, J. P. 1993, A&A, 277, 691, doi: 10.48550/arXiv.astro-ph/9305003
  * Binns et al. (2008) Binns, W. R., Wiedenbeck, M. E., Arnould, M., et al. 2008, New A Rev., 52, 427, doi: 10.1016/j.newar.2008.05.008
  * Blandford & Ostriker (1978) Blandford, R. D., & Ostriker, J. P. 1978, ApJ, 221, L29, doi: 10.1086/182658
  * Blasi (2013) Blasi, P. 2013, A&A Rev., 21, 70, doi: 10.1007/s00159-013-0070-7
  * Blasi & Amato (2012) Blasi, P., & Amato, E. 2012, J. Cosmology Astropart. Phys, 2012, 011, doi: 10.1088/1475-7516/2012/01/011
  * Blazère et al. (2015) Blazère, A., Neiner, C., Tkachenko, A., Bouret, J. C., & Rivinius, T. 2015, A&A, 582, A110, doi: 10.1051/0004-6361/201526855
  * Bloemen et al. (1993) Bloemen, J. B. G. M., Dogiel, V. A., Dorman, V. L., & Ptuskin, V. S. 1993, A&A, 267, 372
  * Breitschwerdt et al. (2002) Breitschwerdt, D., Dogiel, V. A., & Völk, H. J. 2002, A&A, 385, 216, doi: 10.1051/0004-6361:20020152
  * Bronfman et al. (2000) Bronfman, L., Casassus, S., May, J., & Nyman, L. Å. 2000, A&A, 358, 521, doi: 10.48550/arXiv.astro-ph/0006104
  * Bustard et al. (2017) Bustard, C., Zweibel, E. G., & Cotter, C. 2017, ApJ, 835, 72, doi: 10.3847/1538-4357/835/1/72
  * Bykov (2014) Bykov, A. M. 2014, A&A Rev., 22, 77, doi: 10.1007/s00159-014-0077-8
  * Bykov & Fleishman (1992) Bykov, A. M., & Fleishman, G. D. 1992, MNRAS, 255, 269, doi: 10.1093/mnras/255.2.269
  * Bykov et al. (2020) Bykov, A. M., Marcowith, A., Amato, E., et al. 2020, Space Sci. Rev., 216, 42, doi: 10.1007/s11214-020-00663-0
  * Bykov & Toptygin (2001) Bykov, A. M., & Toptygin, I. N. 2001, Astronomy Letters, 27, 625, doi: 10.1134/1.1404456
  * Cao et al. (2021) Cao, Z., Aharonian, F. A., An, Q., et al. 2021, Nature, 594, 33, doi: 10.1038/s41586-021-03498-z
  * Caprioli (2015) Caprioli, D. 2015, in International Cosmic Ray Conference, Vol. 34, 34th International Cosmic Ray Conference (ICRC2015), 8, doi: 10.22323/1.236.0008
  * Caprioli et al. (2011) Caprioli, D., Blasi, P., & Amato, E. 2011, Astroparticle Physics, 34, 447, doi: 10.1016/j.astropartphys.2010.10.011
  * Cesarsky & Montmerle (1983) Cesarsky, C. J., & Montmerle, T. 1983, Space Sci. Rev., 36, 173, doi: 10.1007/BF00167503
  * di Matteo (2015) di Matteo, A. 2015, for the Pierre Auger Collaboration 2015, POS(ICRC2015)249, doi: doi.org/10.22323/1.236.0249
  * Drury et al. (1994) Drury, L. O., Aharonian, F. A., & Voelk, H. J. 1994, A&A, 287, 959, doi: 10.48550/arXiv.astro-ph/9305037
  * Engelmann et al. (1990) Engelmann, J. J., Ferrando, P., Soutoul, A., et al. 1990, A&A, 233, 96
  * Fermi (1949) Fermi, E. 1949, Physical Review, 75, 1169, doi: 10.1103/PhysRev.75.1169
  * Ferrand & Marcowith (2010) Ferrand, G., & Marcowith, A. 2010, A&A, 510, A101, doi: 10.1051/0004-6361/200913520
  * Gabici et al. (2019) Gabici, S., Evoli, C., Gaggero, D., et al. 2019, International Journal of Modern Physics D, 28, 1930022, doi: 10.1142/S0218271819300222
  * Greisen (1966) Greisen, K. 1966, Phys. Rev. Lett., 16, 748, doi: 10.1103/PhysRevLett.16.748
  * Gupta et al. (2018) Gupta, S., Nath, B. B., & Sharma, P. 2018, MNRAS, 479, 5220, doi: 10.1093/mnras/sty1846
  * Gupta et al. (2020) Gupta, S., Nath, B. B., Sharma, P., & Eichler, D. 2020, MNRAS, 493, 3159, doi: 10.1093/mnras/staa286
  * H. E. S. S. Collaboration et al. (2018) H. E. S. S. Collaboration, Abdalla, H., Abramowski, A., et al. 2018, A&A, 612, A6, doi: 10.1051/0004-6361/201629790
  * Heger et al. (2000) Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368, doi: 10.1086/308158
  * Heger et al. (2005) Heger, A., Woosley, S. E., & Spruit, H. C. 2005, ApJ, 626, 350, doi: 10.1086/429868
  * Higdon & Lingenfelter (2013) Higdon, J. C., & Lingenfelter, R. E. 2013, ApJ, 775, 110, doi: 10.1088/0004-637X/775/2/110
  * Higdon et al. (1998) Higdon, J. C., Lingenfelter, R. E., & Ramaty, R. 1998, ApJ, 509, L33, doi: 10.1086/311757
  * High Resolution Fly’S Eye Collaboration et al. (2009) High Resolution Fly’S Eye Collaboration, Abbasi, R. U., Abu-Zayyad, T., et al. 2009, Astroparticle Physics, 32, 53, doi: 10.1016/j.astropartphys.2009.06.001
  * Hillas (1967) Hillas, A. M. 1967, Physics Letters A, 24, 677, doi: 10.1016/0375-9601(67)91023-7
  * Hörandel (2003) Hörandel, J. R. 2003, Journal of Physics G Nuclear Physics, 29, 2439, doi: 10.1088/0954-3899/29/11/001
  * Jokipii & Morfill (1987) Jokipii, J. R., & Morfill, G. 1987, ApJ, 312, 170, doi: 10.1086/164857
  * Jones et al. (2001) Jones, F. C., Lukasiak, A., Ptuskin, V., & Webber, W. 2001, ApJ, 547, 264, doi: 10.1086/318358
  * Kääpä et al. (2023) Kääpä, A., Kampert, K.-H., & Becker Tjus, J. 2023, in European Physical Journal Web of Conferences, Vol. 283, European Physical Journal Web of Conferences, 03006, doi: 10.1051/epjconf/202328303006
  * Kampert & Unger (2012) Kampert, K.-H., & Unger, M. 2012, Astroparticle Physics, 35, 660, doi: 10.1016/j.astropartphys.2012.02.004
  * Kang et al. (1996) Kang, H., Ryu, D., & Jones, T. W. 1996, ApJ, 456, 422, doi: 10.1086/176666
  * Kelner et al. (2006) Kelner, S. R., Aharonian, F. A., & Bugayov, V. V. 2006, Phys. Rev. D, 74, 034018, doi: 10.1103/PhysRevD.74.034018
  * Knödlseder (2013) Knödlseder, J. 2013, in Astrophysics and Space Science Proceedings, Vol. 34, Cosmic Rays in Star-Forming Environments, ed. D. F. Torres & O. Reimer, 169, doi: 10.1007/978-3-642-35410-6_13
  * Knurenko & Sabourov (2011) Knurenko, S. P., & Sabourov, A. 2011, Astrophysics and Space Sciences Transactions, 7, 251, doi: 10.5194/astra-7-251-2011
  * Krymskii (1977) Krymskii, G. F. 1977, Akademiia Nauk SSSR Doklady, 234, 1306
  * Lagage & Cesarsky (1983) Lagage, P. O., & Cesarsky, C. J. 1983, A&A, 125, 249
  * Lerche & Schlickeiser (1982) Lerche, I., & Schlickeiser, R. 1982, A&A, 116, 10
  * Longmore et al. (2014) Longmore, S. N., Kruijssen, J. M. D., Bastian, N., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 291–314, doi: 10.2458/azu_uapress_9780816531240-ch013
  * Mannheim et al. (2000) Mannheim, K., Protheroe, R. J., & Rachen, J. P. 2000, Phys. Rev. D, 63, 023003, doi: 10.1103/PhysRevD.63.023003
  * Mao & Shen (1972) Mao, C. Y., & Shen, C. S. 1972, Chinese Journal of Physics, 10, 16
  * Marcowith et al. (2006) Marcowith, A., Lemoine, M., & Pelletier, G. 2006, A&A, 453, 193, doi: 10.1051/0004-6361:20054738
  * Montmerle (1979) Montmerle, T. 1979, ApJ, 231, 95, doi: 10.1086/157166
  * Morlino (2017) Morlino, G. 2017, in Handbook of Supernovae, ed. A. W. Alsabti & P. Murdin, 1711, doi: 10.1007/978-3-319-21846-5_11
  * Mueller et al. (1991) Mueller, D., Swordy, S. P., Meyer, P., L’Heureux, J., & Grunsfeld, J. M. 1991, ApJ, 374, 356, doi: 10.1086/170125
  * Nath & Eichler (2020) Nath, B. B., & Eichler, D. 2020, MNRAS, 499, L1, doi: 10.1093/mnrasl/slaa151
  * Neiner et al. (2015) Neiner, C., Grunhut, J., Leroy, B., De Becker, M., & Rauw, G. 2015, A&A, 575, A66, doi: 10.1051/0004-6361/201425193
  * Obermeier et al. (2011) Obermeier, A., Ave, M., Boyle, P., et al. 2011, ApJ, 742, 14, doi: 10.1088/0004-637X/742/1/14
  * Oey & Clarke (1997) Oey, M. S., & Clarke, C. J. 1997, MNRAS, 289, 570, doi: 10.1093/mnras/289.3.570
  * Panov et al. (2007) Panov, A. D., Adams, J. H., J., Ahn, H. S., et al. 2007, Bulletin of the Russian Academy of Sciences, Physics, 71, 494, doi: 10.3103/S1062873807040168
  * Rachen (2015) Rachen, J. 2015, Proc. 28th Texas Symp. Relativ. Astrophys., Geneva, Switzerland, 13–18 December 2015, Electronic Proceedings, 230
  * Rachen & Biermann (1993) Rachen, J. P., & Biermann, P. L. 1993, A&A, 272, 161, doi: 10.48550/arXiv.astro-ph/9301010
  * Reville & Bell (2012) Reville, B., & Bell, A. R. 2012, MNRAS, 419, 2433, doi: 10.1111/j.1365-2966.2011.19892.x
  * Roy et al. (2021) Roy, A., Dopita, M. A., Krumholz, M. R., et al. 2021, MNRAS, 502, 4359, doi: 10.1093/mnras/stab376
  * Salpeter (1955) Salpeter, E. E. 1955, ApJ, 121, 161, doi: 10.1086/145971
  * Sarkar et al. (2015) Sarkar, K. C., Nath, B. B., Sharma, P., & Shchekinov, Y. 2015, MNRAS, 448, 328, doi: 10.1093/mnras/stu2760
  * Schure & Bell (2013) Schure, K. M., & Bell, A. R. 2013, MNRAS, 435, 1174, doi: 10.1093/mnras/stt1371
  * Strong & Moskalenko (1998) Strong, A. W., & Moskalenko, I. V. 1998, ApJ, 509, 212, doi: 10.1086/306470
  * Suzuki et al. (2022) Suzuki, H., Bamba, A., Yamazaki, R., & Ohira, Y. 2022, ApJ, 924, 45, doi: 10.3847/1538-4357/ac33b5
  * Taillet & Maurin (2003) Taillet, R., & Maurin, D. 2003, A&A, 402, 971, doi: 10.1051/0004-6361:20030318
  * Tatischeff et al. (2021) Tatischeff, V., Raymond, J. C., Duprat, J., Gabici, S., & Recchia, S. 2021, MNRAS, 508, 1321, doi: 10.1093/mnras/stab2533
  * The Pierre Auger Collaboration et al. (2013) The Pierre Auger Collaboration, Aab, A., Abreu, P., et al. 2013, arXiv e-prints, arXiv:1307.5059, doi: 10.48550/arXiv.1307.5059
  * The Pierre Auger Collaboration et al. (2015) —. 2015, arXiv e-prints, arXiv:1509.03732, doi: 10.48550/arXiv.1509.03732
  * Thoudam & Hörandel (2012) Thoudam, S., & Hörandel, J. R. 2012, MNRAS, 421, 1209, doi: 10.1111/j.1365-2966.2011.20385.x
  * Thoudam & Hörandel (2014) —. 2014, A&A, 567, A33, doi: 10.1051/0004-6361/201322996
  * Thoudam et al. (2016) Thoudam, S., Rachen, J. P., van Vliet, A., et al. 2016, A&A, 595, A33, doi: 10.1051/0004-6361/201628894
  * Tibet AS$\gamma$ Collaboration et al. (2021) Tibet AS$\gamma$ Collaboration, Amenomori, M., Bao, Y. W., et al. 2021, Nature Astronomy, 5, 460, doi: 10.1038/s41550-020-01294-9
  * Unger et al. (2015) Unger, M., Farrar, G. R., & Anchordoqui, L. A. 2015, Phys. Rev. D, 92, 123001, doi: 10.1103/PhysRevD.92.123001
  * Vieu & Reville (2023) Vieu, T., & Reville, B. 2023, MNRAS, 519, 136, doi: 10.1093/mnras/stac3469
  * Vieu et al. (2022) Vieu, T., Reville, B., & Aharonian, F. 2022, MNRAS, 515, 2256, doi: 10.1093/mnras/stac1901
  * Voelk & Biermann (1988) Voelk, H. J., & Biermann, P. L. 1988, ApJ, 333, L65, doi: 10.1086/185289
  * Wandel et al. (1987) Wandel, A., Eichler, D., Letaw, J. R., Silberberg, R., & Tsao, C. H. 1987, ApJ, 316, 676, doi: 10.1086/165233
  * Waxman (1995) Waxman, E. 1995, Phys. Rev. Lett., 75, 386, doi: 10.1103/PhysRevLett.75.386
  * Weaver et al. (1977) Weaver, R., McCray, R., Castor, J., Shapiro, P., & Moore, R. 1977, ApJ, 218, 377, doi: 10.1086/155692
  * Webb et al. (1985) Webb, G. M., Axford, W. I., & Forman, M. A. 1985, ApJ, 298, 684, doi: 10.1086/163652
  * Werner et al. (2015) Werner, M., Kissmann, R., Strong, A. W., & Reimer, O. 2015, Astroparticle Physics, 64, 18, doi: 10.1016/j.astropartphys.2014.10.005
  * Wiedenbeck et al. (1999) Wiedenbeck, M. E., Binns, W. R., Christian, E. R., et al. 1999, ApJ, 523, L61, doi: 10.1086/312242
  * Yoon et al. (2011) Yoon, Y. S., Ahn, H. S., Allison, P. S., et al. 2011, ApJ, 728, 122, doi: 10.1088/0004-637X/728/2/122
  * Zatsepin & Kuz’min (1966) Zatsepin, G. T., & Kuz’min, V. A. 1966, Soviet Journal of Experimental and Theoretical Physics Letters, 4, 78
  * Zirakashvili & Völk (2006) Zirakashvili, V. N., & Völk, H. J. 2006, Advances in Space Research, 37, 1923, doi: 10.1016/j.asr.2005.06.013