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

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 Ne²²/Ne²⁰ 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 (²²Ne/²⁰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 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.

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.

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}$ cm²s^-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.

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.

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

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,

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.

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.

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

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.

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.

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.

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