High-energy cosmic rays from compact galactic star clusters: particle fluxes and anisotropy

TL;DR

This paper explores how compact galactic star clusters could produce high-energy cosmic rays and influence their observed fluxes and anisotropy, suggesting they may significantly contribute to cosmic rays around 100 PeV without conflicting with anisotropy observations.

Contribution

It introduces a model of cosmic-ray production from star clusters and examines their impact on flux and anisotropy, extending understanding of cosmic-ray origins at ultra-high energies.

Findings

Star clusters can produce cosmic rays above 10^{17} eV.

Galactic star clusters may account for a substantial fraction of 100 PeV cosmic rays.

The model aligns with observed cosmic-ray anisotropy levels.

Abstract

It has been shown that supernova blast waves interacting with winds from massive stars in compact star clusters may be capable of producing cosmic-ray (CR) protons to above $10^{17}$ eV. We give a brief description of the colliding-shock-flows mechanism and look at generalizations of the diffusion of ~ 100 PeV CRs in the turbulent galactic magnetic field present in the galactic disk. We calculate the temporal evolution of the CR anisotropy from a possible distribution of young compact massive star clusters assuming the sources are intermittent on time scales of a few million years, i.e., comparable to their residence time in the Milky Way. Within the confines of our model, we determine the galactic/extra-galactic fraction of high-energy CRs resulting in anisotropies consistent with observed values. We find that galactic star clusters may contribute a substantial fraction of ~ 100 PeV CRs without producing anisotropies above observed limits.

Equations

$${\mathbf{B}}_{turb}(\mathbf{r})=\sum _{n=1}^N{A}_n{e}^{i({\mathbf{k}}_{\mathbf{n}}\cdot \mathbf{r}+{\psi }_n)}{\hat{\xi }}_n,$$

$${\hat{\xi }}_n=\cos ({\beta }_n){\hat{\mathbf{e}}}_n^2+i\sin ({\beta }_n){\hat{\mathbf{e}}}_n^3,$$

$$J(\mu )=J_0+J_1\mu$$

$$|A|=\frac{J_{\max }-J_{\min }}{J_{\max }+J_{\min }}=\left|\frac{J_1}{J_0}\right|.$$

$$\overline{u_{z^{\prime }}}=\frac{\int u_{z^{\prime }}J(\mu )d^{3}p}{\int J(\mu )d^{3}p}=\frac{{\int }_{-1}^{1}c{J}_1{\mu }^{2}d(\mu )}{{\int }_{-1}^1{J}_{0}d(\mu )}=\frac{c{J}_1}{3J_0},$$

$$\overline{u_{x^{\prime }}}=\overline{u_{y^{\prime }}}=0.$$

$$D=\frac{{⟨\Delta r⟩}^2}{6\Delta t},$$

$$\lambda =3D/c,$$

$$\Delta {\theta }_{\max }=\pi ,$$

$$\Delta t=2\lambda c.$$

$$\begin{gathered} A_x=3\frac{\overline{v_x}}{c}= \\ 3\frac{{\sum }_{i=1}^N{v}_x^i/c}{N}, \end{gathered}$$

$$A_y=3\frac{\overline{v_y}}{c}=3\frac{{\sum }_{i=1}^N{v}_y^i/c}{N},$$

$$A_z=3\frac{\overline{v_z}}{c}=3\frac{{\sum }_{i=1}^N{v}_z^i/c}{N},$$

$$A_{\mathrm{tot}}=\frac{A_{\mathrm{CSF}}}{1+f_{\mathrm{ex}}},$$

$$\frac{F_{\mathrm{gal}}}{F_{\mathrm{ex}}}=\frac{1}{f_{\mathrm{ex}}}=\frac{A_{\mathrm{tot}}}{A_{\mathrm{CSF}}-A_{\mathrm{tot}}}.$$

Full text

††journal: Advances in Space Research

1 Introduction

The main source of galactic cosmic rays (CRs) with energies in the TeV energy regime is almost certainly supernova remnants (SNRs) which occur on time scales which are short compared to the CR residence time in the disk [see e.g. Berezinskii et al., 1990, Baring et al., 1999, Aharonian et al., 2012, Strong et al., 2007, Drury, 2012, Helder et al., 2012, Amato, 2014, Blandford et al., 2014, Bell, 2015, Grenier et al., 2015]. The source of CRs with energies in the PeV-EeV range is far less certain. Isolated SNRs with typically observed characteristics are unlikely to produce PeV-EeV CRs although relativistic supernovae (SNe) may be able to produce CRs well above the “knee” region at $10^{15-16}$ eV [e.g., Wang et al., 2007, Chakraborti et al., 2011, Bykov et al., 2018a]. The energy budget requirement for the population of CRs above the knee is less stringent than for lower energy CRs and the sources of PeV-EeV CRs may be fundamentally different from isolated SNRs. Here we consider CRs produced in colliding shock flows (CSFs) in galactic compact star clusters as proposed by Bykov [2014], Bykov et al. [2018b].

While many SNe may occur in isolation, it’s clear that SNe also occur in compact star clusters including clusters with masses $>10^{5}M_{\odot}$. These massive clusters may contain thousands of young massive stars within the few parsec scale size of the cluster. The most massive stars will explode within a few million years from the birth of the cluster and it has been shown [e.g. Bykov, 2014, Bykov et al., 2015, Lingenfelter, 2018, Grimaldo et al., 2019] that, if the blast waves from these core-collapse SNe interact with the strong stellar winds from nearby massive stars, efficient Fermi acceleration can produce protons to hundreds of PeV with an extremely hard spectrum near the cutoff at the maximum obtained energy [see Fig. 2 from Bykov et al., 2018b].

Furthermore, if the first-order Fermi shock acceleration mechanism accelerates heavy ions in CSF systems in the same way as believed to be the case in other shocks [e.g., Ellison et al., 1997], the high-energy peak shown in Fig. 2 will be shifted upward by a factor equal to the CR charge. Iron nuclei might be accelerated to EeV energies.

Aharonian et al. [2019] derived an approximate $1/r$ decrement of the CR density with distance from a star cluster from the analysis of very high energy gamma-ray observations of massive star clusters like Westerlund 1 and 2, the Cygnus cocoon, and others. This may indicate a continuous high energy CR injection into the interstellar medium over a few million years. They pointed out that such sources may operate up to the CR knee around 1 PeV. This is consistent with the modeling of high energy CR spectra in compact clusters of young massive stars by Bykov [2014], Bykov et al. [2015, 2018b] discussed above. A critical constraint on the possible contribution of galactic sources to PeV-EeV CRs is the low observed anisotropy. Here we calculate CR diffusion in the Milky Way with a simulated turbulent magnetic field, as described in § 2.1. The CRs are injected from sources mimicking the known massive star clusters and are propagated using a Monte Carlo method described in § 2.3.

2 Model

2.1 Magnetic field and diffusion coefficients

To calculate the galactic magnetic field, we use the most recent model of the turbulent magnetic field suggested by Han [2017]. In this model, it is considered that on small spatial scales both density and magnetic field fluctuations follow a Kolmogorov spectrum, as shown with the dashed line in Fig. 3. By using rotation measures (RMs) and dispersion measures (DMs) of pulsars, Han derived the spatial energy spectrum of the Galactic interstellar magnetic field at scales from $1/k=0.5$ to $15$ kpc, where $k$ is the wavenumber. It was found that the spectrum in this scale range is much flatter than the Kolmogorov spectrum extended to small scales and the turnover occurs between $0.5$ kpc and $80$ pc. Taking this into account, we use the spectrum shown in Fig. 3 for our modeling.

We model the turbulent magnetic field as a sum of plain waves with random directions, phases, and polarizations [e.g., Giacalone and Jokipii, 1994, Casse et al., 2002], i.e.,

[TABLE]

where $\mathbf{r}$ is the position vector, $\mathbf{k}_{n}\equiv k_{n}\hat{\mathbf{e}}^{1}_{n}$, $A_{n}$, $\psi_{n}$ and $\hat{\xi}_{n}$ are the wavevector, amplitude, phase and polarization vector of each mode, respectively. The polarization vector is given by:

[TABLE]

where $\beta_{n}$ is the polarization angle. $(\hat{\mathbf{e}}^{1}_{n},\hat{\mathbf{e}}^{2}_{n},\hat{\mathbf{e}}^{3}_{n})$ are orthogonal, so $\mathbf{k_{n}}\cdot\hat{\xi}_{n}=0$ and $\nabla\cdot\mathbf{B}_{\mathrm{turb}}=0$. The amplitude of each mode is determined with the energy spectrum shown in Fig. 3 and the root-mean-square field is $6$ $\mu$G. We also discuss below a possible effect of the structure of the regular large-scale galactic magnetic field on the anisotropy of $\sim 10^{17}$ eV CRs accelerated by galactic supernovae in the compact star clusters.

2.2 Cosmic-ray anisotropy

To determine the dipole anisotropy at a given angle $\Theta$, we assume the CR intensity $J(\mu)$ is

[TABLE]

where $\mu=\cos\Theta$. The maximum intensity occurs when $\mu=1$ or $\Theta=0$. The dipole anisotropy, $|A|$, is defined as the maximum and minimum intensity contrast, i.e.,

[TABLE]

We note that eqn. (4) only applies when $J_{1}\ll J_{0}$, i.e., when $|A|\ll 1$. The Auger CR observatory, and other experiments [e.g., Góra, 2018, Apel et al., 2019], show that the CR anisotropy at energies of $10^{17}$ eV and above does not exceed $3-5$% so eqn. (4) applies.

In a coordinate system where $z^{\prime}$ is the direction of the maximum of intensity $I(\mu)$, and therefore the direction of the mean velocity, the mean velocity $\overline{u_{z^{\prime}}}$ is

[TABLE]

and

[TABLE]

In the above, we assume the CR velocity $u\simeq c$ so $|A|\simeq|3\overline{u_{z^{\prime}}}/c|$. In the coordinates of our galaxy model, the $x$-$y$ plane is the galactic plane, the $z$-axis is directed upward, and $|\overline{u_{z^{\prime}}}|=\sqrt{\overline{v_{x}}^{2}+\overline{v_{y}}^{2}+\overline{v_{z}}^{2}}$.

2.3 Monte Carlo cosmic-ray propagation

Given the simulated turbulence discussed in § 2.1, we model the particle propagation by solving the equation of motion for a set of particles as they move through the turbulent magnetic field. From the motion, we find the diffusion coefficient, $D$, using

[TABLE]

where $\left<\Delta r\right>$ is the mean displacement of particles during the time $\Delta t$. Our results for CR energies $3\!\times\!10^{15}-3\!\times\!10^{17}$ eV are shown in Fig. 4.

Once the diffusion coefficient is determined, the scattering mean-free-path is given by

[TABLE]

and we run the Monte Carlo simulation with the following parameters. We follow two CR energies, $10^{17}$ and $3\!\times\!10^{17}$ eV, in a galaxy of 30 kpc diameter and 10 kpc thickness. We assume supernova events occur every 2000 yr and each source ejects $10^{7}$ Monte Carlo particles in random directions. The particles undergo large-angle scattering with the mean-free-path determined by eqn. (8). The maximum angle of scattering is set as:

[TABLE]

and the time step is

[TABLE]

At the position of the Earth, we find the anisotropy and particle flux in a small sphere of radius 0.5 kpc centered at the Solar System. We follow each particle for a simulation time of $4\!\times\!10^{6}$ yr or until it leaves the Galaxy. At every point in time, each particle will have coordinates $x$, $y$, and $z$ and velocity components $v_{x}$, $v_{y}$, and $v_{z}$. With this information, we determine the CR concentration and anisotropy near Earth at each time step. The components of anisotropy are given by:

[TABLE]

[TABLE]

[TABLE]

where $N$ is the number of CRs in our region. The total anisotropy is then $|A|=\sqrt{A_{x}^{2}+A_{y}^{2}+A_{z}^{2}}$.

2.4 CR anisotropy from CSF sources

In Fig. 1 we show the distribution of known star clusters from Portegies Zwart et al. [2010]. We propagate CRs from these clusters to the Earth taking into account the large-scale structure of the Galactic magnetic field. There are several models of the large-scale regular field [see e.g. Pshirkov et al., 2011, Jansson and Farrar, 2012, Han, 2017]. Here we use the 21-parameter model by Jansson and Farrar. In our Monte Carlo simulations, CRs spiral in the large-scale regular field while scattering in the small-scale turbulence.

In Fig. 5 we show the simulated anisotropy for $10^{17}$ eV CRs. Our simulations with the parameterized regular magnetic field are very computationally demanding, therefore, these simulations had only 105 particles per SNe compared to $10^{7}$ particles for the model without the regular magnetic field shown in the left panel in Fig. 5. This results in the higher noise level of the anisotropy (blue line in the right panel in Fig. 5). The time history shows that $|A|$ reaches nearly constant minimum values after $\sim 10^{6}$ yr. As seen in the right panel, the mean anisotropy is similar with and without propagation in the regular field of Jansson and Farrar [2012].

Fig. 6 shows the simulated flux at Earth as a function of run time, assuming that the efficiency of converting the SN explosion energy ($10^{51}$ erg) to the energy of 100 PeV protons is about $0.015\%$. The presence of the regular magnetic field doesn’t make any significant change in the flux which becomes nearly constant after $\sim 1.5\!\times\!10^{6}$ yr.

The observed equatorial dipole anisotropy at a few times $10^{17}$ eV is $\lesssim\!0.01$ [e.g., Mollerach and Roulet, 2018]. For the parameters used here, we obtain $|A|\simeq 0.03$ for $E_{\mathrm{CR}}=10^{17}$ eV. An isotropic extra-galactic component of CRs will reduce the anisotropy and we can estimate the ratio of the galactic to extra-galactic flux, $F_{\mathrm{gal}}/F_{\mathrm{ex}}$, in the following way. If the isotropic extra-galactic flux is $f_{\mathrm{ex}}$ times as large as the average flux from the CSF sources, the total anisotropy is

[TABLE]

and

[TABLE]

Taking $A_{\mathrm{tot}}=0.01$, our results suggest that $F_{\mathrm{gal}}/F_{\mathrm{ex}}\sim 1/2$ for $E_{\mathrm{CR}}=10^{17}$ eV.

3 Conclusions

As shown in Bykov et al. [2018b], and references therein, the colliding shock flows (CSFs) model produces an exceptionally hard CR proton spectrum peaking at or above the CR knee (i.e., Fig. 2). A key constraint on the fraction of ultra-high-energy CRs produced by CSFs in galactic star clusters is the anisotropy. To calculate the anisotropy, we have used a simple model for the galactic geometry, interstellar magnetic turbulence, and distribution of CSF sources. The anisotropy depends strongly on the diffusion coefficient, $D(E)$, and we have determined $D(E)$ with a direct calculation of CR trajectories in simulated turbulence (Fig. 4). The model described above considered high energy CR production by a subset of galactic SNe exploding in compact clusters of young massive stars, i.e., CSF systems. The mean time between SNe in these CSF systems is taken to be 2000 yr. The clusters are distributed within the galactic disk as illustrated in Fig. 1. Cosmic rays are propagated until they leave the thick galactic disk.

We estimated the ratio of the high energy CR flux produced by young massive clusters to the total flux needed to provide the observed anisotropy and found that $30-50\%$ of the total at $\sim 100$ PeV may come from clusters without violating observed anisotropy levels.

4 Acknowledgements

A.B. and M.K. were supported in part by the RSF grant 16-12-10225. Some of the modeling was performed at the “Tornado” subsystem of the St. Petersburg Polytechnic University supercomputing center and at the JSCC RAS.