Non-thermal particle acceleration at highly oblique non-relativistic shocks

The non-thermal acceleration of electrons and ions at an oblique, non-relativistic shock is studied using large scale particle-in-cell (PIC) simulations in one spatial dimension. Physical parameters are selected to highlight the role of electron pre-heating in injection and subsequent acceleration of particles at high Mach number shocks. Simulation results show evidence for the early onset of the diffusive shock acceleration process for both electrons and ions at a highly oblique subluminal shock. Ion acceleration efficiencies of $\lesssim 5\%$ are measured at late times, though this is not the saturated value.


INTRODUCTION
High Mach number, non-relativistic shocks are known to be strong sources of non-thermal radiation (e.g. Hinton & Hofmann 2009). Multi-wavelength measurements are thought to be a signature of diffusive shock acceleration (DSA) (see Blandford & Eichler 1987;Bell 2014, for reviews). Crucially, DSA requires frequent scattering on electric/magnetic field fluctuations to maintain near isotropy of the particle distribution. If particles interacting with the shock can escape the thermal pool in sufficient quantity, the mechanism is self-sustaining (Eichler 1979). The plasma processes that initiate injection into DSA are not yet fully understood, though it is clear the ambient conditions, in particular the magnetic field obliquity, play an important role (e.g. Kirk & Heavens 1989;Caprioli & Spitkovsky 2014;Guo & Giacalone 2015). If the shock propagation and magnetic field directions are uncorrelated, it follows that the mean obliquity, i.e. the angle between the upstream magnetic field and shock normal direction, is θ n = 60 • . For a given shock velocity υ sh , particles can not outrun the shock upstream along fieldlines if | cos θ n | < υ sh /c (so-called superluminal shocks). However, for non-relativistic shock velocities this transition occurs at large angles: θ n 90 • . Thus characterising the ability of shocks with θ n > 60 • , (approximately half of all shock surfaces) to inject electrons and ions from the thermal pool and subsequently accelerate to high energies is of interest to the high-energy astrophysics community. That highly oblique shocks can accelerate particles via DSA is not controversial. γ-ray observations of the binary stellar system η-Carinae provide compelling evidence that both ion and electron acceleration are ongoing at the oblique shocks produced in the wind collision region (White et al. 2020). The relative contribution of electrons and ions to the γ-ray emission from supernova remnants is unclear. However, the X-ray synchrotron rims that outline almost the entire circumference of young remnants such as Cassiopeia A (Patnaude & Fesen 2014) and Tycho (Eriksen et al. 2011) indicate not only that electrons are accelerated locally to TeV energies, but that acceleration occurs over a large range of obliquities. A counter-example is SN 1006 for which non-thermal emission is more prominent in the quasiparallel regions (the "polar caps") of the supernova remnant (Reynoso et al. 2013). However, the reasons underlying this asymmetry are not conclusive (Winner et al. 2020), and may result from a combination of effects related to the ambient plasma conditions, the acceleration process itself or an inability to amplify magnetic fields sufficiently, and not necessarily a consequence of injection efficiency (see for example Bell et al. 2011;Reville & Bell 2013).
The physical mechanisms that determine the injection of electrons and ions into the acceleration process at magnetised shocks are distinct. The large difference in their respective masses imply that the relevant wavesparticle resonances fall in different frequency ranges. In particular, for electrons to participate in DSA, a preheating mechanism to energies ∼ m i u 2 sh /2, is required to facilitate resonance with low frequency MHD waves. This problem is often referred to as the "injection problem". Various mechanisms such as electrostatic Buneman instability, oblique whistler waves or Weibel instability, have been proposed as potential solutions (Amano & Hoshino 2009;Muschietti & Lembège 2017;Amano & Hoshino 2010;Riquelme & Spitkovsky 2011;Lembege 1990;Oka et al. 2017;Bohdan et al. 2019a,b).
In this Letter we explore the role of electron and ion kinetic scales separation in a realistic system, while capturing a sufficiently large number of ion-gyration times to demonstrate the onset of non-thermal ion acceleration. The large number of particles-per-cell (ppc) needed to counter the energy loss caused by plasma drag forces (Kato 2013;May et al. 2014) restricts the current study to 1D3V (one dimension in space and three in momentum) simulations. In the following we address the aforementioned issues and study both electron and ion acceleration at an oblique shock. The numerical set-up is summarised in the next section. Results are presented in section 3, while our main conclusions are discussed in section 4.

PIC SIMULATION SETUP
For the simulations presented, we use the open-source PIC code SMILEI (Derouillat et al. 2018). We employ a relaxation method of shock excitation by initializing the simulation in the shock-normal frame (SNF). The Rankine-Hugoniot (RH) jump conditions determine the downstream parameters (superscript d) for our desired upstream conditions (superscript u) (Tidman & Krall 1971;Leroy et al. 1981;Umeda et al. 2014). Unless stated otherwise, time and length are normalized to the inverse of the electron plasma frequency (ω u pe = 4πn u e e 2 /m e ) and skin-depth (c/ω u pe ), where n u e = 1 is the upstream electron plasma density, c is velocity of light in vacuum, m e (−e) is the electron mass (charge). Velocity, momentum, temperatures and fields are normalized to c, m e c, m e c 2 and m e ω u pe c/e, respectively. Making use of the following relations for the plasma beta (ratio of thermal to magnetic pres- , the upstream and downstream parameters are easily expressed in dimensionless unit, again using the RH jump conditions. Here are the Alfvénic and sonic Mach numbers, respectively. The ion mass is denoted by m i , and T u e , T u i are the electron and ion temperatures. The drift velocity of upstream plasma (or equivalently the shock velocity) is in accordance with previous studies on collisionless shocks (Umeda et al. 2014;Leroy et al. 1981). Different downstream temperature ratios were explored and early time results were comparable. Previous numerical shock studies (Leroy et al. 1981;Park et al. 2015;Kato 2015;Fang et al. 2019;Hanusch et al. 2020) have focused on the regimes of high and low plasma betas. The resulting wave analyses in these regimes are more amenable to theoretical investigation. Here we focus on the intermediate regime β ∼ 1, since extreme limits are in general specific to given astrophysical environments, which is beyond the present scope. The generation of plasma waves is sensitive to β, and consequently the electron heating mechanisms are expected to show a complex character in the intermediate plasma beta regime. We choose the upstream magnetic field, B u = B u xx + B u zẑ making an angle of θ n = 75 • with the shock normal. We assume a reduced ion to electron mass ratio m i /m e = 256, and grid and time resolutions of ∆x = 0.1 and ∆t = 0.095, respectively. At t = 0, the shock discontinuity is initialized with a thickness (in normalized units) δx = 2ρ i = 2(υ u /B d )(m i /m e ) with hyperbolic tangent profile. The total length of the simulation domain is L x ≈ 570ρ i . As mentioned previously, plasma drag can inhibit electron energization if ppc number is insufficient (Kato 2013;May et al. 2014). Using Eq.(37) of Kato (2013), we conservatively estimate a minimum requirement for ppc to be 100. Figure 1. Zoomed-in view of the ion density map and the transverse magnetic field. Time and the length are normalized to ion gyro frequency (ωci = (me/mi)B d ) and ion gyro radius (ρi) respectively. The ion number density is normalised to its initial upstream value, and similarly the magnetic field is normalized to the initial upstream value, B u z .
Consequently, we use 512(1024) ppc per species for the upstream (downstream) population. The global energy conservation was satisfied to better that 1% throughout the simulation.
3. RESULTS Figure 1 shows the spatio-temporal evolution of the plasma density and transverse (B z ) magnetic field. It can be seen that the initial magneto-hydrodynamic (MHD) discontinuity evolves into a well-defined shock structure, separating the upstream and downstream plasmas, and the creation of a shock overshoot region. The transverse magnetic field follows a similar pattern to the plasma density. Periodic magnetic field reformation occurs at the shock at approximately the ion gyrofrequency (corresponding to the downstream magnetic field) (Umeda et al. 2014;Leroy et al. 1981;Lyu & Kan 1990;Lembege 1990;Burgess & Scholer 2015). Magnetic field reformation is generally associated to insufficient dissipation in the shock ramp region leading to ion reflection. In the case of quasi-parallel shocks for example, this results in the generation of high-frequency whistler waves (Lyu & Kan 1990). The transverse magnetic field seen in Fig.1 shows strong filamentary structure in the downstream region. In both panels, one may also notice (for ω ci t > 30), the appearance of sharp features extending in the upstream region. These structures are generated in the shock foot region, excited by the reflection of incoming ions, which are however both electromagnetic and electrostatic in nature, in contrast to expectations of quasi-parallel shocks (Lyu & Kan 1990). These waves play an important role in the energization of the electrons as we show later in Fig.4. Both pe and pi are normalized to mec. The insets show the evolution of energy in the non-thermal populations of the respective particle species with time. Note ηe,i is the ratio of integrated kinetic energy in the non-thermal particles relative to thermal particles. Fig.2 shows the evolution of the momenta distribution for both species in the downstream. The spectra is summed from x i ≈ 2ρ i to x f ≈ 16ρ i in the downstream direction and have been multiplied by p 2 to highlight the emergence of flat non-thermal tails at high energies. At the end of the simulation, an appreciable fraction of electrons have acquired sufficient momentum for their gyro radii to exceed the shock thickness (in the range p e ∼ 10−40). These electrons can participate in the DSA process. Once this threshold is reached (at ω ci t ≈ 50) the emergence of a distinct population above this momentum is visible, although this appears to be a transient effect. The growth of a non-thermal ion population is also evident, although the separation from the thermal population is not as developed as that of the electron population, due to the limited number of gyro-periods captured. The onset of the non-thermal ion population develops in earnest only after the nonthermal electrons have acquired gyro-radii comparable to that of the bulk of the ions, emphasising the necessity for long run times. This is also indicated by features in the time evolution of the acceleration efficiency η e,i , defined here as the ratio of the kinetic energy in the thermal and non-thermal components. For convenience we define the non-thermal component to be the integrated kinetic energy density of particles with p e > 4 in the case of electrons and p i > 90 for ions (see Figure 2). At the end of the simulation, the acceleration efficiency of electrons, according to our definition, is η e ≈ 45% while for the ions it is η i ≈ 5% (see insets). The large value for the non-thermal electron energy efficiency is the combined effect of increasing maximum energy, and depletion of the thermal component. Relative to the incoming upstream kinetic energy density, this efficiency would remain small. The ion energy efficiency on the other hand is non-negligible, and crucially shows no sign of saturating at the end of the simulation. Fig.3 shows the trajectories of the five most energetic electrons and ions from the simulation corresponding to Fig.2. Particles are tracked only after crossing a population specific threshold condition. As noted earlier, the ion acceleration is negligible prior to electrons acquiring sufficient energies to have their gyro-radii comparable to that of the drifting ions. A causal connection between the two however, has not been established. To reach these energies, a continuous energization process (or processes) is required, which we address below (see Fig.4). Such pre-heating/energization mechanisms are unnecessary for ions, since their gyro-radii are already of the order of the shock thickness. The highest energy ions are seen to make large excursions into the upstream direction and display features of DSA acceleration. The guiding centre for example undergoes diffusive like behaviour in the turbulent pre-and post-shock regions. Note that by selecting the five highest energy particles, we are biased towards particles that have remained close to the shock, and therefore undergone the largest number of acceleration cycles. At the end of the simulation, many energetic particles penetrate more than a gyroradius upstream. Particles that venture far upstream have large negative z-component of their velocity, as required to escape along the field-lines.
To explore the electron energization, including the effects of ion to electron mass ratio, we tracked the evolution of two populations of test electrons in simulations with different mass ratios. These test electrons are always present in the full simulation run and evolve according to the self-consistent time-dependent electromagnetic fields. However, they do not contribute to the charge and current densities used to update these fields. Each population of test electrons contains 1000 particles, and are initialized with the same plasma parameters as upstream electrons viz. charge, mass, temperature, and velocity. Each group is distributed in a span of δx = 1000 in the upstream region, at a fixed distance from the shock. Fig.4 reveals the dependence of electron energization (γ factor) on the mass ratio, a consequence of the larger cross-shock potential. One may also note a few additional trends: first the popu-  lation of electrons that arrive later in time (upper row) generally display a larger energy gain compared to the population that arrives earlier (lower row). Energization is due to both the cross-shock potential and the interaction with the high frequency electrostatic and electromagnetic waves traveling upstream, as has been previously discussed (Leroy et al. 1981;Ligorini et al. 2021;Umeda et al. 2014;Muschietti & Lembège 2017). Since the Alfvén Mach number in our simulations satisfies the condition M A ≥ (cos θ/2) (m i /m e )β u e , generation of Buneman waves and oblique whistler modes are expected (Amano & Hoshino 2010). These waves have ample time to grow in magnitude and transfer energy to the later population of electrons (upper row). Since the first group of test electrons interact with these waves earlier in the process, they do not gain as much energy. The second trend is that, superimposed on the gyromotion of the electrons, pitch-angle scattering features are observable. The effect is most prominent in the large mass-ratio simulations at later times (upper-left most panel) indicating that, similar to the cross-shock potential, the strength of the waves is dependent on the mass-ratio.

CONCLUSION AND DISCUSSION
We have investigated oblique, subluminal nonrelativistic shocks using PIC simulations, and found they are capable of accelerating both electrons and ions. At the end of our simulation, a moderate energyconversion efficiency ∼ 5% is found for the ions, although this is not its saturated value. Early-time results from the simulation are broadly consistent with those found in previous studies in a comparable regime (e.g. Xu et al. 2020), using a different shock initialization approach (piston-driven shock in the upstream rest-frame). Several previous works find that electron acceleration is efficient for oblique high-Mach-number shocks, though the rate of acceleration decreases above p e ≈ u sh m i /m e . We verify this trend (see Fig.3) that electrons are "rapidly" energized to the same momentum, at which point they can interact with ion-scale fluctuations. However, as the gyration periods (a proxy for the acceleration time) above this momentum are now comparable, longer simulation times reveal that the spectral cut-offs march approximately in tandem, a feature not previously discussed. At the end of the simulation, the ions are trans-relativistic, and thus their increased inertia (and gyration period) result in less frequent shock-crossings. Acceleration to energies beyond what is achieved here thus require significantly longer run-times, which is not currently feasible. This presents a challenge to determining the saturated ion acceleration efficiency. While the ion efficiencies presented are still low compared to predictions from quasi-parallel shocks configurations (e.g. Caprioli & Spitkovsky 2014;Park et al. 2015), the most energetic particles continue to penetrate upstream. In this case, we expect the acceleratedparticle mediated precursor to continue growing in extent and amplitude with time. If this holds, we speculate that the conditions close to the viscous shock (i.e. the position of the MHD discontinuity) become increasingly turbulent, with particle transport in the shock vicinity being largely independent of the initial field orientation (cf. Reville & Bell 2013).
Identification of a causal connection between ion and electron acceleration has not been established. In our simulations, a transition in the acceleration efficiency appears soon after the maximum energy electrons have comparable gyro-radius to the incoming ions. This indicates that the acceleration of the two species are intimately linked. Whether further electron acceleration bootstraps ion acceleration, or vice-versa, is unclear. The initial energization of electrons from the thermal pool appears to be sensitive to the phase at which they first interact with the (reforming) shock. Favourable electrons are reflected and can subsequently interact with the high-frequency waves in the shock-foot region, where they undergo further heating. The stopping and reflection of the ions at the shock ramp is responsible for driving these high frequency electrostatic and electromagnetic perturbations in the upstream incoming plasma (Amano & Hoshino 2010). We confirm that the interactions of test electrons with these waves in the immediate shock upstream is the primary channel for heating. Noticeable differences are visible between test electrons that interact early in the shock evolution, compared with those that arrive later, when these waves are further developed. It is not possible however to attribute the initial electron energization to one single instability or a single physical process for our simulation parameters.
The generation of waves in the upstream region depends not only on the Alfvén/sonic Mach numbers, and mass ratio, but at early times appears also to be sensitive to the initial setup of the numerical simulations, which we discuss below (Amano & Hoshino 2010, 2012Bohdan et al. 2017;Muschietti & Lembège 2017;Umeda et al. 2014;Oka et al. 2017;Riquelme & Spitkovsky 2011;Xu et al. 2020;Ligorini et al. 2021). The threshold Alfvén Mach number for triggering these waves plays an important role in astrophysics, particularly in systems with variable shock conditions. White et al. (2020) recently interpreted the quenching of electron heating close to periastron pas-sage of the binary system η-Carinae to occur when the shock Mach number falls below the critical whistler Mach number M A,crit = m i /m e , although different predictions can be found in the literature; e.g. Amano & Hoshino (2010) suggest M A = (cos θ/2) (m i /m e )β u e , while conversely others find whistler waves favour low Mach shocks M A m i /m e (see the 2D simulations of Riquelme & Spitkovsky 2011, and references therein). Our current simulations exceed both thresholds, although a parameter survey similar to that of Xu et al. (2020) is needed. Irrespectively, our results have important consequences for the interpretation of astrophysical observations.
In this work we focus on simulations in the SNF (see Burgess & Scholer 2015;Lembege 1990, for a summary of different shock initialization methods). Whether simulation results specific to non-thermal particle acceleration are sensitive to the shock launching mechanism is an open question, as the initial development is unavoidably sensitive to this choice. Previous simulations in the downstream frame, where the incoming plasma is reflected off a wall, the shock formation is mediated by different beam-plasma instabilities (Riquelme & Spitkovsky 2011;Kato 2015). Nevertheless, our intermediate time results bear striking similarity with those of Xu et al. (2020), who perform simulations in the upstream rest-frame. This suggests the long-term evolu-tion of the shock dynamics and particle acceleration is the same in different frames using alternative setups. A more direct comparison can be made between the present results, performed in the SNF, with matching simulations (using the relaxation setup) performed in the de-Hoffman and Teller frame (dHT) (Kumar & Reville 2021). This frame is ideal for acceleration studies since both the shock and the zeroth-order field lines are static (in the SNF, field lines are continuously dragged across the shock). These two frames are related by an appropriate Lorentz boost along the shock surface. While the runtimes for results presented in Kumar & Reville (2021) were shorter, compatible spectra are found at these times, when corrected for frame transformations. 1 However, in Lorentz transforming, the 1D simulation probes a different k-vector of the upstream plasma dispersion. This motivates the necessity for simulations in more than one dimension (cf. Riquelme & Spitkovsky 2011). Thus, further large-scale higher dimensional simulations are required to shed more light on the results presented in this Letter.