Chemistry of ion injection in supernova remnant shocks: hybrid simulations in the light of He/C/O data from AMS-02

The high precision spectrometry of galactic cosmic rays (CR), e.g. PAMELA experiment [1], accurately determined an approximately 0.1 difference between the rigidity spectral indices of protons and He ions. Similar deviations have been indicated earlier by other experiments [2] confirmed by the recent high-fidelity AMS-02 measurements [3]. These findings widen the window into a complicated selection process whereby collisionless shocks, such as supernova remnant (SNR) blast waves, extract different chemical elements from an interstellar matter mix. The similarity of He/p, C/p, and O/p rigidity spectra demonstrated by AMS-02 has provided new evidence that injection is a mass-to-charge (A/Z) dependent process. We investigate the particle injection into the diffusive shock acceleration using one- and two-dimensional self-consistent hybrid simulations. Our 1D simulations prove that an SNR shock can modify the chemical composition of accelerated CRs by preferentially extracting them from a homogeneous background plasma without additional, largely untestable assumptions. Our results confirm the earlier theoretical predictions of how the efficiency of injection depends on the shock Mach number MA. They show that selection rate of different ion species increases with A/Z, saturates, and peaks as a function of MA. The 2D simulations also evidence a deviation from the linear selection rate vs A/Z trend, pointing towards a saturation for higher A/Z. The integrated SNR rigidity spectrum for proton-to-helium flux ratio, obtained by the convolution of the time-dependent injection rates of protons and He ions with a decreasing shock strength over the active lives of SNRs, compares well with the AMS-02 and PAMELA data. The numerical results [4] correctly predict the decrease in p/He flux ratio with increasing rigidity at exactly the rate measured in the experiments for ℜ > 10 GV. Only at lower rigidities, ℜ < 10 GV, the difference between the data and our predictions becomes noticeable. Whether this deviation comes from the propagation through the interstellar medium or solar modulation, remains unclear. Except for this uncertainty, the suggested mechanism for A/Z-dependence of the injection fully explains the measured p/He ratio.


Introduction
The difference between rigidity spectral indices of protons and He ions 5 , first measured by the balloon-born experiment ATIC-2 [2], is now thoroughly vetted by the high-fidelity, more than 115-billion cosmic ray AMS-02 database [3], Fig.1 (left frame). These findings are perceived to undermine the leading hypothesis of galactic cosmic ray (CR) origin, whereby they are accelerated in supernova remnant shocks. The specific mechanism behind the acceleration, known as the first order Fermi acceleration (or diffusive shock acceleration, DSA), is electromagnetic in Nature. It has been studied for more than four decades and its fundamental aspect regarding the particle composition is the exclusive rigidity dependence of their spectra. In fact, following [5], consider the equations of motion of a charge particle in arbitrary electromagnetic fields E( r, t), B( r, t) d R dct = E( r, t) + R × B( r, t) with R and r being particle rigidity and position correspondingly and R 0 = (A/Ze)m p c 2 . It follows from Eq.(1) that all particles with R R 0 independent on the particles species have the same orbits in the phase space ( r, R). In other words, protons and He ions with the same rigidity should be accelerated alike, given the shock compression. For this reason, most of the explanations for the "paradox" appeal to SNR environmental factors, such as inhomogeneous p/He mixes in the shock upstream medium, variable ionization states of He, or even a multi-SNR origin of the observed spectra. The question now is whether such special conditions are vital for the explanation or an SNR shock can modify the chemical composition of accelerated CRs by preferentially extracting them from a homogeneous background plasma without additional and untestable assumptions. It has been argued [5] that such an intrinsic, mass-to-charge ratio (A/Z) based, elemental selectivity should indeed operate at the initial (so-called injection) phase of the Fermi acceleration. The coincidence in particle spectral slopes of three different elements (He, C, and O) discovered recently by the AMS-02 [6] rules out incidental selection mechanisms and points to the A/Z-based one.
In this paper we investigate the particle injection into the DSA using self-consistent hybrid simulations. We provide sufficiently detailed Mach number scans of the p/He injection ratio, that is needed to test the injection bias. Reducing the spatial dimensionality allowed us to perform the modeling with adequate particle statistics, which is vital, when exploring the high-energy tails of the distribution functions, and high grid resolution. The possibility of self-consistent treatment of He ions without need to dramatically reduce the He abundance to keep the noise level low, made it possible to study the more realistic plasma composition. It should be noted, that when the species abundance is relatively high (10% of He ions in a number density is indeed a high abundance) the test-particle approximation, well justified for ions with very low abundance, is probably not the best choice, provided the alternative treatment is possible. Indeed, our simulations with He dynamics included self-consistently show the difference in particle spectra not only for the He ions, but also for protons.

Hybrid simulations
The fully consistent description of the nonlinear shock-wave structure belongs to one of the most challenging problems in numerical physics. This is because the most important and interesting phenomena are multiscale in nature and a hydrodynamic treatment is generally not applicable. While the "ab initio" full electromagnetic particle-in-cell (PIC) codes provide the most reliable numerical description of plasmas, they are computationally very expensive. Following the evolution of a collisionless shock over many ion cyclotron times and with realistic electron-ion mass ratio is even not feasible by means of PIC simulations. If the focus is on the time and spatial scales determined by the ions and the electron scales do not need to be resolved, hybrid modeling (see [7] and references therein) has been proven to be a valuable tool.
Within the hybrid approach, the evolution of the ion distribution function f ( x, v) is governed by the kinetic Vlasov equation: where E, B, and J are the electric and magnetic fields, and the current density, q i = Z e and m i = A m p denote the ion charge and mass (e and m p are the proton charge and mass, respectively). The charge neutralizing electronic plasma component in turn obeys the MHD equations. We adopt a massless electron approximation, in which the equation of motion of the electron fluid is reduced to Here −e, m e , n e , and v e are the electron charge, mass, number density and bulk velocity. In the following we assume that the pressure is isotropic and can be calculated using an adiabatic equation of state, p e ∝ n γ e , with adiabatic index γ = 5/3. We allow for the resistive coupling between electron and ions by introducing a phenomenological anomalous resistivity η. It gives rise to electron Ohmic heating and smoothes the fields on the resistive scale-length. A detailed study of the influence of η was performed in [8], where it was found in particular, that the structure of the shock does not change as long as η is sufficiently small, so that the density and magnetic field remain well coupled. The fluid equation (3) can be solved for the electric field under the assumption of quasi-neutrality n e = n i = n, with the electron bulk velocity from For the solution of the equation (2) for the evolution of the distribution function of the ion plasma-component the PIC method is used. In this method the distribution function f is sampled by a large number N max of macro-particles, whose motion is governed by the following non-relativistic equations 6 : The plasma density n = n i and ion current density J i are reconstructed from ion positions and velocities 6 As we are interested in the injection into the DSA | v| c holds during the simulation.
The equations for the electric and magnetic fields are solved using finite differences on the Eulerian grid in space. We use a first order weighting to obtain the fields at the particle position as well as when depositing ion density and current to the grid. For the numerical calculation equations (5) are discretized and centered in time here ∆t is the time-step. The particle positions are updated straightforwardly the Boris algorithm [9] is applied to update their velocities. For the solution of the Maxwell equations a magnetostatic (Darwin) model is employed, in which the displacement current is neglected. This inhibits the propagation of high frequency waves. The magnetic field evolves according to Faraday's law, The equations for the electric field (4) as well as for the propagation of the magnetic field (11) are discretized using second order finite difference stencils. The electric field is advanced from the time level (m) to (m + 1) in following steps: first E is evaluated at (m + 1/2) using the discretized Eq. (4) with B m+1/2 calculated from the discretized Eq. (11): B m+1/2 = B m − c ∆t 2 (∇ × E m ). Then E m+1/2 is used to calculate B m+1 . Finally the electric field value E m+1 at time level (m + 1) is to be obtained from E m+1 = F ( B m+1 , n m+1 i , J i m+1 ). This last step is not straightforward since is requires the knowledge about the ion current at time level (m + 1). Several methods have been proposed in the literature (see e.g. [10] for a review) to overcome the apparent difficulty of solving E m+1 = F ( B m+1 , n m+1 i , J i m+1 ) with an unknown J i m+1 .
In our work we used two of them: the predictor-corrector method and the Bashford-Adams extrapolation of the ion current. In the predictor-corrector method the "predicted" electric field is used to propagate the ions and obtain the source terms n m+3/2 i and J i m+3/2 needed for the calculation of E n+3/2 . The corrected field at time level (m + 1) is obtained then from This algorithm is simple and has good energy conserving properties but it suffers from higher computational costs, since the particles have to be propagated twice per time step. Hence, we also implemented a Bashford-Adams extrapolation of the ion current in our code, It is used by default when performing two-dimensional simulations, the magnetic field in this case is updated with a fourth order Runge-Kutta scheme.
A simulation is initialized by sending a super-alfvénic (and super-sonic) flow of plasma with mean velocity v 0 = (−v 0 , 0, 0) against a reflecting wall, placed at x = 0, Fig.1. A shock wave forms due to the interaction of the counter-propagating streams and is then traveling in x-direction in our simulation set-up. In order to increase the particle statistics the spatial dimension is reduced to one E = E(x), B = B(x) or two E = E(x, y), B = B(x, y) but all components of the velocity and fields are included. In the two-dimensional case periodic boundary conditions are applied in the y-direction. We measure time in units of inverse proton gyrofrequency, ω −1 c = (e B 0 /m p c) −1 , where B 0 is the magnitude of the background magnetic field. Lengths are normalized to the ion inertial length, c/ω p , with ω p = 4π n 0 e 2 /m p being the proton plasma frequency, where n 0 is the upstream density. Finally, the velocity is given in units of the Alfvén velocity, v A = B 0 / 4π n 0 m p . In all simulations presented below the background magnetic field is set to be (quasi-)parallel to the shock normal. The electron fluid is assumed to be initially in thermal equilibrium with the ions with β e = β p = 1.

Simulation results
To investigate the mass-to-charge dependence of the ion injection into the DSA we have performed a series of one-dimensional simulations with large simulation boxes with lengths up to L x = 48 × 10 3 c/ω p and the spatial resolution of ∆x = 0.25 c/ω p . The time step was ∆t = 0.01 Besides protons also helium, carbon and oxygen ions were included in the simulations with charge states Z = 1 and Z = 2. We used 100 particles per each ionspecies per numerical cell 7 . All ion species were treated self-consistently with relative abundances ∼ 10% for helium and ∼ 0.04% for carbon and oxygen. Note, that the self-consistent treatment is especially important for helium, which contributes with approximately 10%, and hence influences significantly the dynamics of the waves excited in the upstream medium. The simulations were performed on 360 ÷ 480 cores of the HLRN Cluster Module (Berlin, Germany). The largest 1D simulation took approximately 40 hours.
We accelerated and in the downstream energy spectra a power-law tail becomes clearly visible, Fig.2 (left frame). After the power-law tail has been certainly formed we compute the number injection efficiency (or selection rate), i.e. fraction of particles in the tail of the distribution function f α (E) for each ion species α. The resulting dependence of the number efficiency on A/Z is depicted in Fig. 2 (right frame) for simulations with far-upstream flow velocities v 0 = 5, 7, 10 v A , which correspond to the (alfvénic) shock Mach numbers M A = v shock /v A = 6.8, 8, and M A = 13. Different colors indicate the times at which the number efficiency is obtained. In all cases the fraction of accelerated particles increases as function of A/Z, saturates and then decreases, recovering a physically correct A/Z → ∞ asymptotic behavior expected for neutral particles.  [4]) Proton-to-helium ratio as a function of particle rigidity: simulation results (brown line), data from PAMELA and AMS-02 (shaded areas), theoretical fit [5] (orange line). The measured p/He ratio is accurately reproduced in the range R > 10 GV.
The position of the maximum depends in an increasing manner on the shock Mach number M A [4]. Thus the observed η sel (M A ) behavior of the selection rate confirms the earlier theoretical predictions [5]. By performing a series of simulations for the two ion species plasma (protons and He ions) we obtained a sufficiently detailed scan of the the injection efficiency for both species over a range of Mach numbers. By combining in [4] the injection efficiency from the simulation with DSA predicted power-law and accounting for the evolution of a SNR during the Sedov-Taylor stage, we calculated the proton-tohelium ratio, extending in energy beyond the capability of any simulation. The resulting rigidity spectra are presented in Fig. 3  description of rigidity spectra reconstruction procedure can be found in [4]. Most importantly, we found that p/He decreases with rigidity for R > 10 GV with almost the same rate ∆q ≈ 0.1 as observed by the high accuracy measurements 8 . A deviation from the measurements is visible at lower rigidities, R R 0 = (A/Ze)m p c 2 , where the equations of motion are different for p and He. The most likely cause of this deviation is particle interaction with the turbulent solar wind in the Heliosphere, but the interaction with the interstellar medium turbulence may also contribute. By contrast, the deviation from the AMS-02 data in the high-rigidity range where the equations of motion for p and He become identical, is insignificant, but whether it comes from a simplified integration over the SNR or it is a mixing effect from different SNRs or spallation in the interstellar medium, remains unclear. Now, to confirm the 1D prediction [4] about the A/Z behavior of η sel we have performed 2D Cartesian simulations with v 0 = 10 v A , including different ion species with A/Z up to 8. The background magnetic field is initialized to form a θ Bn = 20 • angle with the shock normal, i.e., the shock is quasi-parallel. Note, that in the 1D set-up θ Bn = 0 • was used. To follow the shock for a long time, we use a very elongated box with dimensions L x × L y = 20000 × 20 (c/ω p ) 2 and a grid spacing of ∆x = ∆y = 0.5 c/ω p . We have checked that the results do no change significantly with the transverse box size for L y ≥ 20 c/ω p . The simulations were performed on 480 cores of the HLRN Cluster (Berlin, Germany) and took approximately 42 hours.
In the top panel of Fig. 4 the downstream energy spectra obtained at t = 1000 ω −1 c are shown. The energy is measured in terms of E sh = m p v 2 0 /2. The Maxwellian part at low energies (the dotted line in Fig. 4 depicts a Maxwellian distribution) as well as the power-law tails at high energies are clearly visible. Following [4] we compute the fraction of energetic particles in the tail of the distribution function. The results for ion species with charge state Z = 1 are plotted in the lower panel of Fig. 4. They confirm the prediction of the 1D modelling: for small massto-charge values the selection rate increases with A/Z, but deviates from the linear trend for A/Z > 8, pointing toward a saturation for higher mass-to-charge ratios. The decrease of the selection rate at even higher A/Z could not be observed in our 2D simulations, because only ions with mass-to-charge ratios up to 8 were included.

Summary
We have investigated the ion injection into the DSA by means of 1D and 2D self-consistent hybrid simulations. Our results confirm the earlier theoretical predictions as well as the results of numerical modeling [12], that the injection depends on the mass-to-charge ratio as well as on the shock Mach number. Our 1D simulations show that a selection rate increases with A/Z, saturates, and peaks as a function of shock Mach number. The question of the exact position of