Massive star clusters as the an alternative source population of galactic cosmic rays

Extended γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document}-ray emissions in the vicinity of young star clusters are believed to be produced by the interaction of CRs accelerated therein with the ambient gas. The detailed spatial analysis reveals 1/r type distribution of CRs, which indicates a continuous injection of CRs in these objects. The hard, ∝E-2.3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\propto E^{-2.3}$$\end{document} type power-law energy spectra of parent protons continue up to ∼\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim $$\end{document} 1 PeV. The efficiency of conversion of kinetic energy of powerful stellar winds to CRs can be as high as 10%. This implies that the young massive stars can operate as effective proton PeVatrons with a major or even dominant contribution to the flux of highest-energy galactic CRs.


Introduction
The origin of cosmic rays (CRs) remains as a mystery after these relativistic particles were discovered for more than one century. The current paradigm postulates that the bulk of the CRs are accelerated by supernova remnants (SNRs). In this regard, the -ray observations provide us indirect evidence on the CRs acceleration in SNRs. The detection of the so-called 0 -decay bump in the spectra of several midage SNRs (F.L. Collaboration 2013; Giuliani et al. 2011) is considered as a substantial evidence of acceleration of protons and nuclei in SNRs. Furthermore, the detection of more than a dozen of young (a few thousand years old or younger) SNRs in TeV -rays (H. Collaboration 2006;H. Collaboration, H.E.S.S. 2007) highlights these objects as efficient particle accelerators. However, the mid-age SNRs reveal a break of CR spectra at about 100 GeV; while, the very origin of the -rays in young SNRs (hadronic or leptonic) is still not clear.
On the other hand, the recent measurements of 60 Fe abundance in CRs (Binns et al. 2016) indicate that a substantial fraction of Cosmic rays (CRs) could be accelerated in young OB star clusters and related super bubbles. The measurements of the Galactic diffuse -ray emissions also show that the CRs have a similar radial distribution to OB stars rather than SNRs (Acero et al. 2016;Yang et al. 2016). The interacting winds of massive stars have been recognised as potential CR accelerators as early as in the 1980s. The acceleration could take place in the vicinity of the stars (Cesarsky and Montmerle 1983) or in superbubbles, multi-parsec structures caused by the collective activity of massive stars (Bykov and Toptygin 1982;Parizot et al. 2004). The acceleration of multiple shocks can raise the maximum energy of CR protons out of 1PeV (Klepach et al. 2000) which makes the stellar clusters the attractive candidates for cosmic PeVatrons. The accelerated CRs with the surrounding gas should produce secondary -ray emissions. The diffuse GeV -rays detected by Fermi LAT around the compact clusters Cygnus OB2 (Ackermann et al. 2011), NGC 3603 (Yang andAharonian 2017) and Westerlund 2 (Yang et al. 2018) can be naturally interpreted within this scenario.
The spatial distribution of CRs is a powerful tool to diagnose the injection history and acceleration site of CRs (Aharonian and Atoyan 1996;HESS Collaboration 2016). This method relies on the accurate measurement of the spatial distribution of both -ray and gas. It has been successfully applied to the diffuse TeV -ray emission of the Central Molecular Zone (CMZ) in the Galactic Centre (GC) (HESS Collaboration 2016). While the hard spectra of -rays extending to energies of tens of TeV, indicate the presence of a proton PeVatron(s) in CMZ, the 1 / r type radial distribution of parents protons revealed up to ∼ 200 pc, points to the continuous operation of proton PeVatron(s) located within the central 10 pc of our Galaxy. We argue that the compact stellar clusters, Arches, Quintuplet and Nuclear in GC, could be alternative sites for the CR acceleration. In this paper, we also explored the possibility of extraction of spatial distributions of CRs in the proximity of other two clusters, Cygnus OB2 and Westerlund 1.

-ray observations
The extended GeV -ray source around the cluster NGC 3603 (Yang and Aharonian 2017) and the TeV -ray source associated with 30 Dor C (H.E.S.S. Collaboration 2015) are too weak for derivation of statistically significant radial distributions of CRs. On the other hand, the angular size of the diffuse GeV source associated with Westerlund 2 is too large to be detected with the current atmospheric Cherenkov telescopes. Fortunately, in the case of Cygnus Cocoon discovered by the Fermi LAT collaboration as a bright extended -ray source associated with Cyg OB2 (Ackermann et al. 2011), the photon statistics is sufficient for derivation of spectral and spatial distributions of CRs. It is also important that -ray emission of this source extends to TeV energies (ARGO-YBJ Collaboration 2014). The same is true also for the extended TeV -ray emitter HESS J1646-458 apparently linked to Westerlund 1 (H. Collaboration 2012). For Cygnus Cocoon, we analysed Fermi LAT data using the standard LAT software package. For HESS J1646-458, we used the radial profiles published by the H.E.S.S collaboration (H. Collaboration 2012). For the distribution of molecular hydrogen, we applied the data from the CO galactic survey performed by the CfA 1.2m millimetre-wave Telescope; while for the atomic hydrogen, we used the data from the Leiden/Argentine/Bonn (LAB) Survey. We also use the results for Westerlund 2 from (Yang et al. 2018) in this work.
The main conclusion is that the CR density declines as r −1 up to ≈ 50 pc from both stellar clusters. The results are shown in Fig. 1b, together with the earlier published radial distributions of CR protons in CMZ (HESS Collaboration 2016). In Fig. 1a, we show the differential -ray luminosities of the extended sources associated with Cyg OB2, Westerlund 1 and CMZ. The energy distributions of -rays are quite similar; dN∕dE ∝ E − type differential energy spectra with power-law index ≈ 2.2 extend to 10 TeV and beyond without an indication of a cutoff or a break. The -rays are likely to originate from interactions of CRs with the ambient gas through the production and decay of neutral -mesons (see below). Because of the increase of the 0 -meson production cross-section with energy of incident protons and nuclei, the spectrum of secondary -rays appears slightly harder compared to the spectrum of parent protons, ≈ p − 0.1 (Kelner et al. 2006), thus the power-law index of the proton distribution should be p ≈ 2.3.

CR radial distribution
The apparent similarity of the radial ( ∝ r −1 ) and energy ( ∝ E −2.3 ) distributions of CR protons for different stellar clusters is a hint that we observe the same phenomenon. The most natural explanation of the 1 / r dependence of CR radial distribution is that relativistic particles have been continuously injected into and diffused away in the interstellar medium (ISM). The characteristic time scale is determined by the time of propagation of CRs over the typical distances of tens of parsecs, t ∼ R 2 ∕D . Formally, in the case of a large diffusion coefficient D (i.e., fast diffusion), it could be as short as 10 3 year. However, this would imply very low efficiency of conversion of CRs to -rays. Given the tight energy budget, the diffusion time cannot be much shorter than the age of the stellar cluster (see below), t ≥ 10 6 year. On the other hand, the acceleration of multi-TeV CRs in an individual SNR cannot last more than 10 4 year (see, e.g., Bell et al. 2013). Thus, to support the quasi-continuous CR injection, an unrealistically high rate of ∼ 1 SN per 100 year in the cluster is required. This disfavors SNRs and gives preference to massive stellar winds as particle accelerators.
In the case of spherically symmetric diffusion, the CR density at a distance from the central source r depends on the injection rate Q (E) and the diffusion coefficient: w(E, r) ∝Q(E)∕rD(E) , i.e., the 1 / r profile is independent of the absolute value of the diffusion coefficient unless the latter varies dramatically over the scales of tens of parsecs.
The above relations are valid when the energy losses of CRs can be neglected. While for CR protons and nuclei this is a fully justified assumption, relativistic electrons undergo severe energy losses. However, electrons cannot be responsible for the observed -ray images. The leptonic (inverse Compton; IC) origin of -rays is excluded both at GeV and TeV energies. Firstly, the propagation of multi-TeV ( ≫ 10 TeV) electrons in the ISM could hardly exceed 100 pc (Atoyan et al. 1995). Moreover, inside a typical cluster of a radius less than 3 pc and the overall starlight luminosity of L r ≈ 10 40 erg/s , the energy density of optical photons exceeds u r ∼ L∕4 r 2 c ≈ 100 eV/cm 3 . Outside the cluster, u r decreases as 1∕r 2 , thus, up to tens of parsecs, it dominates over the average radiation density in the Galactic plane. Therefore, in the case of IC origin of -rays, we would expect a sharp increase of the -ray intensity towards a bright central source coinciding with the cluster. The brightness distributions of the observed -ray images of objects discussed in this paper do not agree with this prediction. It is convenient to write the radial distribution of the CR density in the form Below, we will adopt r 0 = 10 pc, i.e., normalise the CR proton density w 0 outside but not far from the cluster. Then, the (1) w(r) = w 0 (r∕r 0 ) −1 .
total energy of CR protons within the volume of the radius R 0 is F r o m w 0 l i s t e d i n Ta b l e 1 , w e o b t a i n W p ≈ 1.5 × 10 51 , 3.2 × 10 49 , 1.4 × 10 48 , 2.3 × 10 49 e r g s for Westerlund 2 Cocoon, Westerlund 1 Cocoon, Cygnus Cocoon, and CMZ, respectively. This estimate strongly depends on the value of R obs which is determined by the brightness of the -ray image. The extensions of the large diffuse structure depend on the detector's performance, the level of the background, etc. Thus, the content of CR protons within R obs does not provide information about all CRs injected into ISM. The latter can be calculated by integrating The CR proton radial distributions in Cyg Cocoon and Wd 2 Cocoon above 100 GeV, and in Wd 1 Cocoon and CMZ above 10 TeV. The -ray flux enhancement factor due to the contribution of CR nuclei was assumed = 1.5 . For comparison, the energy densities of CR protons above 100 GeV and 10 TeV based on the measurements by AMS are also shown Aguilar et al. (2015) Eq.
(1) up to the so-called diffusion radius R D , the maximum distance penetrated by a particle of energy E during the time T 0 . In the case of negligible energy losses of propagating particles, where D 30 is the diffusion coefficient of protons in the units of 10 30 cm 2 ∕s , and T 6 is T 0 normalised to 10 6 years. The ages of the individual clusters vary in a narrow range between 2 and 7 Myr (see Table 1). In the source neighbourhood, the diffusion coefficient cannot be very large; otherwise, the demand on the total energy in CRs would exceed the available energy contained in the stellar winds, where L 39 = 10 39 L 0 is the total mechanical power of the stellar winds in units of 10 39 erg/s , and f is the efficiency of conversion of the wind kinetic energy to relativistic protons with energy larger than 10 TeV. Substituting R 0 = R D into Eq.
(2), we obtain The large resolved size (300 pc) and the large CR density in Westerlund 2 ( w 0 = 6 eV/cm 3 ), combined with the well-known age ( 2 × 10 6 year) and the available energy budget in the form of kinetic energy of stellar winds ( 2 × 10 38 erg/cm 2 s ), robustly constrain the CR acceleration efficiency in the cluster and the diffusion coefficient in its Cocoon. Indeed, from the obvious condition R obs ≤ R D , Eq.
(3) gives D 29 ≥ 0.04. Substituting this lower limit into Eq. (5), we obtain f ∼ 0.1 , i.e., the acceleration efficiency should be as large as 10 percent. Actually, depending on the shape of the CR spectrum below 100 GeV, the lower limit for f could be by a factor of few higher. This implies that for any reasonable acceleration efficiency, the upper and lower limits on the diffusion coefficient shrink its value to few times 10 27 cm 2 ∕s , significantly smaller than the diffusion coefficient in the interstellar medium. The requirements to the parameters of other clusters are less stringent because of the smaller values of the (resolved) extensions of -ray sources and/or the lower CR densities. However, the reported angular extensions of -ray sources (4) W tot = fL 0 T 0 = 3 × 10 52 fL 39 T 6 erg, could not be used as unbiased measures of the real physical size of the object. The parameter w 0 experimentally derived from the detected images of -rays is the more objective quantity. In particular, the comparable values of w 0 in Westerlund 2 Cocoon and Westerlund 1 Cocoon (taking into account that the densities in these objects are derived in different energy bands) seem quite natural given the almost identical parameters of characterising these two clusters of massive stars. On the other hand, the significantly low level of the CR density in Cygnus Cocoon and CMZ can be explained either by the low efficiency of conversion of the kinetic energy of the winds to CRs, and/or faster diffusion of relativistic particles in these objects. The case of CMZ is especially interesting, given that the total kinetic energy power in three ultracompact (Arches, Quintuplet and Nuclear) clusters is L ≃ 10 39 erg/s , i.e., exceeds by an order of magnitude the overall stellar wind power in Westerlund 1 and Westerlund 2. One can see from Eq. (5) that the acceleration efficiency could significantly, in principle, exceed 1 percent, provided that the diffusion coefficient in CMZ is much larger than in Westerlund 2 Cocoon. Although this cannot be a priori excluded, the alternative assumption regarding the low efficiency of CR acceleration seems a more likely option given the unusual nature of these ultracompact clusters where tens of massive OB stars are packed within the few pc linear size regions. Actually, the acceleration efficiency exceeding 10%, as derived for Westerlund 2, should not be typical for all star clusters. Otherwise, it would lead to overproduction of CRs, given that the overall kinetic energy power of massive stellar winds exceeds 10 42 erg/s.
The spectra of CR protons inside of all three diffuse -ray sources are described by power-law energy distributions with an index p ≈ 2.3 . It is formed from the initial (acceleration) spectrum of protons, Q (E) ∝ E − 0 , but can be modified due to the energy-dependent diffusion, J p (E, r) ∝Q(E)∕D(E)r −1 . For the Kolmogorov type turbulence, D(E) ∝ E 1∕3 , we arrive at a "classical" E −2 type acceleration spectrum. One cannot, however, exclude that at energies ≥ 10 TeV, the diffusion slightly depends on the energy. Even in this extreme case, the acceleration spectrum of protons would be relatively hard with 0 = p ≈ 2.3 . The hard -ray spectra of  (Figer 2008) 1.5-2.5 3-6 2-7 4-6 L kin of cluster (erg/s) 2 × 10 38 (Rauw et al. 2007;Reimer et al. 2008) 2 × 10 38 (Ackermann et al. 2011) 1 × 10 39 (Muno et al. 2006) 1 × 10 39 (Hußmann 2014) Dist ( Westerlund 1 Cocoon and CMZ continue up to 20-30 TeV without an indication of a cutoff or a break. Correspondingly, the energy spectra of parent protons should not break at least until 0.5 PeV (see Fig. 1a). This makes the clusters of massive stars potential sources of multi-TeV neutrinos with a fair chance to be detected by the cubic-km volume neutrino detectors. In particular, Westerlund 1, which has the highest -ray flux at 20 TeV, seems to be a promising target for neutrino observations (Bykov 2014).

Conclusion
The main contributors to the stellar wind are O stars and WR stars due to their considerable mass loss rate and high wind velocity. The power of a single O star is of the order 10 36 − −10 37 erg/s , depending on the mass; while, the power of WR star can be as high as 6 × 10 37 erg/s (Cesarsky and Montmerle 1983; Parizot et al. 2004). The total power of the O stars and WR stars in our Galaxy is estimated to be 3 × 10 41 erg/s using the initial mass function of O stars and the star counts of WR stars in the solar neighbourhood (Cesarsky and Montmerle 1983). The total CR injection power in our Galaxy is estimated to be between 6 × 10 40 erg/s and 3 × 10 41 erg/s . Thus, the stellar wind can provide 10% to 50 % of the CRs in our Galaxy if 10% of the mechanical energy can be converted to CRs. Thus, we argue that the stellar clusters offer a viable solution to the longstanding problem of the origin of Galactic CRs with massive/luminous stars as major contributors to observed fluxes of CRs up to the knee around 1 PeV. Furthermore, the multi-TeV -ray observations provide evidence that the clusters of massive stars operating as PeVatrons may substantially contributing to the flux of galactic CRs. The extension of spectrometric and morphological -ray measurements up to 100 TeV in the energy spectrum and up to several degrees in the angular size from regions surrounding powerful stellar clusters would provide crucial information about the origin of CRs in general, and the physics of proton PeVatrons, in particular. Such observations with the Cherenkov Telescope Array will be available in coming years.