Kinetic simulations of nonrelativistic perpendicular shocks of young supernova remnants. I. Electron shock-surfing acceleration
Artem Bohdan, Jacek Niemiec, Martin Pohl, Yosuke Matsumoto, Takanobu Amano, Masahiro Hoshino
DDraft version May 1, 2019
Typeset using L A TEX twocolumn style in AASTeX62
Kinetic simulations of nonrelativistic perpendicular shocks of young supernova remnants. I. Electron shock-surfingacceleration.
Artem Bohdan, Jacek Niemiec, Martin Pohl,
1, 3
Yosuke Matsumoto, Takanobu Amano, andMasahiro Hoshino DESY, 15738 Zeuthen, Germany Institute of Nuclear Physics Polish Academy of Sciences, PL-31342 Krakow, Poland Institute of Physics and Astronomy, University of Potsdam, 14476 Potsdam, Germany Department of Physics, Chiba University, 1-33 Yayoi-cho, Inage-ku, Chiba 263-8522, Japan Department of Earth and Planetary Science, the University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan
Submitted to ApJABSTRACTElectron injection at high Mach-number nonrelativistic perpendicular shocks is studied here forparameters that are applicable to young SNR shocks. Using high-resolution large-scale two-dimensionalfully kinetic particle-in-cell (PIC) simulations and tracing individual particles we in detail analyze theshock surfing acceleration (SSA) of electrons at the leading edge of the shock foot. The central questionis to what degree the process can be captured in 2D3V simulations. We find that the energy gain inSSA always arises from the electrostatic field of a Buneman wave. Electron energization is moreefficient in the out-of-plane orientation of the large-scale magnetic field because both the phase speedand the amplitude of the waves are higher than for the in-plane scenario. Also, a larger number ofelectrons is trapped by the waves compared to the in-plane configuration. We conclude that significantmodifications of the simulation parameters are needed to reach the same level of SSA efficiency as insimulations with out-of-plane magnetic field or 3D simulations.
Keywords: acceleration of particles, instabilities, ISM – supernova remnants, methods – numerical,plasmas, shock waves INTRODUCTIONThe current paradigm of cosmic-ray (CR) origin as-sumes that the most part of galactic CR populationis produced at nonrelativistic forward shocks of super-nova remnants (SNRs). The main acceleration mech-anism considered at shocks is diffusive shock acceler-ation (DSA), a first-order Fermi process (e.g., Axfordet al. 1977; Drury 1983; Blandford & Eichler 1987).Astronomical observations give strong support to thisparadigm. In particular, detection of broadband non-thermal emission from SNRs, extending in some objectsto TeV-range gamma rays, proves the presence of ultra-relativistic particles in these sources, though for mostSNRs it is still unclear which parent particle populations
Corresponding author: Artem [email protected] (protons or electrons) generate dominant high-energyemission (Aharonian 2013).Acceleration of particles through DSA comes frommultiple interactions with the shock front, while theybounce between the shock upstream and downstreamplasmas. Particle confinement to the shock vicinity isprovided by elastic scattering off magnetohydrodynamic(MHD) turbulence that renders diffusive particle mo-tions. The critical ingredient and the main unsolvedproblem in the DSA theory is the particle injection. CRsundergoing DSA have Larmor radii much larger than theinternal shock transition width, that is commensuratewith the gyroradius of the incoming protons (with shockspeed v sh ). CRs thus see the shock as a sharp disconti-nuity in the plasma flow. To be fed into the accelerationprocess particles need therefore to be extracted from thethermal pool and pre-accelerated. Since protons have alarger initial momentum and can be easily scattered ei-ther by MHD waves embedded in the ambient plasma a r X i v : . [ a s t r o - ph . H E ] A p r Bohdan et al. or by self-generated turbulence, their injection is rela-tively easy to account for. The problem is more severefor electrons, because of their smaller mass and conse-quently smaller gyroradii and inertial lengths, comparedto protons, and is known as the electron injection prob-lem.Here we study electron injection at young SNR shockwaves using particle-in-cell (PIC) numerical simulationsthat provide a fully self-consistent treatment of the elec-tron scales. Such shocks are characterized by high sonic, M s , and Alfv´enic, M A , Mach numbers. Present obser-vational data do not give clear constraints on the large-scale magnetic-field configuration in portions of SNRshocks from which strong nonthermal emission is de-tected. Radio polarimetry are notoriously difficult tointerpret (e.g., Stroman & Pohl 2009). Different ap-proaches of data modeling for the same source can sug-gest the presence of quasi-perpendicular fields (Petruket al. 2009; Schneiter et al. 2010; West et al. 2016) or theopposite, quasi-parallel configurations (Rothenflug et al.2004; Bocchino et al. 2011; Schneiter et al. 2015). As inour recent studies (Matsumoto et al. 2012, 2013, 2015;Wieland et al. 2016; Bohdan et al. 2017), in this work weexamine perpendicular shocks as the most simple form ofa quasi-perpendicular magnetic-field configuration. Thephysics of such shocks is governed by reflection of ions atthe shock caused by shock potential (Fig. 1), the interac-tion of which with the incoming plasma excites a varietyof instabilities upstream of the shock. The most impor-tant instabilities in the regime of high Mach numbersare the electrostatic two-stream Buneman instability atthe leading edge of the foot, resulting from the interac-tion between cold incoming electrons and reflected ions(Buneman 1958), and the Weibel instability in the shockfoot driven by the interaction of the incoming and re-flected ions (Kato & Takabe 2010; Niemiec et al. 2012;Matsumoto et al. 2015; Wieland et al. 2016).The Buneman instability can mediate the generationof supra-thermal electrons via shock surfing acceleration(SSA). In a 1D picture the Buneman instability pro-duces strong, coherent electrostatic waves that captureelectrons and let them be accelerated by the convec-tive electric field (Hoshino & Shimada 2002), thus pro-viding for efficient electron injection. A number of 2Dsimulations of perpendicular shocks (Amano & Hoshino2009a; Matsumoto et al. 2012, 2013; Wieland et al. 2016)demonstrated that the length of the potential wells islimited to about the ion inertial length. Electrons canthus escape from the trapping region and re-enter itfrom the downstream or the upstream side to experiencemultiple surfing-acceleration events (Amano & Hoshino2009a; Matsumoto et al. 2012). Figure 1.
Perpendicular shock structure. Top panel is theparticle number density profile. The shock transition consistsof a foot, a ramp, an overshoot and the downstream region. E x is the shock potential. v and v sh are the upstream andthe shock velocities. Bottom panel is the x-component of ionphase-space distribution. The Weibel instability generates strong magneticfields with filamentary structure. It was also recentlyshown with 2D simulations that spontaneous turbu-lent magnetic reconnection in the Weibel instabilityregion canlead to electron acceleration (Matsumoto etal. 2015). Thin current sheets (magnetic filaments) be-come unstable and break up into chains of magneticislands and X-points. Particles can be accelerated whileinteracting with these structures.The spectrum of waves generated at the shock is usu-ally at least two-dimensional. Which of the unstablemodes appear in a 2D simulation strongly depends onthe configuration of the mean magnetic field though, asmodes may be artificially suppressed if their wave vectoris not contained in the simulation plane. In Bohdan etal. (2017) we showed that the Weibel instability is bestreproduced with the in-plane setup, whereas the Bune-man modes are considerably stronger and more coher-ent with a strictly out-of-plane orientation. Suprather-mal tails in the electron spectra are found for all sim-ulated shocks, and the initial acceleration of electronsalways occurs through the SSA process in the Bunemanwave region. However, the subsequent stages of injec-tion strongly depend on the field configuration. For out-of-plane field adiabatic heating dominates the spectralevolution. For configurations with an in-plane magnetic-field component particles are non-adiabatically acceler-ated in interactions with turbulent magnetic structuresin the shock, resembling a second-order Fermi process,and magnetic reconnection does also occur. The frac- lectron SSA at SNR shocks M A shockwas recently presented by Matsumoto et al. (2017) foran oblique subluminal configuration, c/ tan Θ Bn > v sh ,where Θ Bn is the angle of the large-scale magnetic fieldwith respect to the shock normal, v sh is the shock ve-locity, and c is the speed of light. Buneman waves andWeibel magnetic turbulence were found to coexist in theshock structure. Energetic electrons that initially expe-rienced SSA underwent pitch-angle diffusion by inter-acting with magnetic turbulence in the shock foot andramp. This provides confinement in the shock transitionregion during which particles gain energy by shock driftacceleration (SDA). The computational cost of 3D ex-periments is still too high to sample the range of plasmaconditions that one may find in SNR shocks. Never-theless, the 3D results indicate which parts of 3D shockphysics can be reliably probed with 2D simulations.In this work we report on new large-scale 2D fullykinetic PIC simulations of nonrelativistic strictly per-pendicular shocks in the regime of high Mach numbers, M A (cid:38)
20 and M s (cid:38)
30, as appropriate for forwardshocks of young SNRs. The simulations are conductedin 2D3V configuration, i.e., we follow two spatial coor-dinates and all three components of the velocity and theelectromagnetic fields. Numerical experiments are per-formed for both in-plane and out-of-plane configurationsof the large-scale magnetic field. These simulations com-plement our previous investigations of 2D perpendicularshocks (e.g., Matsumoto et al. 2012, 2013, 2015; Wielandet al. 2016; Bohdan et al. 2017). The aim of this workis to analyze in detail the initial energization via SSA inthe Buneman-instability region. The successive acceler-ation in the shock foot and ramp on account of, e.g.,inelastic scattering off the Weibel-instability turbulenceis the subject of a separate publication.Conditions for efficient electron energization via SSAwere first investigated by Matsumoto et al. (2012),supported with PIC simulations with out-of-planemagnetic-field configuration. The process occurs in low-temperature (low beta) plasmas, in which the Bunemaninstability can effectively grow. For efficient accelerationthe electrostatic waves should also be strong enough totrap electrons and hold them during acceleration, whichdefines a minimum Alfv´enic Mach number for a shock tobe capable of producing relativistic electrons via SSA, M A ≥ (1 + α ) (cid:18) m i m e (cid:19) , (1) where α is the flux ratio of reflected to incoming ionsand m i and m e are the ion and the electron mass, re-spectively. In the presence of an in-plane magnetic fieldthe motion of the reflected ions is not fully containedin the simulation grid and thus the corresponding com-ponent of the Buneman waves cannot be captured (Bo-hdan et al. 2017). To account for this effect we proposeda modified trapping condition: M A ≥ (cid:115)
21 + sin ϕ (1 + α ) (cid:18) m i m e (cid:19) , (2)where ϕ is the orientation angle of the large-scale per-pendicular magnetic field with respect to the simulationplane, with ϕ = 0 o representing the in-plane configu-ration (see Fig. 2). The earlier 2D simulations of Bo-hdan et al. (2017) all satisfied the trapping condition ofEquation 1 and were performed for a single value of thereduced mass ratio, m i /m e = 100, and a small ( β e (cid:28) β e = 0 .
5) plasma beta. Our present workaugments this analysis with investigations of the trap-ping conditions of Equations 1 and 2 and SSA efficiencyfor different mass ratios in the range m i /m e = 50 − ϕ = 0 o .This is the main hypothesis under discussion here.The paper is organized as follows. We present a de-scription of the simulation setup in Section 2. The re-sults are presented in Section 3. Section 4 contains thesummary and discussion. SIMULATION SETUPThe simulation setup adopted in this work is the sameas that used in Bohdan et al. (2017) and illustrated inFigure 2. As a result of the collision of two counter-streaming electron-ion plasma beams, two shocks areformed that propagate in opposite directions and areseparated by a contact discontinuity (CD). The plasmaflow is set along the x -direction in the xy plane. Plasmaparticles are continuously injected at both sides of thesimulation box with velocities v L = v L ˆ x and v R = v R ˆ x ,where the indices L and R refer, respectively, to the left and right sides of the simulation box. As the two shocks Bohdan et al.
Table 1.
Simulation ParametersRuns ϕ L y ( λ si ) m i /m e ω pe / Ω e M A M s β e Eq. Eq. ∗ ∗ ∗ ∗ α = 0 . α = 0 . . o . . · − . . o
24 100 12 31 . . · − . o
12 100 17 . . · − . o . . · − . o . . · − . o . . . · − . o
12 100 12 35 . . · − . Note —Parameters of simulation runs described in this paper. Listed are: the orientation of the uniform perpendicular magneticfield with respect to the 2D simulation plane, ϕ , the transverse size of the computational box, L y , in units of the ion skin depth, λ si , the ion-to-electron mass ratio m i /m e , the plasma magnetization, ω pe / Ω e , and Alfv´enic and sonic Mach numbers, M A and M s , the latter separately for the left (runs *1) and the right (runs *2) shock. We also list the electron plasma beta, β e , for eachsimulated shock and the critical Alfv´enic Mach number (Eq. 1) for α = 0 .
2, as well as the modified trapping condition (Eq. 2)calculated for α = 0 . α = 0 . λ se = 20∆. Figure 2.
Illustration of the simulation setup. 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 homoge-neous magnetic field, B , that is perpendicular to theshock normal and lies in the yz plane. The magneticfield thus forms an angle ϕ with the y -axis. Initializedwith the flow is a motional electric field E = − v × B ,with v = v L or v = v R , respectively, for the left andthe right beam. We assume that the beams move withequal absolute velocities, v L = v R = 0 . c , and thatthe magnetic field strength in both plasmas is equal, B = B . The motional electric field thus has equalstrength and opposing signs in the two slabs. We use themethod of Wieland et al. (2016) to suppress the artifi-cial electromagnetic transient that results from the ini-tial strong electric-field gradient between the two plasmaslabs.We collide plasma beams of equal density but differ-ent temperatures, thus studying two different shocks inone simulation. The temperature ratio between the twobeams is 1000, so that the sonic Mach numbers, M s , of the two shocks differ by a factor of √ (cid:39)
30. Interms of the electron plasma beta (the ratio of the elec-tron plasma pressure to the magnetic pressure) the leftbeam has β e , L = 5 · − and the right beam β e , R = 0 . β e , R = 0 . down-stream rest frame of the two shocks.The parameters of the simulation runs described inthis paper are listed in Table 1. We have performedseven large-scale numerical experiments (runs A–G),that feature in total fourteen simulated shocks. Herewe refer to each of these shock cases as to a separatesimulation run, and tag the shocks in the left plasma( β e , L = 5 · − ) with *1, and the right shocks with *2( β e , R = 0 . ϕ = 0 o , and run G usesthe out-of-plane magnetic field orientation, ϕ = 90 o .We do not consider simulations with ϕ = 45 o , becausethe shock structure and the acceleration mechanismsobserved in this case are almost identical to those inruns with the in-plane field configuration (Bohdan et al.2017). The runs with the in-plane magnetic field covera wide range of ion-to-electron mass ratios and Alfv´enicMach numbers, as illustrated in Figure 3, which permitsan investigation of the influence of these parameters onthe electron acceleration efficiency and to scale our re- lectron SSA at SNR shocks v A = B / (cid:112) µ ( N e m e + N i m i ), where µ is the vacuum per-meability, N i and N e are the ion and the electron num-ber densities, and B is the far-upstream magnetic-fieldstrength. The sound speed reads c s = (Γ k B T i /m i ) / ,where k B is the Boltzmann constant, Γ is a nonrela-tivistic adiabatic index, and T i is the ion temperature.The Alfv´enic, M A = v sh /v A , and sonic, M s = v sh /c s ,Mach numbers of the shocks in Table 1 are given in theconventional upstream reference frame. As the in-planeand the out-of-plane magnetic field lead to a differentnumber of degrees of freedom, the adiabatic indicesare different with Γ = 5 / ϕ = 0 o and ϕ = 90 o . Thus the resulting expectedshock speeds take values v sh = 0 . c for runs A–F and v sh = 0 . c for runs G. In the simulation frame thespeeds are smaller by the shock compression ratio.To investigate the role of SSA in electron pre-acceleration, we adjust the magnetic-field strength, B ,to establish Alfv´enic Mach numbers that test the trap-ping conditions defined by Equations 1 and 2. A com-parison of the Alfv´enic Mach numbers and the massratio of all runs with trapping limits is offered in Fig- Figure 3.
The Alfv´enic Mach numbers and mass ratiosof the simulation runs. Runs A–F with in-plane magneticfield configuration are depicted with red dots. Run G withthe out-of-plane field is marked with a green dot. The bluesolid line shows the scaling given by the trapping conditionof Eq. 1, calculated for α = 0 .
2. The blue dash-dotted anddotted lines show the modified trapping condition (Eq. 2) for α = 0 . α = 0 .
5, respectively. ure 3. Nevertheless, we always consider weakly mag-netized plasmas with the ratio of the electron plasmafrequency, ω pe = (cid:112) e N e /(cid:15) m e , to the electron gyrofre-quency, Ω e = eB /m e , in the range ω pe / Ω e = 8 . − . e is the electron charge, and (cid:15) is the vacuum per-mittivity. To keep the plasma beta constant we adjustthe plasma temperatures and hence the sound speedsand resulting sonic Mach numbers (see Table 1).In this work we want to verify several hypotheses. Thefirst is the scaling of the SSA efficiency with the ion-to-electron mass ratio for shocks that fulfill the trappingcondition of Equation 1, here applied to the in-planemagnetic field configurations. Runs A, B, E, and Fdefine the set of simulations conducted for m i /m e =50 , , α ≤ .
5. The question to be addressed is whether 2Dsimulations with in-plane magnetic field configurationcan reproduce the SSA efficiency observed in 2D runswith the same m i /m e and the out-of-plane fields, heremarked as runs G.The third set of simulations consists of runs D andE, performed for the same mass ratio m i /m e = 200.The Alfv´enic Mach number in run D clearly violatesEquation 1, and so we expect a very low intensity ofBuneman waves. Nevertheless, particle acceleration canstill occur in the shock foot and ramp, whose structure isdefined by the magnetic filaments, and we are interestedin the nonthermal electron population that forms in theabsence of SSA. Note, that cross-comparison of runs Band D, and C and E can yield the mass-ratio dependencefor shocks having the same Alfv´enic Mach numbers.The electron skin depth in the upstream plasma iscommon for all runs and equals λ se = 20∆, where ∆is the size of grid cells. The ion skin depth, λ si = (cid:112) m i /m e λ se , is used here as the unit of length. Thetime scale and all temporal dependencies are given interms of the upstream ion Larmor frequency, Ω i , whereΩ i = eB /m i . The simulation time is typically t =(6 − − i , which is enough to cover at least a few shockself-reformation cycles (see Bohdan et al. 2017). Thetime-step we use is δt = 1 / ω − .The two plasma beams injected at sides of the simu-lation box are composed of an equal number of ions andelectrons, N ppc = 20. Electron and ion plasma pairsare initialized at the same locations to ensure the initialcharge-neutrality of the system. There is no escape ofparticles from the computational box, and we use injec-tion layers receding from the CD as in Bohdan et al.(2017), which helps alleviating numerical grid-Cerenkov Bohdan et al. effects and saves computational resources. The simu-lation box expands in x -direction during the run. Thefinal size of a simulation box can reach L x ≈ λ si . Thetransverse size of the simulation box, L y = (8 . − λ si ,is large enough to cover several of the magnetic fila-ments, that are typically separated by ∼ λ si , and at thesame time limits the computational expense that growsquadratically with m i /m e . The largest simulation boxof size L x × L y = (3264 × m i /m e = 400. Open boundary conditions are imposedin the x -direction and periodic boundaries are appliedin the y -direction.The numerical code we use is a 2D3V-adapted andmodified version of the relativistic electromagnetic PICcode TRISTAN (Buneman 1993) with MPI-based paral-lelization (Niemiec et al. 2008; Wieland et al. 2016) andthe option to trace individual particles. RESULTSIn Section 3.1 we describe the structure of the Bune-man wave modes in all simulations and also summarizethe findings of Bohdan et al. (2017). Then we discussthe electron acceleration efficiency through SSA in Sec-tion 3.2. 3.1.
The Buneman Instability
Figure 4 presents the maps of the electrostatic fieldamplitude in the foot of the right shocks (runs A2-G2,see Table 1), that propagate in moderate-temperatureplasmas with β e = 0 .
5. Only portions of the simulationboxes are shown to facilitate one-to-one comparison be-tween the runs. Run G2* is run G2 at a different phaseof shock reformation. The electrostatic fields are calcu-lated as | E ES | = | − ∇ φ | , where φ is the electric poten-tial, that is derived directly from the charge distribution.The maps are plotted for simulation times, at which thecyclic shock self-reformation allows the strongest Bune-man modes. Note, that the maps for runs B2 and G2can be compared with Figures 6a3 and 6c3, respectively,in Bohdan et al. (2017), in which results for runs B1 andG1 are presented (marked as runs A1 and C1, respec-tively).The properties of the Buneman instability discussed inBohdan et al. (2017) can be readily observed in Figure 4.The wave vectors are approximately parallel to the shocknormal for the in-plane configurations (runs A2-F2) andoblique for out-of-plane magnetic field (run G2). Thisreflects the motion of shock-reflected ions: for ϕ = 0 o the ions are confined to the xz -plane whereas for ϕ =90 o they stream in the simulation plane. The Bunemanwave region shows a patchy structure for the in-planefield configurations, that can be linked to clumps in the overshoot produced by merging magnetic filaments. Intotal, the Buneman waves occupy a much smaller regionthan for the out-of-plane configuration, for which thewaves are coherent and more intense.The phase velocity of the Buneman modes matchesthe relative speed between shock-reflected ions and in-coming electrons of the upstream plasma. Since for ϕ = 0 o part of the ion motion is outside of the sim-ulation grid, the wavelengths of the Buneman wavesare smaller ( λ ≈ . λ se ) than for out-of-plane field, forwhich λ ≈ . λ se . Note, that Figure 4 shows | E | andhence the wavelength is twice the separation of wavefronts, here provided in units of the ion skin depth. Thesurface area of the Buneman wave region for shocks inmoderate-temperature plasma is 20%-30% larger thanat the corresponding low- β shocks, but the intensity ofthe waves is 20%-50% smaller (compare Fig. 6 in Bo-hdan et al. 2017).For the high- β systems presented in Figure 4, Table 2lists peak amplitude of Buneman waves and the frac-tion of pre-accelerated electrons. The runs A2, B2, E2and F2 satisfy the trapping condition of Equation 1 (seeFig. 3), and both the peak and average strength of theelectrostatic field are similar. Small differences betweenthem arise from shock reformation. We conclude thatirrespective of the mass ratio, the physical conditions atshocks with M A satisfying Equation 1 are similar. How-ever, the electrostatic force is weaker in average than theLorentz force on a γ (cid:38) | E ES | / ( cB ) < Table 2.
Dimensionless peak amplitude of Buneman wavesand fraction of pre-accelerated electronsRun max( | E ES | / ( cB )) N e , BI /N e , tot (%)A2 1.1 0.43B2 1.3 0.46C2 2.3 0.6D2 0.4 0.34E2 1.3 0.49F2 1.1 0.44G2 2.7 6.8G2* 2.3 2.7 Note —For β e = 0 . | E ES | / ( cB )), of the electrostatic waves,calculated as mean | E ES | / ( cB ) for the 100 simulation cellswith the highest | E ES | / ( cB ), and the fraction of electronspre-accelerated to ( γ − > . Considerably larger electrostatic field amplitudes,reaching | E ES | / ( cB ) ∼ .
3, can be observed for run C2.Here, the Alfv´en Mach number of the shock, M A = 46, lectron SSA at SNR shocks Figure 4.
Dimensionless electrostatic field amplitudes in selected regions of the shock foot with the most intense Bunemanwaves for runs 2. The map marked as run G2* is chosen at time moment when the average field strength is the same as in runC2. is much larger than the minimum M A defined by Equa-tion 1 and also satisfies the modified trapping conditionof Equation 2, that for the measured α (cid:39) .
32 gives theminimum M A (cid:39) .
2. The field intensity in run C2 isabout a factor of 2 larger than in both run B2 with thesame mass ratio, m i /m e = 100, and run E2 with massratio m i /m e = 200 but similar Alfv´en Mach number, M A (cid:39)
45. This shows that the strength of the electro-static modes is driven by the value of the Alfv´enic Machnumber in relation to the trapping condition (Eq. 1).The absolute value of M A is not important, as in runD2 we see Buneman waves with amplitudes a factor of3 lower than those in run B2 with the same Alfv´enicMach number. Essentially all observed wave intensities are slightly weaker than the saturation level estimatedby Ishihara et al. (1980).The modified trapping condition (Eq. 2) was expectedto compensate for the effect of the field configuration.Shocks with sufficiently large M A should then repro-duce similar Buneman wave intensities in 2D in-planemagnetic field configurations than in simulations withthe out-of-plane fields. However, the electrostatic fieldin run C2 is weaker by 20% than that in run G2 with ϕ = 90 o . At a different phase of shock reformation runG2, now called G2*, has the same electric-field ampli-tude as C2, but four times the number of pre-acceleratedelectrons. This discrepancy might arise from Equation 2only compensating for the neglect of the z -motion ofions. In out-of-plane simulation we observe that the Bohdan et al. relative speed between electrons and reflected ions canreach ∼ . c , because of acceleration in upstream elec-tric field, which is a factor of ∼ . z -direction. It may be that we need toalso account for this effect by adding a factor of 1.5 tothe modified trapping condition, M A ≥ . (cid:115)
21 + sin ϕ (1 + α ) (cid:18) m i m e (cid:19) . (3)This equation gives M A (cid:39) . M A = 35 . ϕ = 0 o field configuration.This value is much larger than any of the Mach num-bers studied here for m i /m e = 100, and thus requiresattention in the future.3.2. Electron Acceleration in the Buneman Zone
Table 2 lists the fraction of electrons that have beenpre-accelerated in the Buneman wave zone to ( γ − > . N e , BI /N e , tot . This fraction is much larger in run G2than it is in runs A2-F2. Bohdan et al. (2017) arguedthat at least part of this difference is due to differencesin the amplitude of the electrostatic waves and theircoverage area. Figure 5.
Simulation-frame kinetic-energy spectra of elec-trons in the regions of the shock foot selected for Fig. 4 color-coded for run A2 (blue), run B2 (green), run E2 (red) and forrun F2 (orange). The dotted green line indicates the spec-trum of upstream cold plasma electrons (extracted from runB2).
Figures 5 and 6 show kinetic-energy spectra of elec-trons occupying the Buneman wave regions highlightedin Figure 4. Figure 5 shows energy spectra for runsA2, B2, E2 and F2, for which the Alfv´enic Mach num-bers exceed by the similar margin the trapping condition(Eq. 1). The spectra are statistically indistinguishable,and the fraction of pre-accelerated electrons is ∼ . ∼ .
34% of electrons are pre-accelerated.Although run C2 satisfies the modified trappingcondition of Equation 2, the acceleration efficiency, N e , BI /N e , tot (cid:39) .
6, is much less than for run G2. InFigure 6 the spectrum for run C2 is also compared withthe spectrum calculated for run G2 at a different phaseof the shock-reformation (denoted as run G2*), at whichthe strength of the Buneman waves matches that for runC2. Still, the fraction of pre-accelerated electrons in runG2* is four times that in run C2, but the maximum en-
Figure 6.
Spectra of electrons as in Fig. 5 for run B2 (black),run C2 (blue), run D2 (green), run G2 (red) and run G2*(orange). The dotted black line indicates the spectrum ofupstream cold plasma electrons (extracted from run B2). lectron SSA at SNR shocks Figure 7.
Interaction of electrons with Buneman waves for in-plane runs (panels (a1)-(e1), case E2) and out-of-plane runs(panels (a2)-(e2), case G2). Panels (a*): map of E x at the time indicated by the vertical black lines in the lower panels. Overlaidare the position of an electron (black dot) at the same time moment as E x maps, its trajectory history for the past 60 ω − ce andpast positions of the electron for every ω pe t = 10 intervals, designated with red dots. Panels (b*): evolution of electron energy.Panels (c*): evolution of electron momentum. Panels (d*): dimensionless components of electric field at electron position in thesimulation frame. Panels (e*): components of electric field at electron position in the electron rest frame. ergies of the electrons are comparable, max( γ ) ≈ − Bohdan et al.
Figure 7 illustrates the first stage of the SSA pro-cess for the in-plane ( left panels a1-e1) and the out-of-plane case ( right panels a2-e2). For specific electronsextracted from runs E2 and G2, we see the time evolu-tion of the energy (Fig. 7b) and the momentum (Fig. 7c),as well as the electric field at the location of the particlein the simulation frame and in the instantaneous parti-cle rest frame (Fig. 7d and e, respectively). The latteris particularly interesting, because in the electron restframe the electric field is the sole provider of accelera-tion. We refer to the selected electron in the in-planecase (left panels) as the first electron and the other oneas the second electron. Initially both electrons movewith the plasma bulk. To be trapped by electrostaticwaves, the electrons must travel with the waves againstthe upstream plasma flow, and hence be picked-up fromthe thermal pool. Before doing so, the electrons movein the negative x-direction undisturbed through severalelectrostatic wavefronts. Significant energy gain com-mences at time tω pe = 6270 for the first electron and at tω pe = 3745 for the second electron. The particles thenremain trapped by the waves and undergo the first stageof acceleration at time intervals tω pe = (6275 − tω pe = (3745 − p y momentumremains small. The end of the first-stage accelerationis marked by the black vertical line in Figure 7, beyondwhich the electrons resume gyrating.The acceleration of the first electron occurs in thesame way as in 1D geometry (Hoshino & Shimada 2002):the electron is pushed toward the upstream region by theelectrostatic field of a Buneman wave, which for sometime compensates the x -component of the Larmor ac-celeration and thus keeps the electron roughly in phasewith the wave. Consequently the average values of E x and E y electric field components are close to zero inthe particle reference frame (Fig. 7e1). The continuousgradient in p z at tω pe = (6275 − z -direction, and the energy gain terminateswhen the electron loses phase coherence with the Bune-man wave. In reality the energization will terminateearlier. In the upstream flow frame all the energy gaincomes from the field of the Buneman wave though.The second electron displays a similar behaviour, butin the frame of the obliquely propagating waves. At tω pe (cid:39) E PRF ,y ≈ x -direction to the z -direction.Let us estimate the energy gain arising from trappingat an electrostatic wavefront. Equations A5 and A8give the rate of energy gain for the out-of-plane and thein-plane configuration, respectively. The phase speedthat the electrons need to match is v ph, = 0 . c and v ph, = 0 . c for in-plane and out-of-plane configura-tion, respectively, and so ∆ v and hence the energizationrate is twice larger in the out-of-plane case than it isfor in-plane magnetic field. The total energy gain is theproduct of the rate of gain and the time of interaction.The time of interaction is limited by three factors: theintermittency of waves, escape by acceleration perpen-dicular to the wave front, and escape to the side of thewave front.In the in-plane case the wave front is infinitely ex-tended in z direction, and no escape to that side is pos-sible. For an out-of-plane magnetic field and an averagespeed along the wave front of ∼ (0 . − . c , the elec-trons would escape trapping on t esc ≈ (25 − ω − ,as the wave fronts in Figure 7 have a lateral extent ofabout 5 λ se .The escape time perpendicular to the wave front canbe estimated as t esc ≈ πω − ( v Φ + v ) /v e , WRF , where v e , WRF is the velocity of electrons in the wave frame.For the out-of-plane case this gives t esc , (cid:38) ω − , asthe average electron speed v e , WRF (cid:46) . c .The trapping time coming from the wave time inter-mittency can be estimated directly from simulations.The ability to accelerate an electron up to a certain en-ergy depends not only on the instantaneous local elec-trostatic field strength but also on the previous strengthhistory and the ability to trap an electron during thewhole acceleration period. In Figure 8 the time evo-lution of electrostatic field strength ( E ES , black dash-dotted line) at a chosen location is presented. This fieldis able to trap electrons with energies shown with blacksolid line, which is calculated assuming the equality be-tween electrostatic and Lorentz forces at the chosen lo- lectron SSA at SNR shocks Figure 8.
Schematic time evolution of the electrostatic fieldstrength (black dash-dotted line) at a chosen location in theBuneman wave rest frame. The black solid line is the max-imal energy of electron can be trapped by the electrostaticfield. Red lines represent energy histories of electrons, forwhich trapping is possible (red solid line) and impossible(red dashed line). cation. At a time t (cid:48) the electrostatic field is capable totrap an electron with energy (cid:15) . However, taking intoaccount the evolution of E ES and the energy of electrons(red dashed line) this electron cannot be trapped dur-ing the whole acceleration period. Therefore electronswith a final energy (cid:15) and the energy history shown withthe red solid line can be present in the simulation. Ac-cording to these considerations the trapping time reads t tr , ≈ ω − and t tr , ≈ ω − , (4)which is approximately the time that we analytically es-timated based on the acceleration in the direction of thewave motion (see the Appendix). Thus one of the mainlimiting factors for electron acceleration is the intermit-tency of the Buneman waves.Calculated average energy gains are∆ γ ≈ .
18 and ∆ γ ≈ . , (5)which are similar to those for electrons in Figure 7 andaverage energies of accelerated electrons in Figures 5and 6. The analytically expected energy increase can bewritten as ∆ ε ≈ e E ES | F | ∆ v t tr == ∆ v t tr ω pe m e ( m e /m i ) (1 / , (6)where e E ES = m e ∆ vω pe ( m e /m i ) (1 / (Ishihara et al.1980; Amano & Hoshino 2009b; Matsumoto et al. 2012) and | F | is assumed to be about 1. Therefore the maindifference in the acceleration rate comes from velocitydifference, ∆ v = ( v Φ + v ), and the energy gain of elec-trons is still stronger in the out-of-plane case due to alarger phase speed of the Buneman waves.We note that the modified trapping conditions (Eq. 2or 3) refer to reaching a certain strength of the electro-static field that is needed for trapping, while the energygain of electrons is related to the velocity difference be-tween reflected ions and upcoming electrons. This ve-locity difference imposes the main restriction for the in-plane simulations in their applicability to mimic realisticSSA efficiency. Using a higher Mach number can not sig-nificantly change the SSA efficiency in case of the samevelocity difference defined by the magnetic field config-uration. For the same mass ratio the number of pre-accelerated electrons is larger by about (30-40)% in theruns with a higher Alfv´enic Mach number (see Table 2,runs B2-C2 and D2-E2), which is not the factor of 10required to reach the SSA efficiency seen in out-of-plainruns. Therefore significant modifications of the parame-ters of the simulation (not just a change of the Alfv´enicMach number) are needed to reproduce the out-of-plainSSA efficiency by means of in-plane simulations.The energy difference associated with climbing or slid-ing down the potential well of a Buneman wave canbe estimated as ∆ γmc = eE BI λ BI / π . The wave-length of Buneman waves, λ BI = 2 π ∆ v/ω pe , then im-plies an energy change ∆ γ ≈ .
05 in the in-plane caseand ∆ γ ≈ .
17 in the out-of-plane run. This is insuf-ficient to redirect an incoming electron to stationarityin the wave frame. Fluctuations in the Buneman wavefield are clearly needed to trap particles and keep themin phase with the waves.We observe that for ϕ = 90 o a larger number of elec-trons are picked up from the bulk plasma for furtheracceleration than is seen with the in-plane configura-tion, which can be explained by a twice stronger e E ES force in the out-of-plane case. SUMMARY AND DISCUSSIONWe analyse electron injection processes at nonrel-ativistic perpendicular collisionless shocks with highAlfv´enic Mach numbers with 2D3V numerical PIC sim-ulations. Earlier studies indicated that SSA operates atthe leading edge of the foot as first-stage electron pre-acceleration mechanism, provided the Alfv´enic Machnumber satisfies a condition of efficient driving of theelectrostatic Buneman waves (the trapping condition,Matsumoto et al. 2012). In Matsumoto et al. (2015) andBohdan et al. (2017) we showed that in 2D simulationsthat use a field component which lies in the simulation2
Bohdan et al. plane, the downstream nonthermal-electron fraction ismuch lower than with out-of-plane mean field. Notingthat much of this difference results from an incompleteaccount of the Buneman instability in the in-plane ge-ometry, and motivated by results of recent 3D studieswhich demonstrate that the injection physics past theSSA stage can adequately be studied with 2D in-planesimulations (Matsumoto et al. 2017), here we furtherinvestigate electron acceleration by SSA at perpendic-ular high- M A shocks with in-plane magnetic field con-figurations. The aim is to infer the SSA efficiency, inparticular the validity of the trapping condition in itsoriginal form and the variant proposed in Bohdan et al.(2017), and the relation to the SSA efficiency observedin simulations with the out-of-plane fields.Our results can be summarized as follows: • The energy gain in SSA always arises from theelectrostatic field of a Buneman wave with whichthe electron travels for some time. The appar-ent acceleration, ˙ v , reflects the superposition ofelectrostatic acceleration and Larmor accelerationthat might be described as effect of the motionalelectric field in the wave frame. This processis more efficient in the out-of-plane case becauseboth the phase speed and the amplitude of thewaves are higher than for ϕ = 0 o . • As in high- M A shock simulations with out-of-plane magnetic fields, for in-plane magnetic fieldthe strength of the electrostatic wave modes in theshock foot is determined by the Alfv´enic Machnumber in relation to the trapping condition.The more M A exceeds the trapping condition,the stronger the intensity of the Buneman waves.Shocks with Alfv´enic Mach numbers satisfyingthe trapping condition by the similar margin showcomparable wave strengths in simulations for dif-ferent ion-to-electron mass ratios. • Shocks in simulations with in-plane magnetic fielddemonstrate electrostatic wave intensities lowerthan those observed in the out-of-plane case, evenif the modified trapping condition is satisfied. • The trapping time is mostly defined by intermit-tency of, and limited phase-coherence of electrons with, the Buneman waves. This limits the dura-tion of the velocity match between electrons andthe waves. • The number of electrons pre-accelerated via SSAin the shock foot strongly correlates with thestrength of the electrostatic waves. Shocks withthe same physical conditions defined through thetrapping condition show similar SSA efficiency.The latter is proportional to M A for a given massratio. However, SSA always produces larger frac-tions of pre-accelerated electrons in simulationswith the out-of-plane configurations, even if theintensities of the Buneman waves are similar asin the in-plane case. One reason for that is thelarger number of electrons being picked up fromthe bulk plasma for SSA compared to the in-planeconfiguration.We conclude that with an in-plane magnetic-field con-figuration we can not achieve the same level of SSA effi-ciency as in simulations with out-of-plane magnetic fieldor 3D simulations (Matsumoto et al. 2017), unless theparameters and settings of the simulation setup are sig-nificantly modified.This paper is conceived as the first of a series inves-tigating different aspects of electron acceleration pro-cesses at non-relativistic perpendicular shocks using PICsimulations. Interaction with Weibel filaments and mag-netic reconnection in the shock transition, plasma heat-ing, and the generation of turbulent magnetic field willbe covered in forthcoming publications.We thank the anonymous referee for their com-ments. The work of J.N. has been supported by Nar-odowe Centrum Nauki through research project DEC-2013/10/E/ST9/00662. This work was supported byJSPS-PAN Bilateral Joint Research Project Grant Num-ber 180500000671. The numerical experiment was pos-sible through a 10 Mcore-hour allocation on the 2.399PFlop Prometheus system at ACC Cyfronet AGH. Partof the numerical work was conducted on resources pro-vided by the North-German Supercomputing Alliance(HLRN) under projects bbp00003 and bbp00014. lectron SSA at SNR shocks A. ANALYTICAL MODEL OF ELECTRON SSAA.1.
Out-of-plane configuration, ϕ = 90 o In the simulation frame, the large-scale magnetic field of the right plasma slab, B = B ˆ z , induces a motional electricfield, E = − v B ˆ y , where v is the speed of the upstream plasma flowing in − x direction. The entire Larmor orbitof 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 Θ to the x-axis. The electric field carried by the wave is E x = E ES F cos Θ E y = E ES F sin Θ , (A1)where the wave factor is F = sin (cid:18) ω pe v Φ + v [ x cos Θ + y sin Θ − v Φ t ] + Φ (cid:19) . (A2)Here we allow for an arbitrary phase, Φ. The phase speed of the wave, v Φ , 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, ω pe = k ( v Φ + v ).Now consider an electron with velocity components v x and v y . Using non-relativistic kinematics we find the accel-eration of the electron as ˙ v x = − Ω e E ES B F cos Θ − Ω e v y ˙ v y =Ω e v + Ω e v x − Ω e E ES B F sin Θ . (A3)Let us rotate the coordinate system by an angle Θ, so that x (cid:48) is oriented in the direction of motion of the waves and y (cid:48) is perpendicular to it. The corresponding accelerations then read˙ v x (cid:48) =Ω e (cid:18) v sin Θ − E ES B F − v y (cid:48) (cid:19) ˙ v y (cid:48) =Ω e ( v cos Θ + v x (cid:48) ) . (A4)The wave factor, F , is explicitly time-dependent and may induce rapidly oscillating acceleration. The other termsonly describe Larmor gyration in the flow frame and hence no real energy gain. The wave factor must be approximatelyconstant, if continuous energy gain is to be achieved for about 10 plasma times, ω − , as observed. This requires thaton average v x (cid:48) − v Φ (cid:46) . c or roughly acceleration from v x (cid:48) = 0 . c to v x (cid:48) = 0 . c , after which the electron is out ofphase with the wave and commences Larmor motion.The Larmor motion of the reflected ions mandates a wave direction for which sin Θ is negative. Likewise, the wavefactor, F , must be negative to effect energy gain. Equation A4 then indicates that acceleration in y (cid:48) direction followsthat in x (cid:48) direction, and for a fair range of initial conditions ˙ v y (cid:48) is slightly less than ˙ v x (cid:48) and increases with the same rate,at least for up to 1 Ω − (cid:39) ω − . Correspondingly, the momentum component p x increases approximately linearly,and the increase in speed is approximately E ES / (2 B ), whereas p y remains approximately constant.The effective acceleration toward the upstream region arises from the superposition of acceleration in the electrostaticfield of the Buneman waves and the Larmor acceleration, that are oppositely directed in y direction, but both havepositive components in x direction. Acknowledging that both F and sin Θ must be negative, the rate of energy gainin the upstream flow frame is m ddt ( v x + v ) + v y eE ES | F | [( v x + v ) cos Θ − v y | sin Θ | ] (A5)and hence completely independent of the motional electric field. In the simulation frame the velocity component v cos Θ disappears from Equation A5 and a new component of energy-gain rate appears, Ω e v v y , which captures theapparent energy by Larmor motion in this frame.4 Bohdan et al.
A.2.
In-plane configuration, ϕ = 0 o The main impact of the in-plane configuration is that the part of the Larmor motion is perpendicular to the simulationplane, and so the orientation and properties of the Buneman waves are modified, as only wave vectors in the simulationplane can be captured. The wave factor changes to F = sin (cid:18) ω pe v Φ + v [ x − v Φ t ] + Φ (cid:19) . (A6)The acceleration then follows by appropriate rotation of that given in Equation A3,˙ v x = − Ω e E ES B F + Ω e v z ˙ v z = − Ω e v − Ω e v x . (A7)Obviously, there is linear acceleration in − z direction, if the particle can be held at approximately constant phase( F < v x ≈ v Φ ) in the wave. As v Φ (cid:38) v it is the Larmor acceleration that is responsible for the particle’s slidingalong the wavefront, and the electrostatic field of the waves provides slow energy gain at a rate m ddt ( v x + v ) + v z eE ES | F | ( v x + v ) , (A8)which also only involves the electrostatic field of the Buneman waves. The energy gain will be less than that forout-of-plane configuration, because only part of the motion of the back-streaming ions can drive waves that hence havelower amplitude, E ES , and additionally the velocity term in Equation A8 is reduced.REFERENCES Aharonian, F. A. 2013, Astrop. Phys., 43, 71Amano, T., & Hoshino, M. 2007, ApJ, 661, 190Amano, T., & Hoshino, M. 2009, ApJ, 690, 244Amano, T., & Hoshino, M. 2009, Physics of Plasmas, 16,102901Axford, W. I., Leer, E., & Skadron, G. 1977, InternationalCosmic Ray Conference, 11, 132Blandford, R., & Eichler, D. 1987, PhR, 154, 1Bocchino, F., Orlando, S., Miceli, M., & Petruk, O. 2011,A&A, 531, A129Bohdan, A., Niemiec, J., Kobzar, O., & Pohl, M. 2017,ApJ, 847, 71Buneman, O. 1993, in
Computer Space Plasma Physics:Simulation Techniques and Software , Eds.: Matsumoto &Omura, Tokyo: Terra, p.67Buneman, O. 1958, Phys. Rev. Lett,1,8Dahlin, J. T., Drake, J. F., & Swisdak, M. 2014, Physics ofPlasmas, 21, 092304Dahlin, J. T., Drake, J. F., & Swisdak, M. 2015, Physics ofPlasmas, 22, 100704Drury, L. O. 1983, Reports on Progress in Physics, 46, 973Hoshino, M., Mukai, T., Terasawa, T., & Shinohara, I.2001, J. Geophys. Res., 106, 25979Hoshino, M., & Shimada, N. 2002, ApJ, 572, 880 Ishihara, O., Hirose, A., & Langdon, A. B. 1980, PhysicalReview Letters, 44, 1404Kato, T. N., & Takabe, H. 2010, ApJ, 721, 828Matsumoto, Y., Amano, T., & Hoshino, M. 2012, ApJ, 755,109Matsumoto, Y., Amano, T., & Hoshino, M. 2013, PhysicalReview Letters, 111, 215003Matsumoto, Y., Amano, T., Kato, T. N., & Hoshino, M.2015, Science, 347, 974Matsumoto, Y., Amano, T., Kato, T. N., & Hoshino, M.2017, Physical Review Letters, 119, 105101Niemiec, J., Pohl, M., Stroman, T., & Nishikawa, K.-I.2008, ApJ, 684, 1174-1189Niemiec, J., Pohl, M., Bret, A., & Wieland, V. 2012, ApJ,759, 73Northrop, T. G. 1963, Reviews of Geophysics and SpacePhysics, 1, 283Oka, M., Fujimoto, M., Shinohara, I., & Phan, T. D. 2010,Journal of Geophysical Research (Space Physics), 115,A08223Oka, M., Phan, T.-D., Krucker, S., Fujimoto, M., &Shinohara, I. 2010, ApJ, 714, 915Petruk, O., Dubner, G., Castelletti, G., et al. 2009,MNRAS, 393, 1034 lectron SSA at SNR shocks15