Chemical Enhancements in Shock-accelerated Particles: Ab-initio Simulations

TL;DR

This study uses hybrid simulations to explore how ions with different mass and charge ratios are thermalized and accelerated in collisionless shocks, revealing mechanisms behind heavy-ion enhancements in cosmic rays.

Contribution

It provides new insights into ion thermalization and acceleration processes, highlighting the role of ion charge and mass ratios in non-relativistic shocks.

Findings

Ions thermalize to a temperature proportional to their mass number A.

Non-thermal tails develop with extent scaling with Z and normalization with (A/Z)^2.

Heavy ions are preferentially accelerated, explaining cosmic ray observations.

Abstract

We study the thermalization, injection, and acceleration of ions with different mass/charge ratios, $A/Z$, in non-relativistic collisionless shocks via hybrid (kinetic ions-fluid electrons) simulations. In general, ions thermalize to a post-shock temperature proportional to $A$. When diffusive shock acceleration is efficient, ions develop a non-thermal tail whose extent scales with $Z$ and whose normalization is enhanced as $(A/Z)^2$, so that incompletely-ionized heavy ions are preferentially accelerated. We discuss how these findings can explain observed heavy-ion enhancements in Galactic cosmic rays.

Equations

$$f_i(E)=\frac{(\gamma -1)n{\chi }_i{\eta }_i}{E_{inj,i}}{\left(\frac{E}{E_{inj,i}}\right)}^{-\gamma },$$

$$K_{ip}\equiv \frac{f_i(E/Z_i)}{{\chi }_i{f}_p(E)}=\frac{{\eta }_i}{{\chi }_i{\eta }_p}{\left(\frac{E_{inj,i}}{E_{inj,p}}\right)}^{\gamma -1}$$

$$K_{ip}^{\mathrm{GCRs}}={\frac{{\varphi }_i}{{\chi }_i{\varphi }_p}|}_{\mathrm{TeV}}\frac{Z_i^{\gamma -1-\delta }}{A_i^{1/2}}\simeq {\frac{{\varphi }_i}{{\chi }_i{\varphi }_p}|}_{\mathrm{TeV}}\frac{Z_i^{1/6}}{A_i^{1/2}}$$

Full text

Introduction.— Non-relativistic shocks are well-known as sources of energetic particles. Prominent examples of such shocks are the blast waves of supernova remnants (SNRs), which are thought to be the sources of Galactic cosmic rays (GCRs) (e.g., Caprioli et al., 2010; Morlino and Caprioli, 2012), and heliospheric shocks, where solar energetic particles (SEPs) are measured in situ (e.g., Reames, 2015; Desai et al., 2016). Chemical abundances in GCRs and SEPs provide crucial information about their sources and the processes responsible for their acceleration.

At trans-relativistic energies, the chemical composition of GCRs roughly resembles the composition of the solar system (Simpson, 1983), the most evident deviation being the enhancement of secondaries produced by spallation of primary GCRs during their propagation in the Milky Way. A more careful analysis, however, reveals that the GCR composition is controlled by volatility and mass/charge ratios: refractory elements show larger enhancements than volatile ones, and heavier volatile elements are more abundant than lighter ones Meyer et al. (1997); Ellison et al. (1997). Moreover, elemetns with low first ionization potential tend to be overrepresented in GCRs (e.g., Simpson, 1983). At TeV energies, where spallation is negligible, the fluxes of H, He, C-N-O, and Fe do not differ by more than one order of magnitude (e.g., Caprioli et al., 2011, and references therein). Since their typical solar number abundances relative to H are $\chi_{He}=0.0963$, $\chi_{CNO}=9.54\times 10^{-4}$, $\chi_{Fe}=8.31\times 10^{-5}$ (Lodders, 2003), the abundances observed in GCRs suggest that heavy ions must be preferentially injected and accelerated compared to protons.

Diffusive shock acceleration (DSA) (e.g., Bell, 1978; Blandford and Ostriker, 1978) at SNR shocks is likely the mechanism responsible for ion acceleration up to $\sim 10^{17}$ eV (Caprioli, 2015). DSA produces universal power-law momentum spectra $f(p)\propto p^{-3r/(r-1)}$, where $r$ is the shock compression ratio; for strong shocks $r\to 4$ and $f(p)\propto p^{-4}$. For relativistic particles the energy spectrum is then $f(E)=4\pi p^{2}f(p)dp/dE\propto E^{-2}$, while at non-relativistic energies one gets $f(E)\propto E^{-3/2}$ (Caprioli and Spitkovsky, 2014a).

Hybrid simulations.— In order to study ab initio how DSA of ions with different mass/charge ratio works, we performed 2D kinetic simulations with dHybrid, a massively parallel hybrid code, in which ions are treated kinetically and electrons as a neutralizing fluid (Gargaté et al., 2007). Hybrid simulations of non-relativistic shocks have been extensively used for assessing the efficiency of proton DSA (Caprioli and Spitkovsky, 2014a), the generation of magnetic turbulence due to plasma instabilities driven by accelerated particles (Caprioli and Spitkovsky, 2014b), the diffusion of energetic particles in such self-generated magnetic fields (Caprioli and Spitkovsky, 2014c), and the injection of protons into the DSA process Caprioli et al. (2015). In the literature there are few examples of kinetic simulations with heavy ions, e.g., the pioneering 1D hybrid simulations of weak shocks including $\alpha-$particles (Burgess, 1989; Trattner and Scholer, 1991) and the recent hybrid study of the thermalization of weakly-charged ions at shocks Kropotina et al. (2016). However, a self-consistent kinetic characterization of ion enhancement in DSA has never been performed before (see, e.g, Ellison et al., 1990; Caprioli et al., 2011, for Monte Carlo and semi-analytical approaches).

In the hybrid simulations presented in this Letter we include additional ion species characterized by number abundances $\chi_{i}$, atomic mass $A_{i}$, and charge $Z_{i}$ (in proton units), initially in thermal equilibrium with protons and electrons. We fix $\chi_{i\neq H}=10^{-5}$ to effectively make ions other than protons dynamically unimportant. Lengths are measured in units of $c/\omega_{p}$, where $c$ is the speed of light and $\omega_{p}\equiv\sqrt{4\pi ne^{2}/m}$, with $m,e$ and $n$ the proton mass, charge and number density; time is measured in units of $\omega_{c}^{-1}\equiv mc/eB_{0}$, $B_{0}$ being the strength of the initial magnetic field; velocities are normalized to the Alfvén speed $v_{A}\equiv B/\sqrt{4\pi mn}$, and energies to $E_{sh}\equiv mv_{sh}^{2}/2$, with $v_{sh}$ the velocity of the upstream fluid in the downstream frame. We account for the three spatial components of the particle momentum and of the electric and magnetic fields. Shocks are produced by sending a supersonic flow against a reflecting wall and are characterized by their sonic and Alfvénic Mach numbers $M_{s}\equiv v_{sh}/c_{s}$, $M_{A}\equiv v_{sh}/v_{A}$, with $c_{s}$ the sound speed; in this work we consider $M_{s}\simeq M_{A}\equiv M$. The shock inclination is defined by the angle $\vartheta$ between the direction of ${\bf B}_{0}$ and the shock normal, such that $\vartheta\lesssim 45^{\circ}$ corresponds to quasi-parallel shocks.

The time-step is chosen as $\Delta t=0.01/M\omega_{c}^{-1}$ and the computational box measures $2.5\times 10^{4}c/\omega_{p}$ by $2Mc/\omega_{p}$, with two cells per ion skin depth. In order to suppress the numerical heating that can arise in long-term simulations with species of disparate densities, we use 100 protons per cell, and 4 particles per cell for all the other species. We have checked the convergence of our results against 3D simulations, time and space resolution, number of particles per cell, and transverse size of the simulation box (see also Caprioli and Spitkovsky, 2014a). The electron pressure is a polytrope with an effective adiabatic index chosen to satisfy the shock jump conditions with thermal equilibration between downstream protons and electrons (Caprioli et al., 2015).

Our benchmark case comprises ion species with $A=\{1,2,4,8\}$ and $Z=\{1,2\}$ and a quasi-parallel ($\vartheta=20^{\circ}$) shock with $M=10$, which exhibits efficient proton DSA and magnetic field amplification (Caprioli and Spitkovsky, 2014a, b). In our case we find that $\sim 10\%$ of the shock kinetic energy is converted into accelerated protons, and the field is amplified by a factor of $\gtrsim 2$ in the upstream. The downstream spectra of different ion species are shown in Fig. 1, as a function of $E/Z$ and normalized to their abundances $\chi_{i}$. The color code gathers species with the same $A/Z$, while solid and dashed lines correspond to $Z=1$ and 2, respectively. Each of the species shows a thermal peak plus a power-law tail with the universal DSA slope $\gamma\simeq 3/2$; non-thermal spectra roll over at a maximum energy $E_{max,i}$, which increases linearly with time (Caprioli and Spitkovsky, 2014c). For strong shocks, Rankine–Hugoniot conditions return a downstream thermal energy $\mathcal{E}\simeq 0.6E_{sh}$ (Caprioli and Spitkovsky, 2014a). Since half of the post-shock proton energy goes into electron heating by construction, we expect $\mathcal{E}_{H}\simeq\mathcal{E}/2$. Then, since heavier ions have more kinetic energy to convert into thermal energy, their temperature is expected to scale with their masses, i.e., $\mathcal{E}_{i\neq H}=A_{i}\mathcal{E}$. Dotted lines in Fig. 1 correspond to Maxwellian distributions with such expected temperatures: they provide a good fit for the positions of thermal peaks, but only a rough one for the shape of the thermal distributions of heavy ions, whose relaxation is still ongoing 111We checked that the scaling $\mathcal{E}_{i}\propto A_{i}$ is recovered also for protons if the electron pressure is set to zero..

When comparing different ion curves in Fig. 1, we notice three important scalings:

At fixed $Z$, the thermal peaks are shifted to the right linearly in $A$, i.e, each species thermalizes at a temperature proportional to its mass (see also Kropotina et al., 2016); 2. 2.

All the ion spectra rollover at the same $E_{max}/Z$, consistent with the fact that DSA is a rigidity-dependent process 222Rigidity is defined as $p/Z$, not $E/Z$: the two definitions are equivalent only for relativistic particles. However, we showed in (Caprioli and Spitkovsky, 2014c) that the self-generated diffusion coefficient and, in turn, the acceleration time do scale as $E/Z$ in the non-relativistic case (see also Ellison and Ramaty, 1985).; 3. 3.

The normalization of the non-thermal spectra at given $E/Z$ is an increasing function of the mass/charge ratio, which implies that the efficiency of injection into DSA depends on $A/Z$.

The first two results validate the theoretical expectations, while the last one represents the first self-consistent characterization of the parameter that regulates the injection of ions into the DSA process.

Injection enhancement in DSA.— In this section we discuss how the observed boost in ion injection depends on $A/Z$. The ion non-thermal spectra, neglecting the cutoffs, are power laws that can be written as

[TABLE]

where $\eta_{i}$ is the fraction of ions that enter DSA above the injection energy $E_{inj,i}$. We then introduce the ratio

[TABLE]

as a measure of the enhancement in energetic ions with respect to protons at fixed $E/Z$. $K_{ip}$ is promptly read from Fig. 1 by taking the ratio of the power-law spectra at any $E/Z$ between 10 and 100$E_{sh}$. Note that the enhancement has two contributions: one straightforward, $\eta_{i}/\eta_{p}$, which depends on the fraction of particles that enter DSA for each species, and one more subtle that depends on $E_{inj,i}$, which cannot be predicted analytically.

Fig. 2 shows the enhancements obtained for shocks with $\vartheta=20$ and $M=\{5,10,20,40\}$; injection fractions and enhancements are calculated at time $t=10^{3}\omega_{c}^{-1}$, when DSA spectra have been established, by considering the post-shock spectra of species with $A/Z$ up to 8, integrated over $10^{3}c/\omega_{p}$.

For shocks with $M\gtrsim 10$, where accelerated protons generate non-linear upstream magnetic turbulence with $\delta B/B_{0}\gtrsim 1$, the fraction of injected particles is $\eta_{p}\approx 1\%$ for protons and increases linearly with $A/Z$ (top panel); at the same time, $K_{ip}\propto(A/Z)^{2}$, attesting to a very effective enhancement of particles with large charge/mass (bottom panel). The scaling with $A/Z$ is weaker for the lowest-$M$ shock, for which $\delta B/B_{0}\approx 0.2$: $\eta_{i}$ is roughly constant at the percent level and $K_{ip}\propto A/Z$.

Chemical enhancements.— The high-$M$ case is relevant, e.g., for SNR shocks propagating into the warm interstellar medium (ISM), where atoms are typically singly ionized. Ions that are injected into DSA will then be stripped of their electrons while being accelerated up to $\sim$PV rigidities (Morlino, 2009; Caprioli et al., 2010). In the popular scenario in which GCRs are produced at SNR shocks via DSA (Caprioli, 2015), we can compare our findings with the chemical enrichment measured in GCRs (Meyer et al., 1997; Caprioli et al., 2011). In order to compare observations at Earth and shock injection simulations, we take the observed GCR flux ratios at 1 TeV, $\phi_{i}(E)$ (e.g., table 1 in ref. Caprioli et al., 2011), weigh them with the fiducial solar abundances, $\chi_{i}$ (Lodders, 2003), and write the enhancement at a given $E$ as $K_{ip}Z_{i}^{1-\gamma}$ (see Eq. 2). We also account for the rigidity-dependent residence time in the Galaxy $\propto(E/Z)^{-\delta}$, with $\delta\simeq 1/3$ above a few GV (Blasi and Amato, 2012), and extrapolate the enhancements down to the non-relativistic injection energies. Such an extrapolation introduces an additional factor $A_{i}^{-1/2}$, because DSA spectra are power laws in momentum and hence energy spectra flatten by $E^{1/2}$ at $\sim A_{i}$GeV. Finally, we obtain that ion injection into DSA must be enhanced at SNR shocks according to

[TABLE]

in order to explain the abundances observed in GCRs.

We consider a strong quasi-parallel shock with $M=20$ and singly-ionized He, CNO, and Fe atoms with effective $A/Z=\{4,14,56\}$ and calculate $K_{ip}$ in the upstream, since at $t=10^{3}\omega_{c}^{-1}$ ions $A/Z\gtrsim 14$ have already been over-injected but have not yet developed the universal downstream DSA spectrum. The enhancements found in simulations and those in GCR data (Eq. 3) are compared in Fig. 3: the scaling $K_{ip}\simeq(A/Z)^{2}$ found for strong shocks provides a very good fit, with singly-ionized He, CNO, and Fe particles enhanced by a factor of about ten, hundred, and a few thousand, respectively. It is remarkable that such a Fe enhancement requires a very large fraction $\eta_{Fe}\lesssim 50\%$ of the pre-shock Fe ions to enter DSA; this may have implications for the overall ISM chemical composition, since regions processed by shocks may become depleted in heavy elements. Nevertheless, in the ISM many elements are typically trapped in molecules (C,O) and dust grains (Fe, refractory elements), so that fragmentation and sputtering may represent crucial steps in the injection of heavy elements (Meyer et al., 1997). Our results suggest that dust grains with very large $A/Z\gg 1$ should also have no problem of being efficiently energized via DSA-like processes, thereby sputtering pre-energized ions that can be easily injected (Ellison et al., 1997; Slavin et al., 2015).

In the low-$M$ regime relevant to heliospheric shocks, our simulations show that DSA can account for enhancements by factors of a few to ten, which are commonly observed in SEP events (e.g., Mason et al., 1999; Desai et al., 2016, and references therein). However, chemical enhancements in SEPs may be time dependent (e.g., Reames, 2015) and greatly vary from event to event; in addition to shock strength and inclination, they seem to depend on the presence of pre-existing magnetic turbulence and energetic seed particles (produced, e.g., in solar flares) (Tylka et al., 2005), which makes it nontrivial to compare individual SEP events with our simulations where ion injection only occurs from the thermal pool.

Dependence on shock inclination.— Oblique shocks with $\vartheta\gtrsim 50^{\circ}$ cannot inject thermal protons and drive self-generated magnetic turbulence Caprioli and Spitkovsky (2014a); Caprioli et al. (2015). We find that in such shocks ions with large $A/Z$ do thermalize at a temperature $\propto A$, but progressively further in the downstream with respect to the quasi-parallel case. Neither protons nor heavier ions are injected into DSA or develop a non-thermal tail, which confirms that having a large gyroradius ($\propto A/Z$) is not a sufficient condition for being injected into DSA.

The injection mechanism.— Ion injection occurs in a qualitatively different way than proton injection, which is due to specular reflection off the time-dependent potential barrier at the shock and energization via shock-drift acceleration Caprioli et al. (2015). Unlike protons, heavy ions are not halted by the shock barrier and always penetrate in the downstream for at least one gyroradius ($\sim MA/Zc/\omega_{p}$); here, their distribution tends to become more isotropic due to the presence of rapidly-varying fields, an analog of the violent relaxation in stellar dynamics. If isotropization is rapid enough with respect to advection, there arises a population of backstreaming ions that can overrun the shock barrier, which is “tuned” for preventing downstream thermal protons from returning upstream. The fraction of injected heavy ions is thus controlled by how rapid isotropization is, which depends on $A/Z$ and on the strength of the magnetic turbulence in the shock layer.

Fig. 4 shows the $x-p_{x}$ phase spaces for our benchmark run; we consider ions with $Z=1$ and $A=\{1,2,4,8\}$ at times $t=\{200,700\}\omega_{c}^{-1}$. We see that, while protons are promptly isotropized behind the shock (dashed vertical lines), ions with larger $A$ tend to retain their anisotropy further in the downstream. At early times (left column), protons show the characteristic non-thermal, isotropic population of particles diffusing around the shock (Caprioli and Spitkovsky, 2014a); ions with $A/Z=2$ have also started being injected and accelerated. Injection may be quite “bursty” when patches of quasi-parallel magnetic field are advected through the shock (Caprioli et al., 2015), resulting in coherent batches of particles protruding back into the upstream (as for $A/Z=2$ at $t=200\omega_{c}^{-1}$ in Fig. 4). Ions with $A/Z=4$ have just started overrunning the shock, but there are only few particles with $p_{x}<0$ in the upstream, implying that DSA has not yet been established. Finally, ions with $A/Z=8$ isotropize too far downstream to overrun the shock barrier and are not injected yet.

At later times ($t=700\omega_{c}^{-1}$, right column in Fig. 4), instead, all of the species show the typical DSA spectrum comprising non-thermal particles that diffuse on both sides of the shock. From the color code it is also possible to see how the fraction of particles that leak back into the upstream is larger for heavier ions.

Proton injection is controlled by the quasi-periodic reformation of the shock barrier (Caprioli et al., 2015); instead, injection of heavier ions relies on rapid electromagnetic fluctuations larger than those induced by the local shock reformation and happens at later times for heavier species, and always after the onset of non-linear turbulence 333The expected dependence $E_{max,i}(t)\propto Z$ is achieved only when the acceleration time becomes much larger than the injection time, which may explain why in some heliospheric shocks Fe spectra roll over at significantly lower energies compared with O (Desai et al., 2016)..

A more quantitative characterization of the trajectories of the ions that get injected into DSA is beyond the scope of this Letter. Note, however, that the behavior reported here is not equivalent to the so-called thermal leakage scenario for particle injection (e.g., Edmiston et al., 1982; Malkov, 1998; Blasi et al., 2005), in that the injected ions are not those in the tail of the Maxwellian (strictly speaking, they have not yet thermalized). The global shock structure is always controlled by species with the most inertia, so that density and fields jump within one gyroradius of thermal protons. The isotropization length for heavy ions is effectively larger than for protons, but injection is controlled by how rapidly they can be isotropized (i.e., reverse their $p_{x}$), which depends on the local electromagnetic fluctuations and not on the ion energy.

We conclude that, while the injected protons are reflected by the shock barrier and need to be pre-energized via few cycles of shock-drift acceleration (Caprioli et al., 2015), heavy ions reflect off post-shock magnetic irregularities. The enhancement in ions with $A/Z\gg 1$ is then due to the fact that they are not affected by the proton-regulated shock barrier (their kinetic energy being much larger than the barrier potential), so that they do not experience shock-drift acceleration but rather start diffusing right away. Ions with $A/Z\gtrsim 1$ exhibit intermediate properties between protons and heavier ions because their probability of being reflected or transmitted at the shock barrier depends on the actual angle between their momentum and the shock normal: ions with velocity mainly along the shock surface are more proton-like because they can be reflected by the shock barrier (Caprioli et al., 2015).

Conclusions.— We have presented the first ab-initio calculation of ion DSA at non-relativistic shocks, finding that species with large $A/Z$ show enhanced non-thermal tails with respect to protons, in quantitative agreement with the chemical abundances observed in GCRs. In forthcoming publications we will discuss the implications of these findings also for what concerns the discrepant hardening of non-H species in GCRs (e.g., Ahn et al. [CREAM Collaboration], 2010) and for the role of accelerated He in SNR shocks (Caprioli et al., 2011).