Kinetic simulations of nonrelativistic perpendicular shocks of young supernova remnants. I. Electron shock-surfing acceleration

TL;DR

This study uses 2D PIC simulations to analyze electron shock-surfing acceleration at young supernova remnant shocks, revealing the importance of magnetic field orientation and simulation parameters for energization efficiency.

Contribution

It demonstrates how magnetic field orientation affects electron acceleration in shock simulations and highlights the limitations of 2D models compared to 3D.

Findings

Electron energization via Buneman waves is more efficient with out-of-plane magnetic fields.

2D simulations require parameter modifications to match 3D SSA efficiency.

Electrons gain energy primarily from electrostatic fields of Buneman waves.

Abstract

Electron injection at high Mach-number nonrelativistic perpendicular shocks is studied here for parameters that are applicable to young SNR shocks. Using high-resolution large-scale two-dimensional fully kinetic particle-in-cell (PIC) simulations and tracing individual particles we in detail analyze the shock surfing acceleration (SSA) of electrons at the leading edge of the shock foot. The central question is to what degree the process can be captured in 2D3V simulations. We find that the energy gain in SSA always arises from the electrostatic field of a Buneman wave. Electron energization is more efficient in the out-of-plane orientation of the large-scale magnetic field because both the phase speed and the amplitude of the waves are higher than for the in-plane scenario. Also, a larger number of electrons is trapped by the waves compared to the in-plane configuration. We conclude that significant modifications of the simulation parameters are needed to reach the same level of SSA efficiency as in simulations with out-of-plane magnetic field or 3D simulations.

Figures

Tables

Table 1. Table 1: Simulation Parameters

Table 2. Table 2: Dimensionless peak amplitude of Buneman waves and fraction of pre-accelerated electrons

Run	$\max \left(\left|E_{\mathrm{ES}}\right|/\left(c{B}_0\right)\right)$	$N_{e,\mathrm{BI}}/N_{e,\mathrm{tot}}(\%)$
A2	1.1	0.43
B2	1.3	0.46
C2	2.3	0.6
D2	0.4	0.34
E2	1.3	0.49
F2	1.1	0.44
G2	2.7	6.8
G2*	2.3	2.7

Equations

$$M_{\mathrm{A}}\ge (1+\alpha ){\left(\frac{m_{\mathrm{i}}}{m_{\mathrm{e}}}\right)}^{\frac{2}{3}},$$

$$M_{\mathrm{A}}\ge \sqrt{\frac{2}{1+{\sin }^2\phi }}(1+\alpha ){\left(\frac{m_{\mathrm{i}}}{m_{\mathrm{e}}}\right)}^{\frac{2}{3}},$$

$$M_{\mathrm{A}}\ge 1.5\sqrt{\frac{2}{1+{\sin }^2\phi }}(1+\alpha ){\left(\frac{m_{\mathrm{i}}}{m_{\mathrm{e}}}\right)}^{\frac{2}{3}}.$$

$$t_{\mathrm{tr},0}\approx 13{\omega }_{\mathrm{pe}}^{-1}\mathrm{and}{t}_{\mathrm{tr},90}\approx 11{\omega }_{\mathrm{pe}}^{-1},$$

$$\Delta {\gamma }_0\approx 0.18\mathrm{and}\Delta {\gamma }_{90}\approx 0.42,$$

$$\begin{array}{r}\Delta \epsilon \approx e{E}_{\mathrm{ES}}|F|\Delta v{t}_{\mathrm{tr}}=\\ =\Delta v^2{t}_{\mathrm{tr}}{\omega }_{\mathrm{pe}}{m}_e(m_e/m_i{)}^{(1/6)},\end{array}$$

$$E_x=E_{\mathrm{ES}}F\cos \Theta E_y=E_{\mathrm{ES}}F\sin \Theta ,$$

$$F=\sin \left(\frac{{\omega }_{\mathrm{pe}}}{v_{\Phi }+v_0}\left[x\cos \Theta +y\sin \Theta -v_{\Phi }t\right]+\Phi \right).$$

$${\dot{v}}_x=$$

$${\dot{v}}_y=$$

$${\dot{v}}_{x^{\prime }}=$$

$${\dot{v}}_{y^{\prime }}=$$

$$m\frac{d}{dt}\frac{(v_x+v_0{)}^2+v_y^2}{2}=e{E}_{\mathrm{ES}}|F|\left[(v_x+v_0)\cos \Theta -v_y|\sin \Theta |\right]$$

$$F=\sin \left(\frac{{\omega }_{\mathrm{pe}}}{v_{\Phi }+v_0}\left[x-v_{\Phi }t\right]+\Phi \right).$$

$${\dot{v}}_x=$$

$${\dot{v}}_z=$$

$$m\frac{d}{dt}\frac{(v_x+v_0{)}^2+v_z^2}{2}=e{E}_{\mathrm{ES}}|F|(v_x+v_0),$$

Full text

acceleration of particles, instabilities, ISM – supernova remnants, methods – numerical, plasmas, shock waves

††journal: ApJ

1 Introduction

The current paradigm of cosmic-ray (CR) origin assumes that the most part of galactic CR population is produced at nonrelativistic forward shocks of supernova remnants (SNRs). The main acceleration mechanism considered at shocks is diffusive shock acceleration (DSA), a first-order Fermi process (e.g., Axford et al., 1977; Drury, 1983; Blandford & Eichler, 1987). Astronomical observations give strong support to this paradigm. In particular, detection of broadband nonthermal emission from SNRs, extending in some objects to TeV-range gamma rays, proves the presence of ultrarelativistic particles in these sources, though for most SNRs it is still unclear which parent particle populations (protons or electrons) generate dominant high-energy emission (Aharonian, 2013).

Acceleration of particles through DSA comes from multiple interactions with the shock front, while they bounce between the shock upstream and downstream plasmas. Particle confinement to the shock vicinity is provided by elastic scattering off magnetohydrodynamic (MHD) turbulence that renders diffusive particle motions. The critical ingredient and the main unsolved problem in the DSA theory is the particle injection. CRs undergoing DSA have Larmor radii much larger than the internal shock transition width, that is commensurate with the gyroradius of the incoming protons (with shock speed $v_{\mathrm{sh}}$). CRs thus see the shock as a sharp discontinuity in the plasma flow. To be fed into the acceleration process particles need therefore to be extracted from the thermal pool and pre-accelerated. Since protons have a larger initial momentum and can be easily scattered either by MHD waves embedded in the ambient plasma or by self-generated turbulence, their injection is relatively easy to account for. The problem is more severe for electrons, because of their smaller mass and consequently smaller gyroradii and inertial lengths, compared to protons, and is known as the electron injection problem.

Here we study electron injection at young SNR shock waves using particle-in-cell (PIC) numerical simulations that provide a fully self-consistent treatment of the electron scales. Such shocks are characterized by high sonic, $M_{\mathrm{s}}$, and Alfvénic, $M_{\mathrm{A}}$, Mach numbers. Present observational data do not give clear constraints on the large-scale magnetic-field configuration in portions of SNR shocks from which strong nonthermal emission is detected. Radio polarimetry are notoriously difficult to interpret (e.g., Stroman & Pohl, 2009). Different approaches of data modeling for the same source can suggest the presence of quasi-perpendicular fields (Petruk et al., 2009; Schneiter et al., 2010; West et al., 2016) or the opposite, quasi-parallel configurations (Rothenflug et al., 2004; Bocchino et al., 2011; Schneiter et al., 2015). As in our recent studies (Matsumoto et al., 2012, 2013, 2015; Wieland et al., 2016; Bohdan et al., 2017), in this work we examine perpendicular shocks as the most simple form of a quasi-perpendicular magnetic-field configuration. The physics of such shocks is governed by reflection of ions at the shock caused by shock potential (Fig. 1), the interaction of which with the incoming plasma excites a variety of instabilities upstream of the shock. The most important instabilities in the regime of high Mach numbers are the electrostatic two-stream Buneman instability at the leading edge of the foot, resulting from the interaction between cold incoming electrons and reflected ions (Buneman, 1958), and the Weibel instability in the shock foot driven by the interaction of the incoming and reflected ions (Kato & Takabe, 2010; Niemiec et al., 2012; Matsumoto et al., 2015; Wieland et al., 2016).

The Buneman instability can mediate the generation of supra-thermal electrons via shock surfing acceleration (SSA). In a 1D picture the Buneman instability produces strong, coherent electrostatic waves that capture electrons and let them be accelerated by the convective electric field (Hoshino & Shimada, 2002), thus providing for efficient electron injection. A number of 2D simulations of perpendicular shocks (Amano & Hoshino, 2009a; Matsumoto et al., 2012, 2013; Wieland et al., 2016) demonstrated that the length of the potential wells is limited to about the ion inertial length. Electrons can thus escape from the trapping region and re-enter it from the downstream or the upstream side to experience multiple surfing-acceleration events (Amano & Hoshino, 2009a; Matsumoto et al., 2012).

The Weibel instability generates strong magnetic fields with filamentary structure. It was also recently shown with 2D simulations that spontaneous turbulent magnetic reconnection in the Weibel instability region canlead to electron acceleration (Matsumoto et al., 2015). Thin current sheets (magnetic filaments) become unstable and break up into chains of magnetic islands and X-points. Particles can be accelerated while interacting with these structures.

The spectrum of waves generated at the shock is usually at least two-dimensional. Which of the unstable modes appear in a 2D simulation strongly depends on the configuration of the mean magnetic field though, as modes may be artificially suppressed if their wave vector is not contained in the simulation plane. In Bohdan et al. (2017) we showed that the Weibel instability is best reproduced with the in-plane setup, whereas the Buneman modes are considerably stronger and more coherent with a strictly out-of-plane orientation. Suprathermal tails in the electron spectra are found for all simulated shocks, and the initial acceleration of electrons always occurs through the SSA process in the Buneman wave region. However, the subsequent stages of injection strongly depend on the field configuration. For out-of-plane field adiabatic heating dominates the spectral evolution. For configurations with an in-plane magnetic-field component particles are non-adiabatically accelerated in interactions with turbulent magnetic structures in the shock, resembling a second-order Fermi process, and magnetic reconnection does also occur. The fraction of nonthermal electrons is an order of magnitude larger for the out-of-plane configuration than for other field orientations, mainly on account of a higher SSA efficiency.

The first 3D PIC simulation of a high-$M_{\mathrm{A}}$ shock was recently presented by Matsumoto et al. (2017) for an oblique subluminal configuration, $c/\tan{\Theta_{\rm Bn}}>v_{\mathrm{sh}}$, where $\Theta_{\rm Bn}$ is the angle of the large-scale magnetic field with respect to the shock normal, $v_{\mathrm{sh}}$ is the shock velocity, and $c$ is the speed of light. Buneman waves and Weibel magnetic turbulence were found to coexist in the shock structure. Energetic electrons that initially experienced SSA underwent pitch-angle diffusion by interacting with magnetic turbulence in the shock foot and ramp. This provides confinement in the shock transition region during which particles gain energy by shock drift acceleration (SDA). The computational cost of 3D experiments is still too high to sample the range of plasma conditions that one may find in SNR shocks. Nevertheless, the 3D results indicate which parts of 3D shock physics can be reliably probed with 2D simulations.

In this work we report on new large-scale 2D fully kinetic PIC simulations of nonrelativistic strictly perpendicular shocks in the regime of high Mach numbers, $M_{\mathrm{A}}\gtrsim 20$ and $M_{\mathrm{s}}\gtrsim 30$, as appropriate for forward shocks of young SNRs. The simulations are conducted in 2D3V configuration, i.e., we follow two spatial coordinates and all three components of the velocity and the electromagnetic fields. Numerical experiments are performed for both in-plane and out-of-plane configurations of the large-scale magnetic field. These simulations complement our previous investigations of 2D perpendicular shocks (e.g., Matsumoto et al., 2012, 2013, 2015; Wieland et al., 2016; Bohdan et al., 2017). The aim of this work is to analyze in detail the initial energization via SSA in the Buneman-instability region. The successive acceleration in the shock foot and ramp on account of, e.g., inelastic scattering off the Weibel-instability turbulence is the subject of a separate publication.

Conditions for efficient electron energization via SSA were first investigated by Matsumoto et al. (2012), supported with PIC simulations with out-of-plane magnetic-field configuration. The process occurs in low-temperature (low beta) plasmas, in which the Buneman instability can effectively grow. For efficient acceleration the electrostatic waves should also be strong enough to trap electrons and hold them during acceleration, which defines a minimum Alfvénic Mach number for a shock to be capable of producing relativistic electrons via SSA,

[TABLE]

where $\alpha$ is the flux ratio of reflected to incoming ions and $m_{\mathrm{i}}$ and $m_{\mathrm{e}}$ are the ion and the electron mass, respectively. In the presence of an in-plane magnetic field the motion of the reflected ions is not fully contained in the simulation grid and thus the corresponding component of the Buneman waves cannot be captured (Bohdan et al., 2017). To account for this effect we proposed a modified trapping condition:

[TABLE]

where $\varphi$ is the orientation angle of the large-scale perpendicular magnetic field with respect to the simulation plane, with $\varphi=0^{\rm o}$ representing the in-plane configuration (see Fig. 2). The earlier 2D simulations of Bohdan et al. (2017) all satisfied the trapping condition of Equation 1 and were performed for a single value of the reduced mass ratio, $m_{\mathrm{i}}/m_{\mathrm{e}}=100$, and a small ($\beta_{\rm e}\ll 1$) or moderate ($\beta_{\rm e}=0.5$) plasma beta. Our present work augments this analysis with investigations of the trapping conditions of Equations 1 and 2 and SSA efficiency for different mass ratios in the range $m_{\mathrm{i}}/m_{\mathrm{e}}=50-400$. Matsumoto et al. (2017) demonstrated that the SSA process is well reproduced with 2D out-of-plane simulations, but processes in the shock ramp and overshoot are suppressed. On the other hand, the stochastic Fermi-like acceleration in the Weibel-instability-generated turbulence works similar as in 2D in-plane experiments. If the modified trapping condition would define the parameter range, for which we have the same efficiency of electron pre-acceleration for an in-plane configuration as that observed in 2D simulations with out-of-plane magnetic field, it would be possible to reproduce realistic 3D physics with far cheaper 2D experiments with $\varphi=0^{\rm o}$. This is the main hypothesis under discussion here.

The paper is organized as follows. We present a description of the simulation setup in Section 2. The results are presented in Section 3. Section 4 contains the summary and discussion.

2 Simulation Setup

The simulation setup adopted in this work is the same as that used in Bohdan et al. (2017) and illustrated in Figure 2. As a result of the collision of two counter-streaming electron-ion plasma beams, two shocks are formed that propagate in opposite directions and are separated by a contact discontinuity (CD). The plasma flow is set along the $x$-direction in the $xy$ plane. Plasma particles are continuously injected at both sides of the simulation box with velocities $\boldsymbol{v}_{\rm L}=v_{\rm L}\hat{\boldsymbol{x}}$ and $\boldsymbol{v}_{\rm R}=v_{\rm R}\hat{\boldsymbol{x}}$, where the indices L and R refer, respectively, to the left and right sides of the simulation box. As the two shocks move away from the CD in the left and the right plasma, we refer to them as to the left and the right shocks, respectively. The two plasma streams carry a homogeneous magnetic field, $\boldsymbol{B}_{\rm 0}$, that is perpendicular to the shock normal and lies in the $yz$ plane. The magnetic field thus forms an angle $\varphi$ with the $y$-axis. Initialized with the flow is a motional electric field $\boldsymbol{E}_{\rm 0}=-\boldsymbol{v}\times\boldsymbol{B}_{\rm 0}$, with $\boldsymbol{v}=\boldsymbol{v}_{\rm L}$ or $\boldsymbol{v}=\boldsymbol{v}_{\rm R}$, respectively, for the left and the right beam. We assume that the beams move with equal absolute velocities, $v_{\rm L}=v_{\rm R}=0.2c$, and that the magnetic field strength in both plasmas is equal, $\boldsymbol{B}_{\rm 0L}=\boldsymbol{B}_{\rm 0R}$. The motional electric field thus has equal strength and opposing signs in the two slabs. We use the method of Wieland et al. (2016) to suppress the artificial electromagnetic transient that results from the initial strong electric-field gradient between the two plasma slabs.

We collide plasma beams of equal density but different temperatures, thus studying two different shocks in one simulation. The temperature ratio between the two beams is $1000$, so that the sonic Mach numbers, $M_{\rm s}$, of the two shocks differ by a factor of $\sqrt{1000}\simeq 30$. In terms of the electron plasma beta (the ratio of the electron plasma pressure to the magnetic pressure) the left beam has $\beta_{\rm e,L}=5\cdot 10^{-4}$ and the right beam $\beta_{\rm e,R}=0.5$. This choice of plasma beta facilitates a direct comparison with our earlier work (Bohdan et al., 2017) and also with results of previous 2D simulations of perpendicular shocks (Matsumoto et al., 2012, 2013) and a recent 3D simulation of a quasi-perpendicular shock (Matsumoto et al., 2017), in which $\beta_{\rm e,R}=0.5$ is assumed. Note that our system is approximately in ram-pressure balance, and consequently the simulation frame is also the downstream rest frame of the two shocks.

The parameters of the simulation runs described in this paper are listed in Table 1. We have performed seven large-scale numerical experiments (runs A–G), that feature in total fourteen simulated shocks. Here we refer to each of these shock cases as to a separate simulation run, and tag the shocks in the left plasma ($\beta_{\rm e,L}=5\cdot 10^{-4}$) with *1, and the right shocks with *2 ($\beta_{\rm e,R}=0.5$). Simulation runs A–F assume the in-plane magnetic field configuration, $\varphi=0^{o}$, and run G uses the out-of-plane magnetic field orientation, $\varphi=90^{o}$. We do not consider simulations with $\varphi=45^{o}$, because the shock structure and the acceleration mechanisms observed in this case are almost identical to those in runs with the in-plane field configuration (Bohdan et al., 2017). The runs with the in-plane magnetic field cover a wide range of ion-to-electron mass ratios and Alfvénic Mach numbers, as illustrated in Figure 3, which permits an investigation of the influence of these parameters on the electron acceleration efficiency and to scale our results to the realistic ion-to-electron mass ratio. Note, that some aspects of the shock physics in runs B and G have been already discussed in our previous paper (cf. runs A and C in Bohdan et al., 2017).

The derived shock properties are also listed in Table 1. The Alfvén velocity is defined as $v_{\rm A}=B_{\rm 0}/\sqrt{\mu_{\rm 0}(N_{e}m_{\mathrm{e}}+N_{i}m_{\mathrm{i}})}$, where $\mu_{\rm 0}$ is the vacuum permeability, $N_{i}$ and $N_{e}$ are the ion and the electron number densities, and $B_{0}$ is the far-upstream magnetic-field strength. The sound speed reads $c_{\rm s}=(\Gamma k_{B}T_{\rm i}/m_{\mathrm{i}})^{1/2}$, where $k_{B}$ is the Boltzmann constant, $\Gamma$ is a nonrelativistic adiabatic index, and $T_{\rm i}$ is the ion temperature. The Alfvénic, $M_{\mathrm{A}}=v_{\mathrm{sh}}/v_{\rm A}$, and sonic, $M_{\mathrm{s}}=v_{\mathrm{sh}}/c_{\rm s}$, Mach numbers of the shocks in Table 1 are given in the conventional upstream reference frame. As the in-plane and the out-of-plane magnetic field lead to a different number of degrees of freedom, the adiabatic indices are different with $\Gamma=5/3$ and $\Gamma=2$, respectively for $\varphi=0^{o}$ and $\varphi=90^{o}$. Thus the resulting expected shock speeds take values $v_{\mathrm{sh}}=0.263c$ for runs A–F and $v_{\mathrm{sh}}=0.294c$ for runs G. In the simulation frame the speeds are smaller by the shock compression ratio.

To investigate the role of SSA in electron pre-acceleration, we adjust the magnetic-field strength, $B_{0}$, to establish Alfvénic Mach numbers that test the trapping conditions defined by Equations 1 and 2. A comparison of the Alfvénic Mach numbers and the mass ratio of all runs with trapping limits is offered in Figure 3. Nevertheless, we always consider weakly magnetized plasmas with the ratio of the electron plasma frequency, $\omega_{\rm pe}=\sqrt{e^{2}N_{e}/\epsilon_{0}m_{\mathrm{e}}}$, to the electron gyrofrequency, $\Omega_{\rm e}=eB_{0}/m_{\mathrm{e}}$, in the range $\omega_{\rm pe}/\Omega_{\rm e}=8.5-17.3$. Here, $e$ is the electron charge, and $\epsilon_{0}$ is the vacuum permittivity. To keep the plasma beta constant we adjust the plasma temperatures and hence the sound speeds and resulting sonic Mach numbers (see Table 1).

In this work we want to verify several hypotheses. The first is the scaling of the SSA efficiency with the ion-to-electron mass ratio for shocks that fulfill the trapping condition of Equation 1, here applied to the in-plane magnetic field configurations. Runs A, B, E, and F define the set of simulations conducted for $m_{\mathrm{i}}/m_{\mathrm{e}}=50,100,200$, and $400$, respectively.

The second objective is the modified trapping condition of Equation 2. We test this condition by conducting simulation runs C, which satisfy Equation 2 for $\alpha\leq 0.5$. The question to be addressed is whether 2D simulations with in-plane magnetic field configuration can reproduce the SSA efficiency observed in 2D runs with the same $m_{\mathrm{i}}/m_{\mathrm{e}}$ and the out-of-plane fields, here marked as runs G.

The third set of simulations consists of runs D and E, performed for the same mass ratio $m_{\mathrm{i}}/m_{\mathrm{e}}=200$. The Alfvénic Mach number in run D clearly violates Equation 1, and so we expect a very low intensity of Buneman waves. Nevertheless, particle acceleration can still occur in the shock foot and ramp, whose structure is defined by the magnetic filaments, and we are interested in the nonthermal electron population that forms in the absence of SSA. Note, that cross-comparison of runs B and D, and C and E can yield the mass-ratio dependence for shocks having the same Alfvénic Mach numbers.

The electron skin depth in the upstream plasma is common for all runs and equals $\lambda_{\mathrm{se}}=20\Delta$, where $\Delta$ is the size of grid cells. The ion skin depth, $\lambda_{\mathrm{si}}=\sqrt{m_{\mathrm{i}}/m_{\mathrm{e}}}\lambda_{\mathrm{se}}$, is used here as the unit of length. The time scale and all temporal dependencies are given in terms of the upstream ion Larmor frequency, $\Omega_{\rm i}$, where $\Omega_{\rm i}=eB_{0}/m_{\mathrm{i}}$. The simulation time is typically $t=(6-8)\Omega_{i}^{-1}$, which is enough to cover at least a few shock self-reformation cycles (see Bohdan et al., 2017). The time-step we use is $\delta t=1/40\,\omega_{\rm pe}^{-1}$.

The two plasma beams injected at sides of the simulation box are composed of an equal number of ions and electrons, $N_{\rm ppc}=20$. Electron and ion plasma pairs are initialized at the same locations to ensure the initial charge-neutrality of the system. There is no escape of particles from the computational box, and we use injection layers receding from the CD as in Bohdan et al. (2017), which helps alleviating numerical grid-Cerenkov effects and saves computational resources. The simulation box expands in $x$-direction during the run. The final size of a simulation box can reach $L_{x}\approx 280\lambda_{\rm si}$. The transverse size of the simulation box, $L_{y}=(8.2-24)\lambda_{\rm si}$, is large enough to cover several of the magnetic filaments, that are typically separated by $\sim\lambda_{\mathrm{si}}$, and at the same time limits the computational expense that grows quadratically with $m_{\mathrm{i}}/m_{\mathrm{e}}$. The largest simulation box of size $L_{x}\times L_{y}=(3264\times 96000)\Delta$ is used in run F with $m_{\mathrm{i}}/m_{\mathrm{e}}=400$. Open boundary conditions are imposed in the $x$-direction and periodic boundaries are applied in the $y$-direction.

The numerical code we use is a 2D3V-adapted and modified version of the relativistic electromagnetic PIC code TRISTAN (Buneman, 1993) with MPI-based parallelization (Niemiec et al., 2008; Wieland et al., 2016) and the option to trace individual particles.

3 Results

In Section 3.1 we describe the structure of the Buneman wave modes in all simulations and also summarize the findings of Bohdan et al. (2017). Then we discuss the electron acceleration efficiency through SSA in Section 3.2.

3.1 The Buneman Instability

Figure 4 presents the maps of the electrostatic field amplitude in the foot of the right shocks (runs A2-G2, see Table 1), that propagate in moderate-temperature plasmas with $\beta_{\rm e}=0.5$. Only portions of the simulation boxes are shown to facilitate one-to-one comparison between the runs. Run G2* is run G2 at a different phase of shock reformation. The electrostatic fields are calculated as $|E_{\mathrm{ES}}|=|-\nabla\phi|$, where $\phi$ is the electric potential, that is derived directly from the charge distribution. The maps are plotted for simulation times, at which the cyclic shock self-reformation allows the strongest Buneman modes. Note, that the maps for runs B2 and G2 can be compared with Figures 6a3 and 6c3, respectively, in Bohdan et al. (2017), in which results for runs B1 and G1 are presented (marked as runs A1 and C1, respectively).

The properties of the Buneman instability discussed in Bohdan et al. (2017) can be readily observed in Figure 4. The wave vectors are approximately parallel to the shock normal for the in-plane configurations (runs A2-F2) and oblique for out-of-plane magnetic field (run G2). This reflects the motion of shock-reflected ions: for $\varphi=0^{\rm o}$ the ions are confined to the $xz$-plane whereas for $\varphi=90^{\rm o}$ they stream in the simulation plane. The Buneman wave region shows a patchy structure for the in-plane field configurations, that can be linked to clumps in the overshoot produced by merging magnetic filaments. In total, the Buneman waves occupy a much smaller region than for the out-of-plane configuration, for which the waves are coherent and more intense.

The phase velocity of the Buneman modes matches the relative speed between shock-reflected ions and incoming electrons of the upstream plasma. Since for $\varphi=0^{o}$ part of the ion motion is outside of the simulation grid, the wavelengths of the Buneman waves are smaller ($\lambda\approx 1.9\lambda_{\mathrm{se}}$) than for out-of-plane field, for which $\lambda\approx 3.3\lambda_{\mathrm{se}}$. Note, that Figure 4 shows $|E|$ and hence the wavelength is twice the separation of wave fronts, here provided in units of the ion skin depth. The surface area of the Buneman wave region for shocks in moderate-temperature plasma is 20%-30% larger than at the corresponding low-$\beta$ shocks, but the intensity of the waves is 20%-50% smaller (compare Fig. 6 in Bohdan et al., 2017).

For the high-$\beta$ systems presented in Figure 4, Table 2 lists peak amplitude of Buneman waves and the fraction of pre-accelerated electrons. The runs A2, B2, E2 and F2 satisfy the trapping condition of Equation 1 (see Fig. 3), and both the peak and average strength of the electrostatic field are similar. Small differences between them arise from shock reformation. We conclude that irrespective of the mass ratio, the physical conditions at shocks with $M_{\mathrm{A}}$ satisfying Equation 1 are similar. However, the electrostatic force is weaker in average than the Lorentz force on a $\gamma\gtrsim 2$ electron ($|E_{\mathrm{ES}}|/(cB_{0})<1$).

Considerably larger electrostatic field amplitudes, reaching $|E_{\mathrm{ES}}|/(cB_{0})\sim\!2.3$, can be observed for run C2. Here, the Alfvén Mach number of the shock, $M_{\mathrm{A}}=46$, is much larger than the minimum $M_{\mathrm{A}}$ defined by Equation 1 and also satisfies the modified trapping condition of Equation 2, that for the measured $\alpha\simeq 0.32$ gives the minimum $M_{\mathrm{A}}\simeq 40.2$. The field intensity in run C2 is about a factor of 2 larger than in both run B2 with the same mass ratio, $m_{\mathrm{i}}/m_{\mathrm{e}}=100$, and run E2 with mass ratio $m_{\mathrm{i}}/m_{\mathrm{e}}=200$ but similar Alfvén Mach number, $M_{\mathrm{A}}\simeq 45$. This shows that the strength of the electrostatic modes is driven by the value of the Alfvénic Mach number in relation to the trapping condition (Eq. 1). The absolute value of $M_{\mathrm{A}}$ is not important, as in run D2 we see Buneman waves with amplitudes a factor of 3 lower than those in run B2 with the same Alfvénic Mach number. Essentially all observed wave intensities are slightly weaker than the saturation level estimated by Ishihara et al. (1980).

The modified trapping condition (Eq. 2) was expected to compensate for the effect of the field configuration. Shocks with sufficiently large $M_{\mathrm{A}}$ should then reproduce similar Buneman wave intensities in 2D in-plane magnetic field configurations than in simulations with the out-of-plane fields. However, the electrostatic field in run C2 is weaker by 20% than that in run G2 with $\varphi=90^{o}$. At a different phase of shock reformation run G2, now called G2*, has the same electric-field amplitude as C2, but four times the number of pre-accelerated electrons. This discrepancy might arise from Equation 2 only compensating for the neglect of the $z$-motion of ions. In out-of-plane simulation we observe that the relative speed between electrons and reflected ions can reach $\sim 0.6c$, because of acceleration in upstream electric field, which is a factor of $\sim 1.5$ larger then the value assumed in the derivation of the trapping condition (see Matsumoto et al., 2012), while in in-plane case the acceleration is in $z$-direction. It may be that we need to also account for this effect by adding a factor of 1.5 to the modified trapping condition,

[TABLE]

This equation gives $M_{\mathrm{A}}\simeq 60.3$ for the minimum Alfvén Mach number, with which the amplitudes of the Buneman waves observed at shocks with $M_{\mathrm{A}}=35.5$ in 2D simulations with out-of-plane magnetic fields could be reproduced in runs applying $\varphi=0^{o}$ field configuration. This value is much larger than any of the Mach numbers studied here for $m_{\mathrm{i}}/m_{\mathrm{e}}=100$, and thus requires attention in the future.

3.2 Electron Acceleration in the Buneman Zone

Table 2 lists the fraction of electrons that have been pre-accelerated in the Buneman wave zone to $(\gamma-1)>0.1$, $N_{\rm e,BI}/N_{\rm e,tot}$. This fraction is much larger in run G2 than it is in runs A2-F2. Bohdan et al. (2017) argued that at least part of this difference is due to differences in the amplitude of the electrostatic waves and their coverage area.

Figures 5 and 6 show kinetic-energy spectra of electrons occupying the Buneman wave regions highlighted in Figure 4. Figure 5 shows energy spectra for runs A2, B2, E2 and F2, for which the Alfvénic Mach numbers exceed by the similar margin the trapping condition (Eq. 1). The spectra are statistically indistinguishable, and the fraction of pre-accelerated electrons is $\sim 0.45\%$ for all runs. This is again in line with the ion-to-electron mass ratio dependence of the trapping condition.

Electron spectra for runs B2, C2, and D2, that probe different physical conditions at shocks with in-plane magnetic-field configuration, are compared to the spectrum for the out-of-plane case G2 in Figure 6. The fractions of pre-accelerated electrons differ between the in-plane runs (see Table 2), reflecting the different intensities of the Buneman waves. In run C2, the spectrum extends to higher energies and contains more energetic electrons than that for run B2, which arises from the difference in Mach number. The Alfvénic Mach number of the shock in run D2 instead violates the trapping condition, and only $\sim 0.34\%$ of electrons are pre-accelerated.

Although run C2 satisfies the modified trapping condition of Equation 2, the acceleration efficiency, $N_{\rm e,BI}/N_{\rm e,tot}\simeq 0.6$, is much less than for run G2. In Figure 6 the spectrum for run C2 is also compared with the spectrum calculated for run G2 at a different phase of the shock-reformation (denoted as run G2*), at which the strength of the Buneman waves matches that for run C2. Still, the fraction of pre-accelerated electrons in run G2* is four times that in run C2, but the maximum energies of the electrons are comparable, $\max(\gamma)\approx 3-4$. It is clear that the Buneman wave strength is not the only parameter that determines the efficiency of SSA in the shock foot.

SSA consists of two individual processes: (1) interaction with electrostatic waves and (2) magnetic gyration. In the appendix we present a detailed analytical treatment of the equation of motion of electrons in the wave field, demonstrating that the electrostatic field of the waves does the physical work. Here we summarize the conclusions.

Figure 7 illustrates the first stage of the SSA process for the in-plane (left panels a1-e1) and the out-of-plane case (right panels a2-e2). For specific electrons extracted from runs E2 and G2, we see the time evolution of the energy (Fig. 7b) and the momentum (Fig. 7c), as well as the electric field at the location of the particle in the simulation frame and in the instantaneous particle rest frame (Fig. 7d and e, respectively). The latter is particularly interesting, because in the electron rest frame the electric field is the sole provider of acceleration. We refer to the selected electron in the in-plane case (left panels) as the first electron and the other one as the second electron. Initially both electrons move with the plasma bulk. To be trapped by electrostatic waves, the electrons must travel with the waves against the upstream plasma flow, and hence be picked-up from the thermal pool. Before doing so, the electrons move in the negative x-direction undisturbed through several electrostatic wavefronts. Significant energy gain commences at time $t\omega_{\mathrm{pe}}=6270$ for the first electron and at $t\omega_{\mathrm{pe}}=3745$ for the second electron. The particles then remain trapped by the waves and undergo the first stage of acceleration at time intervals $t\omega_{\mathrm{pe}}=(6275-6300)$ for the first electron and $t\omega_{\mathrm{pe}}=(3745-3753)$ for the second electron. During this stage both electrons move in the direction of shock propagation, and their $p_{y}$ momentum remains small. The end of the first-stage acceleration is marked by the black vertical line in Figure 7, beyond which the electrons resume gyrating.

The acceleration of the first electron occurs in the same way as in 1D geometry (Hoshino & Shimada, 2002): the electron is pushed toward the upstream region by the electrostatic field of a Buneman wave, which for some time compensates the $x$-component of the Larmor acceleration and thus keeps the electron roughly in phase with the wave. Consequently the average values of $E_{x}$ and $E_{y}$ electric field components are close to zero in the particle reference frame (Fig. 7e1). The continuous gradient in $p_{z}$ at $t\omega_{\mathrm{pe}}=(6275-6300)$ reflects the transverse Larmor acceleration, which can be described as the effect of the motional electric field in the frame of the electrostatic wave. It is important to note that for the in-plane magnetic field the wave fronts are infinitely extended in $z$-direction, and the energy gain terminates when the electron loses phase coherence with the Buneman wave. In reality the energization will terminate earlier. In the upstream flow frame all the energy gain comes from the field of the Buneman wave though.

The second electron displays a similar behaviour, but in the frame of the obliquely propagating waves. At $t\omega_{pe}\simeq 3745$ it starts moving in the x-direction, but the electrostatic field of the waves roughly compensates the Larmor acceleration in y-direction, as $E_{\mathrm{PRF},y}\approx 0$ (Fig. 7e2). Instead, the electron is accelerated in x-direction by the electrostatic field of the Buneman waves. We conclude that in all cases the energy gain arises from the electrostatic field of the waves, while the formal acceleration reflects the competition of Larmor acceleration and that imposed by the waves. The in-plane configuration captures only part of the Buneman waves, as only wave vectors in the simulation plane are allowed, and so there is a lower rate of energy gain compared to the out-of-plane case. In addition, the restriction of the wave phase velocity to the simulation plane changes the direction of sliding along a wave front from effectively the $x$-direction to the $z$-direction.

Let us estimate the energy gain arising from trapping at an electrostatic wavefront. Equations A5 and A8 give the rate of energy gain for the out-of-plane and the in-plane configuration, respectively. The phase speed that the electrons need to match is $v_{ph,0}=0.1c$ and $v_{ph,90}=0.4c$ for in-plane and out-of-plane configuration, respectively, and so $\Delta v$ and hence the energization rate is twice larger in the out-of-plane case than it is for in-plane magnetic field. The total energy gain is the product of the rate of gain and the time of interaction. The time of interaction is limited by three factors: the intermittency of waves, escape by acceleration perpendicular to the wave front, and escape to the side of the wave front.

In the in-plane case the wave front is infinitely extended in $z$ direction, and no escape to that side is possible. For an out-of-plane magnetic field and an average speed along the wave front of $\sim(0.1-0.2)\,c$, the electrons would escape trapping on $t_{\mathrm{esc}}\approx(25-50)\omega_{\mathrm{pe}}^{-1}$, as the wave fronts in Figure 7 have a lateral extent of about $5\,\lambda_{\mathrm{se}}$.

The escape time perpendicular to the wave front can be estimated as $t_{\mathrm{esc}}\approx\pi\omega_{\mathrm{pe}}^{-1}(v_{\Phi}+v_{0})/v_{\mathrm{e,WRF}}$, where $v_{\mathrm{e,WRF}}$ is the velocity of electrons in the wave frame. For the out-of-plane case this gives $t_{\mathrm{esc},90}\gtrsim 18\omega_{\mathrm{pe}}^{-1}$, as the average electron speed $v_{\mathrm{e,WRF}}\lesssim 0.1\,c$.

The trapping time coming from the wave time intermittency can be estimated directly from simulations. The ability to accelerate an electron up to a certain energy depends not only on the instantaneous local electrostatic field strength but also on the previous strength history and the ability to trap an electron during the whole acceleration period. In Figure 8 the time evolution of electrostatic field strength ($E_{\mathrm{ES}}$, black dash-dotted line) at a chosen location is presented. This field is able to trap electrons with energies shown with black solid line, which is calculated assuming the equality between electrostatic and Lorentz forces at the chosen location. At a time $t^{\prime}$ the electrostatic field is capable to trap an electron with energy $\epsilon_{1}$. However, taking into account the evolution of $E_{\mathrm{ES}}$ and the energy of electrons (red dashed line) this electron cannot be trapped during the whole acceleration period. Therefore electrons with a final energy $\epsilon_{2}$ and the energy history shown with the red solid line can be present in the simulation. According to these considerations the trapping time reads

[TABLE]

which is approximately the time that we analytically estimated based on the acceleration in the direction of the wave motion (see the Appendix). Thus one of the main limiting factors for electron acceleration is the intermittency of the Buneman waves.

Calculated average energy gains are

[TABLE]

which are similar to those for electrons in Figure 7 and average energies of accelerated electrons in Figures 5 and 6. The analytically expected energy increase can be written as

[TABLE]

where $e\,E_{\mathrm{ES}}=m_{e}\Delta v\omega_{\mathrm{pe}}(m_{e}/m_{i})^{(1/6)}$ (Ishihara et al., 1980; Amano & Hoshino, 2009b; Matsumoto et al., 2012) and $|F|$ is assumed to be about 1. Therefore the main difference in the acceleration rate comes from velocity difference, $\Delta v=(v_{\Phi}+v_{0})$, and the energy gain of electrons is still stronger in the out-of-plane case due to a larger phase speed of the Buneman waves.

We note that the modified trapping conditions (Eq. 2 or 3) refer to reaching a certain strength of the electrostatic field that is needed for trapping, while the energy gain of electrons is related to the velocity difference between reflected ions and upcoming electrons. This velocity difference imposes the main restriction for the in-plane simulations in their applicability to mimic realistic SSA efficiency. Using a higher Mach number can not significantly change the SSA efficiency in case of the same velocity difference defined by the magnetic field configuration. For the same mass ratio the number of pre-accelerated electrons is larger by about (30-40)% in the runs with a higher Alfvénic Mach number (see Table 2, runs B2-C2 and D2-E2), which is not the factor of 10 required to reach the SSA efficiency seen in out-of-plain runs. Therefore significant modifications of the parameters of the simulation (not just a change of the Alfvénic Mach number) are needed to reproduce the out-of-plain SSA efficiency by means of in-plane simulations.

The energy difference associated with climbing or sliding down the potential well of a Buneman wave can be estimated as $\Delta\gamma mc^{2}=eE_{BI}\,\lambda_{BI}/2\pi$. The wavelength of Buneman waves, $\lambda_{BI}=2\pi\Delta v/\omega_{\mathrm{pe}}$, then implies an energy change $\Delta\gamma\approx 0.05$ in the in-plane case and $\Delta\gamma\approx 0.17$ in the out-of-plane run. This is insufficient to redirect an incoming electron to stationarity in the wave frame. Fluctuations in the Buneman wave field are clearly needed to trap particles and keep them in phase with the waves.

We observe that for $\varphi=90^{\rm o}$ a larger number of electrons are picked up from the bulk plasma for further acceleration than is seen with the in-plane configuration, which can be explained by a twice stronger $e\,E_{\rm ES}$ force in the out-of-plane case.

4 Summary and discussion

We analyse electron injection processes at nonrelativistic perpendicular collisionless shocks with high Alfvénic Mach numbers with 2D3V numerical PIC simulations. Earlier studies indicated that SSA operates at the leading edge of the foot as first-stage electron pre-acceleration mechanism, provided the Alfvénic Mach number satisfies a condition of efficient driving of the electrostatic Buneman waves (the trapping condition, Matsumoto et al., 2012). In Matsumoto et al. (2015) and Bohdan et al. (2017) we showed that in 2D simulations that use a field component which lies in the simulation plane, the downstream nonthermal-electron fraction is much lower than with out-of-plane mean field. Noting that much of this difference results from an incomplete account of the Buneman instability in the in-plane geometry, and motivated by results of recent 3D studies which demonstrate that the injection physics past the SSA stage can adequately be studied with 2D in-plane simulations (Matsumoto et al., 2017), here we further investigate electron acceleration by SSA at perpendicular high-$M_{\mathrm{A}}$ shocks with in-plane magnetic field configurations. The aim is to infer the SSA efficiency, in particular the validity of the trapping condition in its original form and the variant proposed in Bohdan et al. (2017), and the relation to the SSA efficiency observed in simulations with the out-of-plane fields.

Our results can be summarized as follows:

• •

The energy gain in SSA always arises from the electrostatic field of a Buneman wave with which the electron travels for some time. The apparent acceleration, $\dot{\mathbf{v}}$, reflects the superposition of electrostatic acceleration and Larmor acceleration that might be described as effect of the motional electric field in the wave frame. This process is more efficient in the out-of-plane case because both the phase speed and the amplitude of the waves are higher than for $\varphi=0^{\rm o}$.

• •

As in high-$M_{\mathrm{A}}$ shock simulations with out-of-plane magnetic fields, for in-plane magnetic field the strength of the electrostatic wave modes in the shock foot is determined by the Alfvénic Mach number in relation to the trapping condition. The more $M_{\mathrm{A}}$ exceeds the trapping condition, the stronger the intensity of the Buneman waves. Shocks with Alfvénic Mach numbers satisfying the trapping condition by the similar margin show comparable wave strengths in simulations for different ion-to-electron mass ratios.

• •

Shocks in simulations with in-plane magnetic field demonstrate electrostatic wave intensities lower than those observed in the out-of-plane case, even if the modified trapping condition is satisfied.

• •

The trapping time is mostly defined by intermittency of, and limited phase-coherence of electrons with, the Buneman waves. This limits the duration of the velocity match between electrons and the waves.

• •

The number of electrons pre-accelerated via SSA in the shock foot strongly correlates with the strength of the electrostatic waves. Shocks with the same physical conditions defined through the trapping condition show similar SSA efficiency. The latter is proportional to $M_{\mathrm{A}}$ for a given mass ratio. However, SSA always produces larger fractions of pre-accelerated electrons in simulations with the out-of-plane configurations, even if the intensities of the Buneman waves are similar as in the in-plane case. One reason for that is the larger number of electrons being picked up from the bulk plasma for SSA compared to the in-plane configuration.

We conclude that with an in-plane magnetic-field configuration we can not achieve the same level of SSA efficiency as in simulations with out-of-plane magnetic field or 3D simulations (Matsumoto et al., 2017), unless the parameters and settings of the simulation setup are significantly modified.

This paper is conceived as the first of a series investigating different aspects of electron acceleration processes at non-relativistic perpendicular shocks using PIC simulations. Interaction with Weibel filaments and magnetic reconnection in the shock transition, plasma heating, and the generation of turbulent magnetic field will be covered in forthcoming publications.

We thank the anonymous referee for their comments. The work of J.N. has been supported by Narodowe Centrum Nauki through research project DEC-2013/10/E/ST9/00662. This work was supported by JSPS-PAN Bilateral Joint Research Project Grant Number 180500000671. The numerical experiment was possible through a 10 Mcore-hour allocation on the 2.399 PFlop Prometheus system at ACC Cyfronet AGH. Part of the numerical work was conducted on resources provided by the North-German Supercomputing Alliance (HLRN) under projects bbp00003 and bbp00014.

Appendix A Analytical model of electron SSA

A.1 Out-of-plane configuration, $\varphi=90^{\rm o}$

In the simulation frame, the large-scale magnetic field of the right plasma slab, $\mathbf{B}=B_{0}\,\hat{\boldsymbol{z}}$, induces a motional electric field, $\mathbf{E}=-v_{0}\,B_{0}\,\hat{\boldsymbol{y}}$, where $v_{0}$ is the speed of the upstream plasma flowing in $-x$ direction. The entire Larmor orbit of all particles with low temperature is leveled in the simulation plane, as are the acceleration imposed by the waves.

Suppose an electrostatic wave propagates at an angle $\Theta$ to the x-axis. The electric field carried by the wave is

[TABLE]

where the wave factor is

[TABLE]

Here we allow for an arbitrary phase, $\Phi$. The phase speed of the wave, $v_{\Phi}$, is measured in the simulation frame. The wave number is related to the velocity of reflected ions through the resonance condition of the Buneman modes, $\omega_{\mathrm{pe}}=k(v_{\Phi}+v_{0})$.

Now consider an electron with velocity components $v_{x}$ and $v_{y}$. Using non-relativistic kinematics we find the acceleration of the electron as

[TABLE]

Let us rotate the coordinate system by an angle $\Theta$, so that $x^{\prime}$ is oriented in the direction of motion of the waves and $y^{\prime}$ is perpendicular to it. The corresponding accelerations then read

[TABLE]

The wave factor, $F$, is explicitly time-dependent and may induce rapidly oscillating acceleration. The other terms only describe Larmor gyration in the flow frame and hence no real energy gain. The wave factor must be approximately constant, if continuous energy gain is to be achieved for about 10 plasma times, $\omega_{\mathrm{pe}}^{-1}$, as observed. This requires that on average $v_{x^{\prime}}-v_{\Phi}\lesssim 0.2\,c$ or roughly acceleration from $v_{x^{\prime}}=0.2\,c$ to $v_{x^{\prime}}=0.6\,c$, after which the electron is out of phase with the wave and commences Larmor motion.

The Larmor motion of the reflected ions mandates a wave direction for which $\sin\Theta$ is negative. Likewise, the wave factor, $F$, must be negative to effect energy gain. Equation A4 then indicates that acceleration in $y^{\prime}$ direction follows that in $x^{\prime}$ direction, and for a fair range of initial conditions $\dot{v}_{y^{\prime}}$ is slightly less than $\dot{v}_{x^{\prime}}$ and increases with the same rate, at least for up to $1\,\Omega_{\mathrm{e}}^{-1}\simeq 12\,\omega_{\mathrm{pe}}^{-1}$. Correspondingly, the momentum component $p_{x}$ increases approximately linearly, and the increase in speed is approximately $E_{\mathrm{ES}}/(2B_{0})$, whereas $p_{y}$ remains approximately constant.

The effective acceleration toward the upstream region arises from the superposition of acceleration in the electrostatic field of the Buneman waves and the Larmor acceleration, that are oppositely directed in $y$ direction, but both have positive components in $x$ direction. Acknowledging that both $F$ and $\sin\Theta$ must be negative, the rate of energy gain in the upstream flow frame is

[TABLE]

and hence completely independent of the motional electric field. In the simulation frame the velocity component $v_{0}\cos\Theta$ disappears from Equation A5 and a new component of energy-gain rate appears, $\Omega_{\mathrm{e}}v_{0}v_{y}$, which captures the apparent energy by Larmor motion in this frame.

A.2 In-plane configuration, $\varphi=0^{\rm o}$

The main impact of the in-plane configuration is that the part of the Larmor motion is perpendicular to the simulation plane, and so the orientation and properties of the Buneman waves are modified, as only wave vectors in the simulation plane can be captured. The wave factor changes to

[TABLE]

The acceleration then follows by appropriate rotation of that given in Equation A3,

[TABLE]

Obviously, there is linear acceleration in $-z$ direction, if the particle can be held at approximately constant phase ($F<0$; $v_{x}\approx v_{\Phi}$) in the wave. As $v_{\Phi}\gtrsim v_{0}$ it is the Larmor acceleration that is responsible for the particle’s sliding along the wavefront, and the electrostatic field of the waves provides slow energy gain at a rate

[TABLE]

which also only involves the electrostatic field of the Buneman waves. The energy gain will be less than that for out-of-plane configuration, because only part of the motion of the back-streaming ions can drive waves that hence have lower amplitude, $E_{\mathrm{ES}}$, and additionally the velocity term in Equation A8 is reduced.