Supernovae in compact star clusters as sources of high-energy cosmic rays and neutrinos

TL;DR

This paper proposes that supernovae in compact star clusters are key sources of high-energy cosmic rays, neutrinos, and gamma rays, with simulations showing they can accelerate protons to energies above the cosmic-ray knee.

Contribution

It introduces a novel model of supernovae in compact clusters as powerful PeVatrons capable of producing very hard spectra of ultra-high-energy protons, linking them to observed gamma-ray and neutrino signals.

Findings

Simulations show acceleration of protons up to 10^{16}-10^{17} eV.

Supernova remnants in clusters may explain TeV gamma-ray and neutrino observations.

A light cosmic-ray component at 10^{17} eV could originate from these sources.

Abstract

We discuss a specific population of galactic PeVatrons which may be the main source of the galactic cosmic-ray (CR) component well above PeV energies. Supernovae in compact clusters of massive stars are proposed as powerful sources of CRs, neutrinos, and gamma-ray emission. Numerical simulations of non-linear Fermi acceleration at converging shock flows have revealed that these accelerators can provide very hard spectra of protons up to $10^{16}-10^{17}$ eV which is well above the "knee" in the all-particle CR spectrum at about $3\times10^{15}$ eV. We suggest that known supernova remnants interacting with stellar winds in the compact clusters of young massive stars Westerlund I and Cl*1806-20 can be associated with the sources of the TeV gamma-ray emission detected by H.E.S.S. and may be responsible for a fraction of the high-energy neutrinos detected with the IceCube observatory. A recent CR composition measurement with the LOFAR array has revealed a light-mass component possibly dominating the all-particle spectrum at energies around $10^{17}$ eV. Such a strong light component (mainly protons and helium) may require specific galactic CR sources such as supernovae interacting with compact clusters of massive stars in addition to isolated supernova remnants.

Equations

$$F(E)=\frac{c}{4\pi }S(E,Z)\frac{\exp \left(-\frac{d^2}{4R_d^2}\right)}{4\pi H{R}_d^2}\times$$

$$\times \sum _{n=0}^{\infty }\exp \left[-\frac{(2n-1{)}^2{\pi }^2{R}_d^2}{4H^2}\right],$$

$$\nabla (D\nabla N)+nv(p)\sigma (p)\delta (z)N(r,z,p)=-Q(p)\delta (z).$$

$$N(z,p)=\nu R{\int }_0^{\infty }dk\frac{\sinh [k(L-z)]}{\sinh (kL)}\cdot \frac{J_1(kR)}{B(p)}\cdot Q(p),$$

$$B(p)=2D(p)\coth (kL)+nv(p)\sigma (p),$$

Full text

††journal: Advances in Space Research

1 Introduction

Measurements of cosmic ray (CR) composition at energies just above the spectral knee [see e.g., Apel et al., 2013, Berezhnev et al., 2012a, Buitink et al., 2016] are of fundamental importance for deducing the origin of CRs [e.g., Hillas, 2005, Strong et al., 2007, Blasi, 2013, Amato, 2014, Blandford et al., 2014]. Recent analysis of the atmospheric depth of CR shower maxima measured with the Low Frequency Array (LOFAR) [Buitink et al., 2016] have provided evidence for a predominantly light-mass CR composition at energies around $10^{17}$ eV, i.e., well above the spectral steepening at $\sim 3\!\times\!10^{15}$ eV referred to as the “knee.”

The knee in the all-particle CR spectrum is believed to result from a rigidity (i.e., charge) dependent cut-off in shock acceleration occurring in isolated supernova remnants (SNRs). If this is the case, heavy elements obtain a higher energy with the iron cut-off reaching 26 times that of protons. The overall composition in the knee region should become heavier with energy and be dominated by heavy elements around $10^{17}$ eV. If the LOFAR observations showing a large light-mass fraction are confirmed, an additional source of light nuclei may be needed as we consider here.

To date several models of particle accelerators which are able to accelerate CRs beyond the knee have been suggested. Some of them consider CR acceleration at the galactic wind termination shock [e.g., Jokipii and Morfill, 1987], while others deal with individual isolated supernovae with varying explosion energies [e.g., Sveshnikova, 2003] and in different circumstellar environments including exploding Wolf-Rayet (WR) stars in stellar cavities [Voelk and Biermann, 1988], core-collapse supernovae (SNe) with dense circumstellar winds [Ptuskin et al., 2010], or type Ibc supernovae with trans-relativistic shocks [e.g., Ellison et al., 2013b].

Models assuming the collective effects of multiple clustered SNe and strong winds of young massive stars in galactic superbubbles produced by OB-star associations have also been considered by Bykov and Toptygin [2001], Bykov [2001], Parizot et al. [2004], Ferrand and Marcowith [2010], Bykov [2014]. These models all allow the acceleration (or re-acceleration) of CRs to beyond the knee. In any case, recent $\gamma$-ray observations of the Cygnus Superbubble [Ackermann et al., 2011] and a candidate young massive OB association/cluster G25.18+0.26 [Katsuta et al., 2017] in the Milky Way, as well as in starburst galaxies [see e.g., Bykov, 2014, Ohm, 2016], have provided confirming evidence that SNe can be efficient CR accelerators in such an environment.

It’s clear that many SNe occur in compact star clusters. Such clusters or superclusters of mass $>10^{5}\,M_{\bigodot}$ may contain hundreds or thousands of young massive stars within a parsec size core [see e.g., Clark et al., 2008, Portegies Zwart et al., 2010]. Such systems are often observed in external starburst or interacting galaxies while a few of them have been found in Milky Way-type galaxies. The most massive stars in such a cluster must begin to explode just a few million years after the cluster is born. Interacting with multiple stellar winds in a cluster environment, core-collapse SNe can accelerate particles to well above PeV energies in less than a thousand years [Bykov, 2014, Bykov et al., 2015]. A specific feature of such systems, predicted by non-linear modeling, is that the particles accelerated therein have very hard spectra at PeV energies.

2 Supernovae with colliding shock flows as efficient galactic PeVatrons

The most massive stars in compact clusters explode as core collapse supernovae after a few million years. Their shocks will travel through regions filled with fast winds from other massive young stars, or with the global cluster wind, forming colliding shock flow (CSF) systems. These CSF systems can result in an efficient conversion of the kinetic power of a supernova shell and stellar winds into the power of high-energy CRs. The first-order Fermi mechanism accounting for the non-linear backreaction of the accelerated CRs on the colliding shock structure was discussed in detail in Bykov et al. [2013, 2015]. In these parsec-sized systems, the CSF stage starts a few hundred years before the SN shock collides with the wind termination shock. At this time, the maximum energy particles that have been accelerated at the SNR shock are able to reach the fast wind termination shock and can be scattered back to the SNR by magnetic fluctuations carried by the fast stellar wind. Therefore, high-energy particles with mean free paths larger than the distance between the two shocks start to be accelerated by the converging fast flows. This is the most favorable circumstance for efficient Fermi acceleration. It results in a high pressure of accelerated particles at the highest energy end of the CR distribution, which decelerates the colliding flows.

Simulations predict non-power-law spectra with a strong upturn above several TeV becoming extremely hard at PeV energies. The upturn occurs when particles are accelerated to an energy allowing them to scatter between the SNR shock and the stellar wind shock thereby gaining energy in a fashion more efficient than in a single shock.

It was shown that during the period of the most efficient acceleration (up until several hundred years after the SN explosion) the intensity of CR production in a CSF system exceeds that of an isolated SNR in the 1014 – 1017 eV regime. Due to the fast magnetic field amplification by CR driven instabilities and efficient CR acceleration in the colliding shock flows, the spectrum of accelerated particles has a much shorter evolution time than that in an isolated SNR [Bykov et al., 2015]. In Fig. 1 both proton and electron spectra are shown at the end of the CSF acceleration phase, just before the shocks collide. The hard spectrum expected from CSFs puts most of the energy into the highest energy protons. This feature of the spectral energy distribution leads to high fluxes of $\gamma$-rays in the range between 1014–1016 eV [see red (lowest) curve in Fig. 1]. Also, CSF systems, along with the remnants of IIb and IIn [Ptuskin et al., 2010] or Ibc SNe [Ellison et al., 2013b], are unique galactic accelerators which are able to contribute importantly to the overall CR spectrum at energies above the “knee” ($3\cdot 10^{15}$ eV).

The maximal energy of a particle accelerated at a CSF system depends on the velocities of the colliding plasma flows, on details of the system geometry, and on the magnitude of the magnetic field. The proton acceleration time in a CSF system was estimated in Bykov et al. [2015] as $T_{\mathrm{acc}}\approx 2\,\cdot 10^{10}\,\epsilon_{\rm PeV}\,(\eta_{b}n)^{-0.5}\,u_{s3}^{-2}\,u_{w3}^{-1}$ s. Where $\eta_{b}<$1 is the magnetic field amplification efficiency, the CR proton energy $\epsilon_{\rm PeV}$ is in PeV, and the speeds of supernova ejecta $u_{s3}$ and cluster wind $u_{w3}$ are in units of $10^{3}$ km s*-1*. For typical parameters of a young SN in a cluster this results in a maximal energy of accelerated CRs $E_{\mathrm{max}}\sim Z\cdot 2\cdot 10^{16}\,\mathrm{eV}$, where $Z$ is the charge of the particle (see Fig. 1). The value of $E_{\mathrm{max}}$ can be much higher than the above estimate if the CSF system was initiated by a hypernova in a cluster of young massive stars. In any case, $E_{\mathrm{max}}$ exceeds the maximal CR energy for typical isolated SNRs.

3 The contribution of supernovae in stellar clusters to the CR spectrum

The unique features of CSF systems of having particularly hard spectra and short acceleration lifetimes ($\tau_{a}\leq 500$ yrs) make it possible for these sources to dominate at PeV energies without violating observational constraints at lower CR energies.

3.1 Single CR source contribution

Consider a system where a single supernova shell is colliding with the stellar wind of a nearby massive star. For a few hundred years the system will be an efficient accelerator: it will produce high-energy particles bouncing between the converging flows. Fig.1 shows trapped CR spectra and escaping fluxes for such a system where non-linear effects and magnetic field amplification have been modeled.

In the terms of the halo-diffusion model, the flux of nuclei with charge $Z$ coming from a single source located at the distance $d$ can be evaluated as follows [see e.g. Sveshnikova et al., 2013]:

[TABLE]

where $S(E,Z)$ is the flux escaping from the source (see Fig. 1), $R_{d}=\sqrt{D(E/Z)t}$ is the diffusion radius of a particle, $H=$4 kpc is the width of the Galactic halo, $t$ is the age of the source, and $D=5\cdot 10^{29}(E/\mathrm{PeV})^{0.33}$ cm2 s*-1* is the diffusion coefficient.

The total energy of a CSF source contained in CRs can be as high as 2$\cdot$1049 erg. A typical lifetime of a CSF system is $\sim 300$ years, so the $\gamma$-ray emission related to the SNR and its neighborhood is time-variable on scales of tens of years. Thus, if CSF systems actually contribute to the observed CR flux at energies beyond the knee, the $\gamma$-ray emission from these systems will be time-variable over observable epochs.

In Fig. 2 we show simulated all-particle spectra from CSF sources at 1.0, 1.2, and 1.4 kpc distances for young (3,000 yr) and relatively old (10,000 yr) systems. The modeling assumes the sources accelerate particles up to $Z\cdot$ 40 PeV. Since the acceleration mechanism depends on the CR rigidity, $R=pc/Ze$, the maximal energy of the nuclei is proportional to charge. The velocity of the SNR shock is assumed to be 6000 km s*-1*, the velocity of the stellar wind is assumed to be 4000 km s*-1*, and the magnetic field is assumed to be amplified at the shock up to 900 $\mu$G.

Fig. 2 shows that a single, relatively young CSF source can make a significant contribution to the observed CR spectrum if it is within $\sim 1.2$ kpc. A young SNR (age $\lesssim\!10$ kyr) is required both to accelerate nuclei up to the highest energies and to provide the highest possible CR fluxes. Our simulations have shown that a source within 1.4 kpc could contribute $5-10\%$ of the total CR flux between $\sim 3\!\times\!10^{16}$ and $10^{17}$ eV, while one within 1.2 kpc could contribute $\gtrsim\!30\%$ and dominate the observed bumps in this energy range.

3.2 The contribution of an ensemble of sources to the GCR spectrum

In order to estimate the total contribution of CSF systems to the overall CR spectrum we have followed a model of CR diffusion in the Galaxy [see e.g., Ginzburg and Syrovatskii, 1964, Thoudam and Hörandel, 2014]. The analytical solution of the transport equation was used to describe a steady CR distribution in the Galactic plane, which was assumed to be flat and thin [e.g., Thoudam and Hörandel, 2014]. Only primary accelerated particles in the $10^{15}-10^{17}$ eV energy range were taken into account.

A steady-state transport equation with diffusion and interaction losses for CR nuclei in cylindrical geometry can be written as

[TABLE]

Here $r$ and $z$ are the radial and vertical coordinates, with z = 0 representing the galactic plane, $p$ is the momentum of the particle, $n$ is the average number density, $v(p)$ is the velocity of a particle, $D(p)$ is the diffusion coefficient, $\sigma(p)$ is the inelastic collision cross-section, and $Q(p)$ is the source spectral energy distribution shown in Fig. 1.

The solution of (2) can be written as

[TABLE]

where

[TABLE]

$J_{1}(x)$ is the Bessel function of index 1, and $\nu$ = 2.5 sources/Myr/kpc2 is the rate of SNe in clusters with CSF systems in the Galaxy. The rate $\nu$ was assumed to be 10% of the total galactic SNe rate, which, in turn, is taken to be 3 SNe per 100 years. Thus we assume that there are $\sim 10$ CSF systems accelerating CRs in the Milky Way at any given time. These data were chosen to illustrate a possible contribution of the CSF systems to the observed CR spectrum. More detailed observational information on the galactic SNe in young massive star clusters is needed to justify or reject the model assumptions.

In Fig. 3 we compare our galactic CSF model predictions to LOFAR data as presented by Buitink et al. [2016]. We assume that the maximum energy of CR protons produced by a single source varies from $20$ to $60$ PeV. As discussed above, such energies can be reached within $300-400$ yr.

The best fit to the LOFAR data is achieved for a four-component model of protons, helium, nitrogen, and iron nuclei with a fraction of protons and helium of 0.8 [Buitink et al., 2016]. In order to illustrate the possible contribution of the light component accelerated at CSF sites to the overall CR spectrum in Fig. 3, simulated spectra for a mixture of $50\%$ proton and $50\%$ helium cumulatively produced by galactic CSFs are shown for different distributions of $E_{\mathrm{max}}$. The simulations assume the efficiency of energy conversion into CRs is $\eta=0.3$ and with the above parameters, our results show that a set of galactic CSF systems in massive stellar clusters can dominate the observed CR flux in the energy range $10^{16.5}-10^{17}$ eV. We propose this as the second galactic component which is suggested by LOFAR measurements.

4 Neutrino and $\gamma$-ray emission from supernovae in stellar clusters

Previously we have modeled CRs and neutrinos produced by colliding shocks in the compact star cluster Westerlund 1. Using parameters for the presumed supernova progenitor of the magnetar CXO J164710.2-455216, we concluded that this galactic source could explain a fraction of the recently detected IceCube neutrinos [Bykov et al., 2015].

We now consider a supernova which produced the 650 year old magnetar SGR 1806-20, which is a member of the massive stellar cluster Cl* 1806-20 [Fuchs et al., 2000, Corbel and Eikenberry, 2004]. The soft-gamma-repeater SGR 1806-20 is known for its giant flare in 2004 and is one of the most interesting and active soft gamma repeaters (SGRs). The age of the cluster Cl* 1806-20 is estimated to be about 3-4 Myr and it contains many massive stars including four Wolf-Rayet (WR) stars, five O-type stars, and a luminous blue variable candidate LBV 1806-20 [van Kerkwijk and Kulkarni, 1995]. The supernova shock that must have been launched after the birth of SGR 1806-20 has had enough time to cross the stellar cluster and may provide the conditions required for efficient particle acceleration according to the CSF model discussed above.

Given the similarity of the $\gamma$-ray source HESS J1808-204 with previously detected H.E.S.S. sources in massive star associations Cyg OB2, Westerlund 1, and Westerlund 2  [Aharonian et al., 2002, Abramowski et al., 2012], we have applied the CSF particle acceleration model to the stellar cluster Cl* 1806-20 taken to be at a distance of 8.7 kpc. In Fig. 4 we compare predicted $\gamma$-ray and neutrino spectra with H.E.S.S. and IceCube observations. The predicted spectra include emission from CRs still trapped in the source as well as those that have escaped from the accelerator and diffused into the surrounding cloudy medium over the last 650 yr. The trapped and escaped CR protons and electrons used to generate Fig. 4 are shown in Fig. 1. The upward curvature in the neutrino spectrum above $\sim 10$ TeV reflects the transition from CR acceleration at the single SNR shock to the more efficient acceleration of high-energy particles as they scatter back and forth between the SNR shock and the cluster wind.

The diffusion coefficient within $\sim$ 30 pc of the CSF accelerator was taken to be $D=5\cdot 10^{29}(E/\mathrm{PeV})^{0.33}$ cm2 s*-1* which is broadly consistent with the ISM values as discussed by Strong et al. [2007]. The number density in the cloudy medium in the dense region around the accelerator was assumed to be $\sim 90$ cm*-3*. This value seems to be justified even if the cluster Cl* 1806-20 is located away from the Galactic center region because the young cluster likely originated from a molecular cloud which should have a much longer life time and may extend for $\sim$ 20-30 pc.

A possible association of Cl* 1806-20 and SGR 1806-20 with the Galactic center region is still a matter of a dispute. To construct the spectra in Figure 4 we adopted a distance of $8.7$ kpc following the analysis by Bibby et al. [2008]. However, Svirski et al. [2011] constrained the distance to SGR 1806-20 with a lower limit of 9.4 kpc and an upper limit of 18.6 kpc (90 per cent confidence). Tendulkar et al. [2012], who measured the proper motion of SGR 1806-20, has adopted the distance to SGR 1806-20 to be 9 $\pm$ 2 kpc, which is consistent with all the measurements available.

The IceCube data points in Fig. 4 show the neutrino energy flux consistent with the position of Cl* 1806-20 at the angular resolution of IceCube [Aartsen et al., 2014]. Thus, the CSF model applied to HESS J1806-204 can explain a subset of the observed IceCube neutrinos in the inner Galaxy. As we show in Fig. 5, one of the PeV events (14, “Bert”) from recent neutrino detections is within $1\sigma$ angular distance from the position of Cl* 1806-20 as determined from $2D$ Gaussian statistics and the median angular errors and positions given by Aartsen et al. [2014].

Primary electrons accelerated in the CSF system and secondary electrons from $\pi^{0}$-decay are both contributing to the inverse Compton (IC) spectrum in Fig. 4. The ISRF field near the Galactic Center needed for the simulation of the IC spectrum was adopted from a GALPROP simulation [Porter et al., 2008]. Summing the lepton emission with the hadron emission leads to the softening of the resulting spectrum in the 100 GeV – 1 TeV energy range, fitting the H.E.S.S. data. For completeness we also show in Fig. 4 the synchrotron emission from primary and secondary electrons accelerated in the CSF system.

For illustration in Fig. 5 we show a galactic map for a subset of neutrinos detected by IceCube in the vicinity of Westerlund 1 and Cl1806-20. Two-$\sigma$ contours (applying Gaussian statistics) are shown for PeV events and one can see that two PeV events 14 (“Bert”) and 35 (“Big Bird”) are within 2$\sigma$ from both Westerlund 1 and Cl 1806-20. We note that while no time or space clumping has been confirmed so far, it has been reported that not less than $25-30\%$ of the IceCube events should originate from the Galactic plane [IceCube Collaboration et al., 2017]. In another analyzes, Denton et al. [2017] pointed out that the maximum likelihood function for the origin of neutrino events corresponds to a $6.6\%$ portion of Galactic neutrinos. This analyzes shows that among the full set of IceCube neutrinos, events 2, 14 and 52 have the highest probabilities for being galactic with 0.42, 0,75 and 0.19 probabilities respectively.

5 Conclusions

The expanding young supernova shells in compact stellar clusters can be very efficient CR accelerators to energies well above the standard limits of diffusive shock acceleration in an isolated SNR. During a relatively short time of a few hundred years after the supernova explodes, the SNR shock can collide with the fast wind from a massive star, or with the cluster wind, and produce Fermi accelerated CRs to energies above $10^{16}$ eV. The energy spectrum of protons in such a colliding shock flow accelerator is hard, approaching the spectral shape $dN/dp\propto p^{-1}$ before a break, providing a large energy flux in the high-energy end of the spectrum. While the strongest acceleration stage lasts for only a few hundred years, high-energy CRs escape the CSF system and diffuse through the ambient matter producing $\gamma$-rays and neutrinos in inelastic nuclear collisions over a much longer period. We have modeled a galactic population of CSF sources and conclude CSF systems would contribute a significant light-mass component to galactic CRs above $\sim 10^{16}$ eV. We compare our predictions to recent LOFAR observations [Buitink et al., 2016] which suggest a need for such a component. According to the analyses of Buitink et al. [2016], the total fraction of light elements (p and He) in LOFAR measurements lies between 0.38 and 0.98 for a significance level $p>0.01$. The light component mix is likely to have both protons and He with an excess of He nuclei. We have simulated a 50/50 mix of protons and He, but a slight change in this ratio will not strongly affect our result.

Considering neutrinos, the recent data release of the IceCube Collaboration has confirmed that no individual steady-state or transient source of very high-energy neutrinos has been identified so far [IceCube Collaboration et al., 2017]. An analysis of the data obtained during 6 years suggests that AGN jets or GRBs are likely to be responsible for only a small fraction of the observed neutrinos [IceCube Collaboration et al., 2017]. Thus, there is a substantial need for alternative models of efficient sources of very high-energy cosmic neutrinos. We have proposed that CSF systems in galactic compact clusters of massive stars may produce CRs and neutrinos efficiently and that two such clusters, Westerlund 1 and Cl1806-20, are spatially consistent with a fraction of the IceCube neutrinos. Events 2, 14, and 52 being within a 2-$\sigma$ distance from both Westerlund 1 and Cl1806-20 have the highest probabilities for being galactic with 0.42, 0.75, and 0.19 probabilities respectively. Consistently with the neutrino production, we have shown that the CRs in these two clusters can account for the TeV $\gamma$-ray emission detected by H.E.S.S.

Recently, Neronov and Semikoz [2016] showed that the galactic latitude distribution of the astrophysical neutrino events with energies above 100 TeV detected with IceCube is apparently inconsistent with an isotropic distribution. This supports our conjecture that SNe in galactic compact clusters, namely Westerlund 1 and Cl*1806-20, may explain $\gtrsim\!10\%$ of the observed neutrino flux. Moreover, compact clusters are numerous in starburst galaxies. The starburst galaxies have been proposed by Loeb and Waxman [2006] as important sources of extra-galactic neutrinos and these are widely discussed in the context of the IceCube data [see also Tamborra et al., 2014, Emig et al., 2015, Bykov, 2014, Bechtol et al., 2017]. We believe it is possible that the colliding shock flow model we propose is the main PeV CR generator in starburst galaxies and therefore responsible for a substantial fraction of extra-galactic IceCube neutrinos as well. The development of the next generation instrument, IceCube-Gen2, due to its five times better sensitivity for individual sources and not less than ten times larger surface area, will allow detailed studies of high-energy cosmic neutrinos, identify their powerful sources, and investigate the energy conversion processes therein [IceCube-Gen2 Collaboration et al., 2014].

6 Acknowledgements

We thank the anonymous referee for constructive comments. The results of the work were obtained using computational resources of Peter the Great Saint-Petersburg Polytechnic University Supercomputing Center.