Between the Cosmic-Ray “Knee” and the “Ankle”: Contribution from Star Clusters

We show that massive, young star clusters may be possible candidates that can accelerate Galactic cosmic rays (CRs) in the range of 107–109 GeV (between the “knee” and “ankle”). Various plausible scenarios, such as acceleration at the wind termination shock and supernova shocks inside these young star clusters, have been proposed, since it is difficult to accelerate particles up to the 107–109 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 reacceleration 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 × 107 Z GeV (50 Z PeV) from these clusters, together with a CR injection fraction of ∼5% of the wind kinetic energy. We discuss the possibility of achieving these requirements in star clusters, as well as the associated uncertainties, in the context of considering star clusters as the natural accelerator of the “second component” of Galactic CRs.


INTRODUCTION
Cosmic rays (CRs hereafter) are high-energy particles that span an extensive range of energy from 1 GeV to ∼ 10 11 GeV.Lower energy CRs up to ∼ 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 ∼ 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 ∼ 10 5 − 10 6 years, most of the particle acceleration occurs during the un-decelerated blast wave phase, which lasts for ≤ 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 ∼ µG, Lagage & Cesarsky (1983) estimated the upper limit of CR energy in supernova remnants to be ∼ 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 γ-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γ 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 ∼ 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 ∼ 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 ∼ 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 × 10 8 GeV, the surface magnetic field of a Wolf-Rayet star is required to be ≈ 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 ∼ 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 ∼ 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 ≲ 10s 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 Ne 22 /Ne 20 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 γ-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 starforming regions have been discussed as potential candidates for CR accelerators (e.g.Bykov 2014 ); these γ-ray observations have strengthened the hypothesis of CR acceleration in these environments.Recently, Gupta et al. (2020) has shown that the excess ( 22 Ne/ 20 Ne) 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 bub-ble 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 ∼ 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 γ-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 ∼ 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 ∼ 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.

First Galactic component: SNR-CRs
As mentioned in the Introduction, supernova remnants are the most likely candidate for cosmic-ray acceleration up to ∼ 10 6 GeV energy (Lagage & Cesarsky 1983).The diffusive shock acceleration at strong shocks produces a power-law spectrum with an index of ∼ −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 × 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 × E c = 6.5 × 10 7 GeV.This result shows that SNR-CRs can contribute only ∼ 30% of the total observed intensity above ∼ 2 × 10 7 GeV (Thoudam et al. 2016).Therefore, additional components are required to explain the all-particle spectrum in the ≳ 10 7 GeV range.

Extragalactic component
Various previous works have already pointed out that the 'ankle'-like feature of the CR spectrum at ≳ 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. 2012Aloisio et al. , 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 photodisintegration 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 ultrahigh 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 ∼ 3 × 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 ab- sorption (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.

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 × 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 µG in the Galaxy, cosmic rays with energy Z ×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 cutoff 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 ≤ 10 Myr can contribute significantly to the observed total flux of CRs (Gupta et al. 2020).Recently LHAASO has observed γ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.

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 where r is the Galactocentric distance and σ = 3 kpc and Σ 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πr dr), as where, Σ 0 ∼ 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-α 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 ∼ 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.

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 Here we include spatial diffusion (first term on the lefthand side), re-acceleration (terms with coefficient ζ), and interaction losses (∝ σ, 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, σ(p) is the cross-section of inelastic collision, N is the differential number density (number per unit volume per momentum) and ζ 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 powerlaw distribution with an index of s ∼ 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 = ±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 ∼ 5 kpc (Taillet & Maurin 2003).
The term on the right side, Q(r, p)δ(z), represents the injection rate of cosmic rays per unit volume in the momentum bin [p, p+dp] by the sources.The δ(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 momentumdependent part, i.e., where ν(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 lowmomentum 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 ∼ 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 ∼ 4.6 kpc (0.55R 0 , where R 0 is the distance of the Earth from the Galactic center ∼ 8.5 kpc) and then decreases rapidly at large distances (for details see section 3.1).
The expression for surface density of star clusters ν has been calculated in section 3.1, and the power-law source spectrum is described in section 3.3.The energydependent diffusion coefficient as a function of particle rigidity follows where D 0 is the diffusion constant, ρ = Apc/Ze is the particle rigidity, β = v(p)/c where v(p) is the CR particle velocity and c is the speed of light, δ = 0.33 is the diffusion index, and ρ 0 = 3 GV is a normalisation constant.
In this injection-diffusion-reacceleration1 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 ζ can be expressed as, ζ = ηV ν SN , where V = 4πR 3 /3 is the volume occupied by each SNR of radius R re-accelerating the cosmic rays.Here, η is a correction factor that takes care of the actual unknown size of the remnants, and ν SN is the rate of supernova explosions per unit surface area in the Galactic disk.The values of R and ν SN have been taken as 100 pc and 25 SNe Myr −1 kpc −2 respectively (Thoudam & Hörandel 2014).
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 ′ , z, p, p ′ ) (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 Substituting the obtained G(r, r ′ , z, p, p ′ ) (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 where R = 20 kpc is the radial boundary of the Galaxy, Σ 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), Equation 4 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 4 are discussed in sections 3.2-3.5.

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, 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 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.

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.

Acceleration at wind termination shock (WTS):
In equation 7, p 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 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, κ = pc 2 /(ζqB), the maximum energy is then (Vieu et al. 2022): Here R WTS is the radius of WTS.In the above equation, V w is the velocity of stellar wind, B WTS is the value of the magnetic field at the WTS position, ζ = 3r g /λ, with λ the mean free path due to the magnetic field.The Bohm diffusion, which is the most optimistic scenario, corresponds to the limit ζ = 3.
We follow the arguments advocated by Vieu et al. (2022) to estimate the magnetic field in the cluster core, Here, n c is the core density, η T is the efficiency of generation of turbulence, N 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 WTS R WTS .Therefore, equation 8 can be expressed as, This leads to a maximum estimate: Equation 11 gives a conservative estimate of E max = 6 PeV (6 × 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 η T , wind velocity v w ) that can conceivably increase the maximum energy by a factor of a few.

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: PeV . ( For a typical young cluster R WTS /R c ∼ 5 − 30, which gives 1 − (R c /R WTS ) 1/7 ∼ 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 ≥ 2×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.
The recent detection of γ-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 γ-ray energy can be approximated as E cr ≈ 10E γ .Some of those sources possibly are young massive clusters (see extended table 2 of (Cao et al. 2021).The γ-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.

Elemental abundances in star cluster winds
We consider a simple stellar population formed at time t = 0 with an initial mass function (IMF) dn dm ∝ 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, is the cumulative mass of element X ejected in winds by a star of initial mass m up to age t where, We use the mass loss rate for each nucleus ṁ(X, m, t ′ ) 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, 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, Using this method, we have calculated the massweighted 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.

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 ϕ(L) ∝ L −2 .We use this distribution to calculate the average luminosity of clusters with kinetic luminosity in the range L min = 10 37 erg s −1 (corresponds to N OB = 10) to L max = 10 39 erg s −1 (corresponds to N OB = 1000) as following, 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 ⟨L w ⟩ on L max is weak, but there is the sensitive dependence on L min , the implications of which we discuss later.

MODEL PREDICTION FOR THE SECOND COMPONENT OF GALACTIC COSMIC RAYS
The values of cosmic ray propagation parameters (D 0 , δ; the normalization of the diffusion coefficient and its power-law dependence on momentum) and reacceleration parameters (η , 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 × 10 28 cm 2 s −1 , η = 1.02, s = 4.5, and δ = 0.33 (Thoudam & Hörandel 2014).We have also used these values in our model.For the 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 .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 × 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 (σ(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 Se-quence (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.
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×10 7 GeV (50 PeV) and the injection fraction of ∼ 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)  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 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.

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 ∼ 2.5 × 10 6 GeV.For the extragalactic component, we have adopted the UFA model (Unger et al. 2015), which considers a significant contri-bution of extragalactic CRs below the ankle to reproduce the observed CR energy spectrum as well as X max (the depth of the shower maximum) and the variance of X 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 × 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) Several ground-based experiments such as KAS-CADE, TUNKA, LOFAR, and the Pierre Auger Obser-vatory have provided measurements of the composition of CRs at energies above ∼ 10 6 GeV.Heavier nuclei interact at a higher altitude in the atmosphere, which results in smaller values of X max as compared to lighter nuclei.For comparison with the theoretical predictions, ⟨ln A⟩, the mean logarithmic mass of the measured cosmic rays, is of utmost importance.This can be obtained from the measured X max values using the following relation mentioned in Hörandel (2003), Here X p max and X Fe max represent the average maximum depths of the shower for protons and iron nuclei, respectively, and 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, 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 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 × 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 ∼ 5% and an energy cutoff of 5 × 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.

VARYING THE EXTRAGALACTIC COMPONENT
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.

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 ∼ 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 × 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 γ-rays from a few SNRs have also identified a few SNRs as cosmic ray PeVatrons that can accelerate particles up to a few times ∼ 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 × 10 7 GeV energy.For other elements with atomic number Z, the maximum energy is 5×10 7 Z GeV in these young compact star clusters, and a cosmic ray injection fraction of ∼ 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 max = 50 PeV (5 × 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 max , as required to fit the observed spectrum.Also note that the recently detected γ-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 ∆ can be calculated from our model using the diffusion approximation given by Mao & Shen (1972), where N is the CR number density, D is the diffusion coefficient and c is the velocity of light.From our model, we get ∆ ∼ 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) × 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 ∼ 50 pc size, and the injected mechanical power can reach ∼ 10 38 erg s −1 over 10 Myr-the lifetime of a superbubble.This process is supplemented by the formation of multiple shocks, largescale 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.

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 ∼ 10 −3 × 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, ∼ 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(Heger et al. , 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 ⟨ln A⟩.
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 µG in the Galactic plane then the transition from diffusion to drift will occur at ∼ Z×10 17−18 eV (Kääpä et al. 2023).The maximum energy for proton in WTS is ∼ 50 PeV (5 × 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.

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 ∼ 10 7 GeV, and the extragalactic component, which dominates above ∼ 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 ∼ 5 × 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 ∼ 5% (depending on the choice of the extragalactic component).Furthermore, the varia-tion of the mean logarithmic mass with CR energy (especially in the energy range of around 10 7 -10 9 GeV) supports the argument that the suggested CR component from star clusters can be considered as the second Galactic component of CRs.

Figure 1 .
Figure 1.Schematic diagram of a stellar wind bubble.The position of termination shock is Rs; R cd and R fs are contact discontinuity and forward shock positions, respectively.

Figure 2 .
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 (Σ; see Eq. 2) as a function of distance from the Galactic center.

Figure 3 .
Figure 3. Model prediction for the star cluster model as a second galactic component considering an injection fraction ∼ 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 Ec = 5 × 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.

Figure 4 .
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 ∼ 5%, and the cutoff is at 5 × 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π)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.

Figure 5 .
Figure 5. Mean logarithmic mass ⟨ln A⟩ 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 Xmax values to ⟨ln A⟩ (see equation 18).

Figure 6 .
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.

Table 1 .
Roy et al. (2021)he SNR shocks inside a Source spectral indices q and fractional abundances of different elements in the wind material.The elemental abundances are calculated followingRoy et al. (2021).