© ESO 2020

1. Introduction

Collisionless shocks are ubiquitous in astrophysical environments and are proven to be efficient particle accelerators, with supernova remnant (SNR) shocks being the most probable source of galactic cosmic rays (Gaisser 1991). The acceleration of charged particles at these shocks is accurately described by the theory of diffusive shock acceleration (DSA, for a recent review, see (Schure & Bell 2012). While this mechanism is conceptually simple, its precise outcome for the energy spectra and chemical composition of accelerated particles is difficult to determine under realistic shock conditions (see e.g., Ohira et al. 2016; Hanusch et al. 2019b). In particular, the “injection” of different species into the DSA remains largely unsolved (Caprioli et al. 2017; Hanusch et al. 2019a; Evoli et al. 2019), with the electron injection being a notoriously difficult problem (see e.g., West et al. 2017; Aharonian et al. 2017). An important question is whether protons that are accelerated to high energies create the sufficient wave turbulence required for the electron injection into the DSA. Electrons, in contrast to ions, have long been thought to be incapable of injecting themselves on self-generated waves, especially in quasi-parallel shocks, so some ion “assistance” appears critical.

The acceleration of electrons at quasi-perpendicular collisionless shocks has been investigated by means of numerical simulations by a number of authors. Riquelme & Spitkovsky (2011) performed fully kinetic simulations and have shown that whistler waves are crucial for electron injection. Shimada & Hoshino (2000), Hoshino & Shimada (2002), Amano & Hoshino (2008), and Matsumoto et al. (2012, 2017) report that electrostatic waves excited by the Buneman instability and accompanied by the particles trapping are important for efficient electron acceleration. Even though Matsumoto et al. (2012) clarify electron-to-ion mass ratio dependence of Buneman instability, it is unclear, whether this mechanism is robust if a real electron-to-ion mass ratio is considered. Artificially enhanced electron-to-ion mass ratios, which are often used in fully kinetic simulations in order to obtain a converged result within finite time, may lead to physically questionable results, as discussed, for example, in Matsukiyo & Scholer (2006), and Bohdan et al. (2019a,b). An alternative approach was chosen by Burgess (2006), Guo & Giacalone (2010), and Trotta & Burgess (2018). These authors investigated the effect of the shock surface and magnetic field fluctuations on electron acceleration by following the trajectories of test-particle electrons in the fields obtained from hybrid simulations. Kinetic simulations that were performed to investigate the injection of electrons at quasi-parallel shocks (e.g., Park et al. 2015), have revealed the injection into the DSA with the scattering of both ions and electrons. Waves that are excited via Bell instability (Bell 2004) mediate the scattering process. Yet another method was chosen by Guo & Giacalone (2015), where the authors followed the trajectories of the test electrons in the prescribed, kinematically-defined electromagnetic fields in the shock region and they studied the dependence of the electron acceleration efficiency on the shock inclination and wave variance. In particular, Guo & Giacalone (2015) have shown that the acceleration of these particles does not strongly depend on the shock inclination unless the upstream turbulent magnetic field is weak. On the contrary, the electrons are found to get efficiently accelerated at quasi-parallel shocks as well, provided the upstream magnetic field fluctuations are strong enough. This finding is consistent with NASA’s earlier Cassini spacecraft observation of Saturn’s bow shock (Masters et al. 2013), which for the first time provided in situ evidence of a significant electron acceleration in a quasi-parallel shock geometry¹. The kinetic model of electron injection and acceleration at quasi-parallel supercritical collisionless shocks, which was developed by Bykov & Uvarov (1999), demonstrated that strong MHD fluctuations generated by ion kinetic instabilities are important for heating and the pre-acceleration of suprathermal electrons on very small scales.

For the analysis of observed X-ray spectra as well as in order to understand the energy partitioning between the energetic cosmic-ray (CR) and the thermal populations, electron acceleration is important along with electron thermalization. As the transition between the unshocked and shocked medium is much shorter than the collisional mean free path, an equilibration of ion and electron temperatures may only occur on long timescales. Observations of Balmer-dominated shocks have shown a dependence of the electron-to-ion temperature ratio on the shock velocity (Ghavamian et al. 2013). Fully kinetic particle-in-cell simulations of low Mach number quasi-perpendicular shocks also indicate that the electron-to-proton temperature ratio of the shocked medium is a function of the shock Mach number (Guo et al. 2018).

In this paper, we study the electron-ion temperature relaxation in a quasi-parallel shock geometry by introducing electrons as test particles in hybrid simulations and we investigate their thermalization in the proton-driven turbulence. We compute the trajectories of these electrons, but we assume that their contribution to charges and currents is negligible, so the electromagnetic fields are not affected. It should be emphasized that although hybrid simulations neglect the contribution of suprathermal electrons, assuming that their number is small enough so as to not affect the simulation results significantly, this low-density electron component, which is negligible for the simulations dynamics, is very important observationally. On the other hand, if hybrid modeling provides realistic field distributions, and there exist reasons to believe this, the study of the electron behavior in these fields is also worthwhile. However limited, this approach elucidates aspects of electron heating that are not accessible to hybrid simulations.

Our choice of a one-dimensional (1D) simulation is motivated by computational economy, spatial resolution reasoning, and physical considerations. First of all, a realistic alternative to 1D simulations would be a two-dimensional (2D) simulation. However, the 2D simulations also have drawbacks that we have recently discussed (Hanusch et al. 2019a). We argue that when the high resolution and particle statistics are the priority, a 1D code may be a better choice. Of course, when the 2D effects are crucial, as in the case of acceleration at the variable shock obliquity (Hanusch et al. 2019b), the 1D simulation setup is not a possibility.

The problems associated with the 2D option are primarily due to inverse turbulent cascades, which are not present in the real three-dimensional (3D) systems. Additionally, and more importantly for this study, the conserved particle canonical momentum component in the direction of the ignorable coordinate, in combination with strong magnetic eddies produced by the inverse cascade, result in a protracted interaction of particles with the eddies. This interaction is akin to the shock drift (or the so-called surfatron) acceleration, occurring when a particle “surfs” on the edge of an eddy. It can strongly modify the particle transport, both in momentum and coordinate space. This phenomenon has recently been studied by a direct comparison of 2D and 3D simulations by Trotta & Burgess (2018).

2. Electron preheating mechanisms

The investigation of electron energization by the electromagnetic fields generated in hybrid simulations implies the existence of independent preheating mechanisms for this electron population. Indeed, the hybrid simulations treat electrons as a fluid so that their interaction with relatively long waves that are generated by nonequilibrium ion populations (such as shock-reflected ions) within the hybrid simulations can only occur adiabatically. These waves cannot heat the electrons appreciably. However, these are not the only waves that are generated by these types of ions in the real shocks. Much shorter waves with higher frequencies, which are not accessible to hybrid simulations, may also be generated. They can tap into thermal electrons and preheat them. Then, electrons start to interact with the waves generated in the hybrid simulations.

The preheating mechanisms are not straightforward, and their thorough description is outside the scope of this short paper, while a brief overview is in order. They have been considered in many earlier publications, starting perhaps from those in the magnetic fusion research and general plasma physics, (e.g., Shapiro & Shevchenko 1968). These types of mechanisms invoke, almost universally, a combination of Cerenkov and cyclotron resonances impacting the same particle populations. The wave modes associated with these resonances may or may not be the same. The salient aspects of the wave-particle interaction are briefly explained in Fig. 1.

Fig. 1. Electron preheating in high frequency plasma waves, primarily the lower-hybrid and oblique Langmuir waves, Vph = ω/k∥ ≃ ωce/k, that are not accessible to the hybrid simulations.

Suppose that the initial population of electrons is stable, for example, a Maxwellian, but there are waves accelerating electrons in the velocity component in the magnetic field direction, v_∥, by a quasilinear diffusion. The lower-hybrid waves are, perhaps, most potent. They have an approximately constant frequency, (in the dense plasma limit ). Here ω_ce and ω_pe are the electron cyclotron and plasma frequencies, respectively. When propagating at large angles to the local magnetic field (k_⊥ ≫ k_∥), they accelerate electrons over a broad range in v_∥ via a Cerenkov resonance, ω − k_∥v_∥ = 0. However, these waves are damped by the same electrons and cannot tap into their distribution deeper than, typically, v_∥ ≃ 3V_Te because |∂f_e/∂v_∥| increases at lower v_∥ and leads to strong wave damping. Energy-wise, however, these waves, which are driven by powerful ion populations which also happen to be observed in hybrid simulations, could pull more electrons out of the core Maxwellian. They are just not excited for resonant velocities of v_∥ = ω/k_∥ <  3V_Te.

Nevertheless, there are still mechanisms whereby the low-energy electrons with v_∥ ≲ 3V_Te can also be accelerated. One such mechanism has been discussed in conjunction with the electron injection into the DSA (Galeev et al. 1995). By this mechanism, a macroscopic electric field is generated in response to the acceleration of an initially small fraction of electrons by the lower-hybrid waves in the region of v_∥ >  3V_Te, for instance. While these electrons tend to escape the spatial region of their acceleration, the emerging charge imbalance must be neutralized by the electric field in the shock precursor. This field also accelerates electrons with lower energies, placing them into the region of v_∥ >  3V_Te where they, being picked up by the waves, add up to the already heated electron population.

The second mechanism is also based on the electron distribution that is already stretched along v_∥. It may easily become unstable with respect to the wave generation via cyclotron resonance, ω − k_∥v_∥ + nω_ce = 0, where n is an integer (Shapiro & Shevchenko 1968). A quasilinear diffusion then follows along the lines on the v_∥, v_⊥ plane that are determined by the relation . Here, the wave phase velocity V_ph = ω/k_∥ along the field should be expressed using the above resonance relation. It may be seen from Fig. 1 that this diffusion leads to energy losses of electrons in the wave frame, thus making the wave unstable. An important aspect of this diffusion is that electrons are swept up to lower v_∥ and higher v_⊥, so that their distribution that is integrated over v_⊥ no longer strongly contributes to the wave damping on the Cerenkov resonance in the critical region v_∥ ≃ 3V_Te. As a result, more thermal electrons fall into the Cerenkov resonance interaction with the lower-hybrid waves and they get accelerated. To conclude this section, a significant part of thermal electrons can be made accessible to the ion-driven waves observed in hybrid simulations. Having discussed it as a proof of principle, we stop short of making specific predictions about the exact number of these types of electrons for two reasons. First, some useful information can be found in the cited papers and references therein. Nevertheless, to our knowledge, there are no such calculations broadly applicable to shocks. Otherwise, a comprehensive theory of electron injection into the DSA would have been developed. Second, a test-particle treatment of energetic electrons in this paper places limits on their number, as is discussed in the Introduction.

3. Model and simulation setup

The investigation of particle dynamics at collisionless shocks, in the context of CR acceleration at SNR shocks, has largely relied on numerical simulations. Although the in situ measurements made by spacecrafts can potentially shed light on the physics of these shocks, they are only attainable for shocks in the Solar System, such as Earth’s bow shock (by satellite missions; e.g., Sundberg et al. 2016; Feldman et al. 1983; Amano et al. 2020) or Saturn’s bow shock (Masters et al. 2013). Moreover, as Solar System shocks are much less powerful and smaller than SNR shocks, it is not clear whether the dependence of the acceleration efficiency on the magnetic field orientation with respect to the shock normal, which is observed for the electrons in Solar-wind plasma shocks, applies to SNR shocks as well.

While the use of fully kinetic simulations is the most fundamental approach, it is computationally very expensive and, in multidimensional geometry, it can even be unfeasible to follow the evolution of a collisionless shock over many ion cyclotron times. In order to overcome these difficulties, unrealistic and strongly increased electron-to-ion mass ratios are often used in kinetic simulations.

When focusing on the acceleration of ions, the hybrid approach has been proven to be a valuable tool (Lipatov 2013). In these simulations, the electrons are treated as a charge neutralizing fluid. If moreover, one neglects the electron mass, the equation of motion of the electron fluid reduces to

(1)

where −e, n_e, and v_e are the electron charge, density, and bulk velocity, respectively, and J is the total current. The last term on the right-hand side of (1) describes the resistive coupling between electrons and ions. A phenomenological anomalous resistivity η gives rise to electron Ohmic heating and smooths the fields on the resistive scale-length. Both the resistivity and pressure, p_e, are assumed to be scalar and an adiabatic equation of state with an adiabatic index of γ = 5/3 is used for the fluid electrons. In the hybrid model, the ions are treated kinetically, and their motion is governed by the following nonrelativistic equations

(2)

In the simulations, lengths are normalized to the ion skin depth, c/ω_p, with the proton plasma frequency . Here, n₀ denotes the plasma density that is far upstream and e and m_p are the proton charge and mass, respectively. Time is measured in units of inverse proton gyrofrequency, , and velocity is measured in units of the Alfvén velocity, . Here, B₀ denotes the magnitude of the background magnetic field, which is set to B₀ = B₀ e_x, which is parallel to the shock normal (x-axis in our convention).

In the following, we investigate the electron kinetics in the fields generated by the ions, which is beyond Eq. (1). If stochastic fields generated in hybrid simulations are realistic, electron orbits in these fields also deserve attention. To this end, we add a population of electrons as test particles in our simulation. These electrons, by definition, do not generate electric or magnetic fields, nor do they exert pressure on the background plasma. The idea is in the spirit of earlier simulations of quasi-perpendicular shocks, where the fields obtained from hybrid simulations were used to propagate test-particle electrons (Burgess 2006; Guo & Giacalone 2010; Trotta & Burgess 2018).

Due to the separation of scales, we have introduced a sub-cycling routine in our hybrid numerical code to properly resolve the trajectories of the test-particle electrons. The propagation of the electrons is performed in N_cyc sub steps, reducing the effective time step for the electron propagation to Δt_e = Δt/N_cyc. A linear interpolation between the fields known at the time-steps of the ion propagation t = n Δt and t = (n + 1) Δt is used to obtain the fields at the sub steps. To reduce the numerical costs, we use a moderately increased electron-to-proton mass-ratio of m_e/m_p = 1/400 and update the electron positions and velocities N_cyc = 20 times during one propagation step Δt of the ions. We note that the use of the guiding center approximation or a gyro-kinetic treatment (Frieman & Chen 1982; Littlejohn 1983) of the electrons would also be possible; however, in starting from some energy, adiabaticity can stop working well for the electrons. The electron fluid is initially assumed to be in thermal equilibrium with the ions with β_e = β_p = 1. The simulation is initialized by sending a super-sonic and superalfvénic hydrogen plasma flow with velocity v₀ against a reflecting wall, placed at x = 0, as is seen in Fig. 2. A shock forms upon the interaction of the counter-propagating plasma streams and propagates in the positive x-direction. Since we are not sure about the heating mechanism, we considered different far-upstream distributions for the test-particle electrons as follows: (T) – the mean velocity of the test-electron population equals the plasma flow upstream speed of v_0e = v₀ and its temperature equals the temperature of the electron core distribution on which the hybrid simulation operates; (B) “beam” – the mean velocity of the test-electron population is significantly bigger than the far upstream plasma flow speed of v_0e ≫ v₀ and v_0e ‖ v₀, and its temperature equals the temperature of the electron fluid; (S) “shell” – the energy of test electrons is significantly higher than the directional energy of the electrons in the plasma flow v_0e = 100 v_A  ≃ 5V_Te ≫ v₀. The test-particle electrons are injected at the right boundary starting from t = 0 (T) and from tω_c = 50 for v_0e/v_A = 30 and v_0e/v_A = 100, (B), respectively, with a velocity distribution according to a Maxwellian flux with a drift velocity of v_0e (see Cartwright et al. (2000) for details). The simulations were performed with a temporal resolution of Δt = 0.01(c/ω_p)/v₀ and a cell size of Δx = 0.25 c/ω_p. At least 1000 ions per cell were used to keep the numerical noise² in the electromagnetic hybrid fields at a low level. To strongly improve particle statistics and avoid problems inherent in 2D simulations, as is briefly discussed in the Introduction, only one spatial dimension (x in our convention), but all of the components of the velocity and fields, are included.

Fig. 2. Simulation setup: a shock is created by sending a super-sonic plasma flow with a velocity v0 against a reflecting wall. The shock propagates to the right, parallel to the background magnetic field.

4. Results

We have performed hybrid simulations for different initial upstream flow velocities, giving rise to the formation of the shocks with different Mach numbers, and we followed the evolution of the shock for several hundreds of ion cyclotron times. In the following, we present the results of the simulations investigating the behavior of test-particle electrons moving in the turbulent fields that are created by the ion plasma component.

Figure 3a) shows the temporal evolution of the spatially dependent B_y for a simulation with v₀/v_A = 10. The shock propagates to the right and the compression of the magnetic field is clearly visible. The dashed line denotes the time at which the density and components of the magnetic fields are plotted in Fig. 3b). It can be seen that the ion density increases upon the shock crossing. Circularly polarized Alfvén waves are excited by the streaming protons. As the magnetic field of these waves is almost frozen into the plasma (v_A ≪ v₀), they are advected downstream and compressed, leading to large amplitudes of the magnetic field behind the shock front.

Fig. 3. (a) Space-time plot of the y-component of the magnetic field. (b) Ion density ni and the components of the magnetic field, By and Bz, for a simulation with an upstream flow velocity of v0/vA = 10 at t ωc = 300.

Detailed information about the accelerated particles and their temperature can be extracted from the particle phase spaces. Figure 4 shows the distributions of protons (top) on the (x, v_x) plane and of the test electron (bottom) on the (x, v_∥) plane at t ω_c = 400. Here, v_∥ is the electron velocity component parallel to the local magnetic field. In the proton phase, space accelerated particles are visible upstream and downstream. A large increase in the proton temperature can be inferred from the width of the proton distribution. For the test-particle electrons, the initial Maxwellian flux in the far upstream widens in the precursor (between x = 1000 c/ω_p and x = 2000 c/ω_p), indicating an increase in temperature. Furthermore, a population of counter-propagating particles with positive v_∥ is present upstream. This indicates a reflection due to magnetic mirroring near the shock transition. The inset shows a region close to the shock transition, where this effect is more pronounced. Upon a shock crossing, the width of the distribution in velocity changes only slightly and only a minor increase in the temperature of the test-particle electron population is expected.

Fig. 4. Phase-space fp(x, vx) of protons (top) and fe(x, v∥) of test-particle electrons (bottom) at t ωc = 400. The inset shows fe(x, v∥) to be close to the shock transition on a linear scale for a simulation (T) with v0/vA = 10.

The multiple reflection of electrons from the shock may significantly increase their energy. The number of these electrons, experiencing repeated reflections at the shock front and getting scattered by the upstream ion-generated turbulence, is considerable, as suggested by the histogram of the test-electron reflection in the time interval 275 <  t ω_c <  300, which is shown in Fig. 5 (top frame). With an increasing number of reflections, the velocity of the test particles increases on average as well; this can be seen in Fig. 5 (bottom frame). We have traced some reflected particles and plotted the results in Fig. 6. The background shows the amplitude of the magnetic field |B(x, t)|. The line color denotes the energy of the traced electrons in terms of their initial energy at t ω_c = 275. It is clearly visible that all the traced test particles gain energy upon reflection and in the interaction with the proton-driven turbulence in the shock precursor³. This can also be inferred from the phase space plot, shown in Fig. 4, where particles with large positive v_∥ are present. This process might be important for the injection of electrons into the DSA, though the downstream energy spectrum of the test-particle electrons does not show a clear power-law tail. We note that injection is not necessarily equivalent to the full DSA acceleration and the goal of our paper is to show the possibility of the first, most critical part of the two. Moreover, a clear power-law tail is not a prerequisite for injection. In the DSA context, the electron injection is understood to be an emergence of an electron population that can potentially be scattered by the ion-generated waves. In Fig. 7, we follow an exemplary test electron which, after being reflected off of the shock because of magnetic mirroring, remains trapped for a while upstream, close to the shock front. The relatively high velocity of this electron at tω_c = 275 indicates that it has already been energized in the proton-driven turbulence developed in the precursor. The results of a recent study by Guo & Giacalone (2015), where the test-particle electrons were propagated in a predefined turbulent electromagnetic field in the shock region, suggest that electrons return back to the shock due to mirroring reflection off of the large amplitude waves upstream. The exemplary trajectory, which is shown in Fig. 7, indicates that there is indeed a mirroring process upstream. The wave spectra upstream of the shock transition, seen in Fig. 8, show that strong forward propagating waves, with the group and phase velocities that are higher than the Alfvén velocity with 0.05 <  k_x ⋅ (c/ω_p) < 0.2 and 0.2 <  ω/ω_c <  1, are excited so that the resonance condition ω − k_xv_x = 0 can be fulfilled for the electrons having v_x ∼ 20v_A.

Fig. 5. Top: reflection on the shock histogram in the time interval 275 <  t ωc <  300. Bottom: velocity of the reflected electrons after the last reflection in the time interval 275 <  t ωc <  300 for a simulation (T) with v0/vA = 10. The color denotes the fraction of reflected test electrons.

Fig. 6. Trajectories of the electrons that are reflected off of the shock in the magnetic field |B(x, t)| for a simulation (T) with v0/vA = 10. The line color denotes the energy of the test-particle electrons in terms of their initial energy at the start of the tracing at tωc = 275.

Fig. 7. Trajectory of an electron in the (x, t) and (x − xsh, v) spaces in the time interval 275 <  t ωc <  300 for a simulation (T) with v0/vA = 10. The background (left frame) shows the amplitude of the magnetic field |B(x, t)|. The test-electron thermal velocity that is far upstream is The relatively high velocity of the traced particle (v/vA ∼ 200) at t ωc = 275 indicates that it has already been energized in the precursor and possibly during encounters with the shock at earlier times.

Fig. 8. Spectrum of the magnetic field upstream for the simulation with v0/vA = 5.

The downstream distribution f_e(v_∥) at t ω_c = 400 for a simulation (T) with v₀/v_A = 10 is plotted in Fig. 9. This is a low-density addition to the electron core distribution on which the hybrid simulation operates. Since both the high temperature part, which is described by a Maxwellian⁴ well (dashed line), and the energetic tail are absent in the core, this correction, however small in density, is physically important. We emphasize that the tail-to-core content in the test-particle electron population is vastly different from that of the fluid electrons. Behind the shock transition, the relative amplitude f^TP(150 v_A)/f^TP(0) of the tail of the test-electron distribution at v = 150 v_A is approximately 0.05; whereas, for the core electron distribution, it takes the value f^fluid(150 v_A)/f^fluid(0)≃9 × 10⁻⁶. We note that the velocity distribution of the suprathermal tail population fits with the Kappa (κ-) distribution best

Fig. 9. Velocity distribution of fe(v∥) for the test-particle electron population downstream of the shock transition at t ωc = 400 for a simulation (T) with v0/vA = 10. The center of the distribution is described by a Maxwellian well (dashed-line). A suprathermal tail is clearly visible. The relative amplitude of the tail at v = 150 vA is significantly higher, approximately 5500 times, than the value obtained for the Maxwellian core electron distribution (electron fluid).

which is shown by a dotted line in Fig. 9. Here, and Γ is the Gamma function. Kappa distributions are frequently employed to describe the velocity distribution of collisionless plasmas out of thermal equilibrium (Lazar et al. 2016). These include space and astrophysical plasmas from solar wind and planetary magnetospheres to the heliosheath and beyond to interstellar and intergalactic plasmas (see e.g., Feldman et al. 1983; Pierrard & Lazar 2010; Oka et al. 2018; Livadiotis et al. 2018; Wilson et al. 2019; Livadiotis 2017, and references thererein).

The spatially dependent temperature profiles of protons (green), test-particle electrons (orange), and electron fluid (blue) are plotted in Fig. 10 for a simulation (T) with an upstream flow velocity of v₀/v_A = 10. It is apparent that the temperature profile of the test particles differs from the temperature of the electron fluid, which in our hybrid model have to follow T_e(x)/T₀ = (n_e(x)/n₀)^γ − 1 because of the adiabatic closure between the density and pressure. The downstream temperature of the test-particle electrons is approximately two times higher than the fluid temperature. A gradual increase in the temperature of the test particles toward the shock transition in the region between x = 1200 c/ω_p and x = 2600 c/ω_p points to a considerable potential of the precursor wave field, which is supplied by ions, to preheat the electrons before they are shocked at the subshock⁵. Independent of the initial test-electron mean energy (distributions (T), (B), or (S)), the obtained temperature space profiles confirm the heating of the test electrons in the precursor. We emphasize here that this considerably heated test-particle electron population, by definition, does not generate electric or magnetic fields, nor exerts additional pressure on the background plasma.

Fig. 10. Spatial dependence of the temperature for the electron fluid (blue), test-particle electrons (orange), and protons (green) at t ωc = 400 for a simulation (T) with v0/vA = 10. Independent of the initial test-electron distribution, (T), (B), or (S), their temperature space profiles show significant heating in the precursor.

To obtain the dependence of the suprathermal (with respect to the core-electron population) electron temperature on the shock Mach number, M_A = v_sh/v_A, we performed simulations for a range of initial flow velocities, v₀ (i.e., different shock velocities, v_sh), and we calculated the downstream test-electron temperatures for different initial test-electron distributions (T), (B), and (S). The resulting proton, T_i, and test-electron temperatures, , in terms of the far-upstream plasma temperature, T₀, are summarized in the Table 1. We note that the downstream electron velocity distribution for all the simulations (T) and for the simulations (B) and (S) and v₀/v_A ≥ 10, is properly described by a Maxwellian, see Fig. 9, and the temperature is extracted from the Maxwellian fit. For simulations (B) – “beam” and (S) – “shell” and relatively low shock Mach numbers (v₀/v_A = 5, 7 corresponding to M_A = 6.9, 9.6), the downstream test-electron distribution function is “shell-like”, as demonstrated in Fig. 11, so that the “thermal” velocity of the test electrons (indicated by a ^* symbol in Table 1) is approximately equal to the radius of the downstream velocity “shell”. Additionally, we calculated the velocity distribution for test electrons that reside downstream for a certain time: t_dsω_c = 175 for v₀/v_A = 5 and t_dsω_c = 100 for v₀/v_A = 7. These electrons had time to thermalize, and their distribution is close to a Maxwellian. The numbers in brackets in the Table 1 refer to the temperatures that were obtained accordingly.

Fig. 11. Downstream velocity distribution of the test electrons f(v∥, v⊥) along field line for a simulation (B) with v0/vA = 5 averaged over the whole downstream region.

The results, which are displayed in Fig. 12 with blue diamonds, provide evidence that the test-electron-to-proton temperature ratio, , is a decreasing function of the shock Mach number with a tendency for saturation at high M_A. For comparison purposes, in the same graph, we show the electron-to-proton temperature ratios that were determined from observations of Balmer-dominated shocks (orange symbols) as a function of the shock velocity⁶v_sh (Ghavamian et al. 2013, 2007; van Adelsberg et al. 2008). The values of v_sh and the temperature ratio T_e/T_i were extracted from the width of the broad and narrow components of the Hα line profile (see e.g., Fig. 1 in Ghavamian et al. (2013)).

Fig. 12. Test-electron-to-proton temperature ratio extracted from the simulations (for the upstream distribution (T) – blue filled diamonds – three time moments for every MA are shown) together with the temperature ratios, (Te/Ti)(vsh), determined from observations of Balmer-dominated shocks (Ghavamian et al. 2013, 2007; van Adelsberg et al. 2008) (orange).

The observational data, as can be seen in Fig. 12, show that T_e/T_i decreases with increasing shock velocity; additionally, at high v_sh, it seems to saturate on a level which is higher than the mass proportional heating, expected from the Rankine-Hugoniot jump conditions For an analysis of the observations, Ghavamian et al. (2007) applied a model of electron heating in which a constant level of electron heating over a wide range of shock velocities (see also Bykov & Uvarov 1999) is implied. It is also assumed that the ion heating increases with the shock speed. A scaling was found to fit to the observational data best (Ghavamian et al. 2013) (orange-dashed line in Fig. 12). Instead, a function is best fitted to our simulation results for the suprathermal electrons. A relatively high Alfvén velocity of v_A = 90 km s⁻¹ has to be assumed to make a comparison of the measurements with the simulation data. This is about four times larger than the velocity one would expect when using the standard parameters of the interstellar medium (for B = 3 μG and n​ = ​0.1 cm⁻³, the Alfvén velocity equals v_A ≃ 20 km s⁻¹). However, SNR environments are diverse and if, in addition, the large-scale field is amplified, one might have to consider the Alfvén velocity in the amplified field, which can be as high as δB/B = 4−10.

5. Discussion and conclusions

Our simulations indicate that an equilibration of test electron and ion temperatures does not occur. A decrease in the shock Mach number electron-to-ion temperature ratio, , is observed instead. We note that our test electrons correspond to the suprathermal, with respect to the core-electron population, which is part of the distribution function and the simulations do not tell how abundant this electron population might be. The shock velocity dependence of the ratios of the core electron-to-proton temperature, (T_e/T_i)(v_sh), that were experimentally determined for Balmer-dominated shocks has a similar trend. For strong shocks, the observational data suggest a saturation or even an upturn of the temperature ratio. The observed scaling, , is supported by theoretical predictions (Vink et al. 2015) that a dependence of , where M_s = v_sh/c_s is the sonic Mach number, can be obtained by solving the Rankine-Hugoniot jump conditions with the assumption that the enthalpy flux is conserved for each particle species separately. This yields an behavior for shocks with Mach numbers in the range

(3)

(see Eq. (14) in Vink et al. 2015) with μ = (m_i + m_e)/2 being the average mass and r being the shock compression ratio. For the increased electron-to-ion mass ratio and parameters used in our simulations, Eq. (3) translates to 1.4 <  M_A <  23. In the case of efficient CR acceleration, that is, when considering the CR pressure in the precursor, Vink et al. (2015) also predict a higher level of saturation of the temperature ratio toward a high M_s⁷. The predicted range for the scaling, with a dependence on m_e and m_i, disfavors fully kinetic simulations, where the reduced mass ratios of m_p/m_e = 64−100 are used regularly (Park et al. 2015). In this case, scaling may only occur in a limited range of Mach numbers.

To conclude, different models have been proposed to heat electrons in front of SNR shocks. All of these models have to rely on numerical simulations since the particle distributions in SNR shocks cannot be measured in situ. When focusing on nonrelativistic collisionless shocks, two scenarios have mainly been considered: heating due to lower hybrid waves in the precursor (Laming 2000) or a mechanism based on counterstreaming instabilities in front of the shock (Cargill & Papadopoulos 1988). Both mechanisms work well for perpendicular shocks. In Rakowski et al. (2008) and Malkov et al. (1995), it has been argued that due to the amplification of the magnetic field ahead of the shock, the perpendicular component might be large enough for the models to also be applicable for quasi-parallel shocks. Our simulation shows that the well developed ion turbulence in the precursor is able to influence the dynamics of the electron population and is responsible for the electron preheating in the quasi-parallel shocks as well.

Acknowledgments

The research was supported by DFG Grant 278305671. M.A.M. acknowledges the NASA ATP-program support within Grant 80NSSC17K0255, and the National Science Foundation under grant NSF PHY-1748958. A.H. & T.V.L. acknowledge the North-German Supercomputing Alliance (HLRN) for providing the computational resources for the simulations (project mvp00015). T.V.L. acknowledges the support by the state contract with ICMMG SB RAS (0315-2019-0009) and by the Ministry of Science and Higher Education of the Russian Federation (agreement 075-03-2020-223/3 within FSSF-2020-0018) in the part related to the analysis of the numerical results.

References

All Tables

All Figures

	Fig. 1. Electron preheating in high frequency plasma waves, primarily the lower-hybrid and oblique Langmuir waves, Vph = ω/k∥ ≃ ωce/k, that are not accessible to the hybrid simulations.
In the text

	Fig. 2. Simulation setup: a shock is created by sending a super-sonic plasma flow with a velocity v0 against a reflecting wall. The shock propagates to the right, parallel to the background magnetic field.
In the text

	Fig. 3. (a) Space-time plot of the y-component of the magnetic field. (b) Ion density ni and the components of the magnetic field, By and Bz, for a simulation with an upstream flow velocity of v0/vA = 10 at t ωc = 300.
In the text

	Fig. 4. Phase-space fp(x, vx) of protons (top) and fe(x, v∥) of test-particle electrons (bottom) at t ωc = 400. The inset shows fe(x, v∥) to be close to the shock transition on a linear scale for a simulation (T) with v0/vA = 10.
In the text

	Fig. 5. Top: reflection on the shock histogram in the time interval 275 <  t ωc <  300. Bottom: velocity of the reflected electrons after the last reflection in the time interval 275 <  t ωc <  300 for a simulation (T) with v0/vA = 10. The color denotes the fraction of reflected test electrons.
In the text

	Fig. 6. Trajectories of the electrons that are reflected off of the shock in the magnetic field |B(x, t)| for a simulation (T) with v0/vA = 10. The line color denotes the energy of the test-particle electrons in terms of their initial energy at the start of the tracing at tωc = 275.
In the text

Fig. 7. Trajectory of an electron in the (x, t) and (x − xsh, v) spaces in the time interval 275 <  t ωc <  300 for a simulation (T) with v0/vA = 10. The background (left frame) shows the amplitude of the magnetic field |B(x, t)|. The test-electron thermal velocity that is far upstream is The relatively high velocity of the traced particle (v/vA ∼ 200) at t ωc = 275 indicates that it has already been energized in the precursor and possibly during encounters with the shock at earlier times.

In the text

	Fig. 8. Spectrum of the magnetic field upstream for the simulation with v0/vA = 5.
In the text

Fig. 9. Velocity distribution of fe(v∥) for the test-particle electron population downstream of the shock transition at t ωc = 400 for a simulation (T) with v0/vA = 10. The center of the distribution is described by a Maxwellian well (dashed-line). A suprathermal tail is clearly visible. The relative amplitude of the tail at v = 150 vA is significantly higher, approximately 5500 times, than the value obtained for the Maxwellian core electron distribution (electron fluid).

In the text

	Fig. 10. Spatial dependence of the temperature for the electron fluid (blue), test-particle electrons (orange), and protons (green) at t ωc = 400 for a simulation (T) with v0/vA = 10. Independent of the initial test-electron distribution, (T), (B), or (S), their temperature space profiles show significant heating in the precursor.
In the text

	Fig. 11. Downstream velocity distribution of the test electrons f(v∥, v⊥) along field line for a simulation (B) with v0/vA = 5 averaged over the whole downstream region.
In the text

	Fig. 12. Test-electron-to-proton temperature ratio extracted from the simulations (for the upstream distribution (T) – blue filled diamonds – three time moments for every MA are shown) together with the temperature ratios, (Te/Ti)(vsh), determined from observations of Balmer-dominated shocks (Ghavamian et al. 2013, 2007; van Adelsberg et al. 2008) (orange).
In the text