# Diffusive Shock Acceleration at Oblique High Mach Number Shocks

Allard Jan van Marle, Artem Bohdan, Paul J. Morris, Martin Pohl, Alexandre Marcowith

► To cite this version:

Allard Jan van Marle, Artem Bohdan, Paul J. Morris, Martin Pohl, Alexandre Marcowith. Diffusive Shock Acceleration at Oblique High Mach Number Shocks. Astrophys.J., 2022, 929 (1), pp.7. 10.3847/1538-4357/ac5962. hal-03644494

HAL Id: hal-03644494

https://hal.science/hal-03644494

Submitted on 13 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Diffusive Shock Acceleration at Oblique High Mach Number Shocks

Allard Jan van Marle<fnref n="1" />, Artem Bohdan<fnref n="2" />, Paul J. Morris<fnref n="2" />, Martin Pohl<fnref n="2" />,<fnref n="3" />, and Alexandre Marcowith<fnref n="1" />

Received 2021 October 25; revised 2022 February 25; accepted 2022 February 27; published 2022 April 8

## Abstract

The current paradigm of cosmic-ray (CR) origin states that the greater part of galactic CRs is produced by supernova remnants. The interaction of supernova ejecta with the interstellar medium after a supernova's explosions results in shocks responsible for CR acceleration via diffusive shock acceleration (DSA). We use particle-in-cell (PIC) simulations and a combined PIC-magnetohydrodynamic (PIC-MHD) technique to investigate whether DSA can occur in oblique high Mach number shocks. Using the PIC method, we follow the formation of the shock and determine the fraction of the particles that gets involved in DSA. With this result, we use PIC-MHD simulations to model the large-scale structure of the plasma and the magnetic field surrounding the shock and find out whether or not the reflected particles can generate upstream turbulence and trigger DSA. We find that the feasibility of this process in oblique shocks depends strongly on the Alfvénic Mach number, and the DSA process is more likely to be triggered at high Mach number shocks.

*Unified Astronomy Thesaurus concepts: Shocks (2086); Supernova remnants (1667); Plasma astrophysics (1261); Interstellar medium (847)*

## 1. Introduction

Collisionless shocks are commonly observed when plasmas move with supersonic velocities, for example, in planetary systems, supernova remnants (SNRs), and jets of active galactic nuclei. In such flows, collective particle interactions form a shock layer on kinetic plasma scales much smaller than the collisional mean free path. In the shock transition, a part of the bulk kinetic energy is converted into energies of thermal particles and electromagnetic fields through wave-particle interactions. The microphysics of these processes is still not fully understood, even though it has recently been the object of several numerical studies, especially in the high Mach number regime (e.g., Marcowith et al. 2020; Pohl et al. 2020). This regime will be partly covered in this study.

Diffusive shock acceleration (DSA), a process entering the category of first-order Fermi acceleration, is the process by which astrophysical shocks accelerate charged particles to relativistic speeds. DSA requires the magnetic field near the shock front to reflect the particle, leading to repeated shock-crossings with the particle gaining energy at each crossing (e.g., Bell 1978; Blandford & Ostriker 1978; Drury 1983). DSA is a self-sustaining process because the presence of high-energy particles triggers instabilities in the magnetic field, which in turn allow the magnetic field to reflect the particles more efficiently (e.g., Bell 1978, 2004). DSA has been explored numerically using the particle-in-cell (PIC) method, as well as the PIC-hybrid method, whereby the ions are treated as particles and the electrons as a fluid. Simulations of this type seem to indicate that although this process is effective in the case of (quasi-)parallel shocks where the magnetic field is aligned with the direction of motion, it becomes ineffective

once the angle between the magnetic field and the shock exceeds approximately 50° (Caprioli & Spitkovsky 2014a, 2014b, 2014c).<fnref n="4" /> However, different results were obtained by van Marle et al. (2018), using a combined PIC-MHD approach. This method, which treats the thermal gas as a fluid but nonthermal ions as individual particles, showed that, given a sufficient injection rate, even for magnetic obliquity up to 70°, shocks can accelerate particles through the DSA process. Several factors were believed to contribute to this discrepancy between PIC-MHD and PIC or PIC-hybrid methods:

1. Time-period covered by the simulation. The DSA occurring in PIC-MHD simulations is only effective after more than 200 ion-gyrotimes ($\omega_{ci}^{-1}$) and it takes more than 500 ion-gyrotimes to fully develop whereas published hybrid simulations where limited to scales of the order of 200 ion-gyrotimes. Notice that at the time $\sim 200 \omega_{ci}^{-1}$ both PIC-hybrid and PIC-MHD simulations find similar results, namely that particles are accelerated by the shock-drift acceleration (SDA).

2. Effective number of particles per cell. DSA requires an upstream current of CRs in order to develop. Because only a small fraction of the particles that cross the shock is reflected, a large particle population is required to produce sufficient nonthermal particles in the upstream medium.

3. Simulation box size. van Marle et al. (2018) showed that the relevant upstream instabilities operate on long

<figure><img src="image_2.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Laboratoire Univers et Particules de Montpellier (LUPM) Université Montpellier, CNRS/IN2P3, CC72, Place Eugène Bataillon, F-34095 Montpellier Cedex 5, France</footnote>

<footnote n="2">²Deutsches Elektronen-Synchrotron DESY, Platanenallee 6, D-15738 Zeuthen, Germany; artem.bohdan@desy.de</footnote>

<footnote n="3">³Institute of Physics and Astronomy, University of Potsdam, D-14476 Potsdam, Germany</footnote>

<footnote n="4">⁴ Notice that most of the studies on DSA including the present one assume shocks are propagating into a homogeneous medium; effects of a nonmonotonous shock are discussed in Hanusch et al. (2019). Most simulations also assume that the pre-shock medium is fully ionized. Simulations for a partially ionized medium (Ohira 2013, 2016) show an increased injection rate, allowing for DSA even in oblique shocks. Also, studies by Kumar & Reville (2021) demonstrate that DSA can be triggered by the shock emitted waves at highly oblique shocks.</footnote>
wavelengths, larger than the box size of hybrid simulations.

In order to address these discrepancies, Haggerty & Caprioli (2019) present an improved version of the PIC-Hybrid code including relativistic ions dynamics. Long simulations, with larger box sizes and equivalent nonthermal particle statistics as in PIC-MHD simulations were performed for Mach number 30 shocks with a magnetic obliquity of 70°. The authors confirm that DSA does not develop at this obliquity angle. This leaves a final possibility: the injection rate of nonthermal particles. Because the PIC-MHD method treats the thermal plasma as a fluid, it cannot model the process by which particles start to deviate from thermal equilibrium as they pass through the shock. Instead, it relies on an ad-hoc description by which a fixed fraction of the gas crossing the shock is converted into nonthermal particles. van Marle et al. (2018) assumes that the injection rate for a high-obliquity shock was equal to that of a quasi-parallel shock. However, as studies shown in Caprioli & Spitkovsky (2014a) and this paper shows, this is not the case. The injection rate for high-obliquity shocks is substantially lower than for a quasi-parallel shock with a similar Mach number. The question is whether this lower fraction is sufficient to trigger the kind of upstream instability that is required to initiate DSA. Another aim is to test, using a full PIC approach, the injection rate obtained with the PIC-Hybrid method as its level has a strong impact on the development of DSA and streaming instabilities. An analytical approximation of the relevant timescales in Appendix A of van Marle (2020) shows that the ability of particles to distort the upstream magnetic field depends, among other factors, on the local particle density. If it is too low, the thermal gas and the magnetic field can resist the force that the particle are exerting on the field. Notice that this is not the same as the criterion for triggering the streaming instability. It is possible to distort the magnetic field without causing the thermal gas to become unstable. Instead, the magnetic field is distorted by the presence of nonthermal particles through the Lorentz force, which reduces the particle energy. Subsequent waves of particles are then accelerated by the distortions at the expense of reducing the distortion of the magnetic field. This creates a cyclic pattern by which one wave of particles distorts the field while losing energy and the next smooths the field while accelerating. A similar pattern was observed for high-$\beta$ shocks at sonic Mach numbers below 2.85 (van Marle 2020).

In this paper, we further explore the parameter space of high-obliquity shocks to investigate which, if any, semi-perpendicular shocks can generate a sufficient number of nonthermal particles to trigger upstream distortions of the magnetic field and to eventually initiate DSA. We do so by performing 2D PIC simulations to determine the injection rate of supra-thermal particles at the shock as a function of shock velocity, Mach number, and obliquity. Once the injection rate is properly defined, we perform corresponding 2D PIC-MHD simulations to follow the large-scale, long-term evolution of the plasma and the magnetic field near the shock. The combination of the two techniques is used to cover a wider range of plasma beta and its impact on the shock acceleration efficiency. Then, PIC-MHD simulations are used to investigate high Alfvénic Mach number shocks, a regime barely addressed in hybrid simulations.

We present a short description of simulation setup in Section 2. The results are presented in Section 3. The

discussion and conclusions are given in Sections 4 and 5, correspondingly.

## 2. Numerical Methods and Setups

### 2.1. PIC Simulations

PIC simulations are performed using the modified version of the relativistic electromagnetic TRISTAN code (Buneman 1993) with MPI-openMP hybrid parallelization (Niemiec et al. 2008) and the particle sorting optimization (Dorobisz et al. 2018). The Vay solver (Vay 2008) is used to update particle positions. The triangular-shape cloud particle shapes (the second-order approximation) and Friedman et al. (1990) filter for electric and magnetic fields are used to suppress the numerical grid-Cherenkov short-wave radiation. The numerical model used in our simulations has been extensively tested, and all individual upstream energy components are conserved to better than 5% accuracy over the entire time of simulation.

We initialize shocks using the reflecting-wall setup, where the left boundary is reflective, the right boundary is open, and periodic boundary conditions are applied in the y-direction. The electron-ion upstream plasma moves in the negative x-direction with the absolute velocity $v_0$ (Figure 1(a)). It collides with the conducting wall, and further interaction of the reflected and the upstream plasma flows provides a shock propagating in positive x-direction. To prevent an initial transient caused by this artificial contact discontinuity between the reflecting wall and the plasma slab, we introduce a drift current to ions close to the reflecting wall in the same way it is discussed in Wieland et al. (2016).

The large-scale upstream magnetic field, $B_0$, makes an angle $\theta_B$ with the shock normal (or positive x-direction) and an angle $\varphi$ with the simulation plane. Simulations where $\varphi = 0°$ and $\varphi = 90°$ are referred to as in-plane and out-of-plane, correspondingly. Shock simulations are performed in a 2D3V configuration, which follows two spatial coordinates and all three components of the particle velocities and electromagnetic fields. Therefore the adiabatic index is $\Gamma = 5/3$ and the expected compression ratio is about 4. The resulting shock speed either equals $v_{sh}' = 0.067c$ or $0.033c$ in the simulation (downstream) reference frame and $v_{sh} = 0.263c$ or $0.133c$ in the upstream reference frame.

The Alfvén velocity is defined as $v_A = B_0/\sqrt{\mu_0\rho}$ with $\rho = N_e m_e + N_i m_i$, where $\mu_0$ is the vacuum permeability, $m_i$ and $m_e$ are the ion and the electron masses, and $N_i$ and $N_e$ are the ion and the electron number densities, which are 40 particles per cell for each species. The sound speed reads $c_s = (\Gamma k_B T_i / m_i)^{1/2}$, where $k_B$ is the Boltzmann constant and $T_i$ is the ion temperature that initially is equal to the electron temperature defined as $k_B T_i = m_i v_{th}^2 / 2$, where $v_{th}$ is defined as the most probable speed of the upstream plasma particles in the upstream reference frame. The plasma beta is $\beta_p = 2\mu_0 N_i k_B T_p / B^2$. The Alfvénic, $M_A = v_{sh} / v_A$, and sonic, $M_s = v_{sh} / c_s = \sqrt{2/\Gamma\beta_p} M_A$, Mach numbers of the shocks are defined in the conventional upstream reference frame. We list in Table 1 the parameters $\beta_p$ and $M_A$, then $M_s$ can easily be deduced.

The ratio of the electron plasma frequency, $\omega_{pe} = \sqrt{e^2 N_e / \epsilon_0 m_e}$, to the electron gyrofrequency, $\Omega_e = eB_0 / m_e$, is in the range of $\omega_{pe} / \Omega_e = 10.6–31.9$. Here, $e$ is the electron charge, and $\epsilon_0$ is the vacuum permittivity. The electron skin depth in the upstream plasma is constant for all runs and equals
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Sketch of the setups for PIC (a) and PIC-MHD (b) shock simulations.</figcaption></figure>

$\lambda_{se} = 8\Delta$, where $\Delta$ is the size of a grid cell. The ion skin depth is defined as $\lambda_{si} = \sqrt{m_i/m_e} \lambda_{se}$ and $m_i/m_e = 50$ for all PIC simulations. The spatial scale is given in terms of the ion upstream gyro-radius, $r_{gi} = m_i v_{sh} / (eB_0)$. The simulation time-step is $\Delta t = 1/16 \omega_{pe}^{-1}$. The timescales are given in terms of the upstream ion Larmor frequency, $\Omega_{ci} = eB_0/m_i$. The simulation time is $T_{sim} \approx 30\Omega_i^{-1}$ for all simulations, except run PIC15 where $T_{sim} \approx 68\Omega_i^{-1}$. The simulation time of PIC simulations is limited by the computational expense, especially at high Alfvénic Mach numbers and low shock velocities. Also, if a simulation is too long, artificial heating can substantially change the upstream conditions.

## 2.2. PIC-MHD Simulation Setup

The PIC-MHD setup is different from the PIC model. In this model, the simulations are run in the frame of reference of the shock rather than the frame of reference of the downstream medium, using a similar approach as used in van Marle et al. (2018) and van Marle (2020). The setup is initialized by starting from the analytical solution for a standing shock as derived through the Rankine–Hugoniot conditions for a magnetized thermal gas. At the start of the simulation, the grid is filled with a thermal gas according to the specific shock conditions of that simulation. Once the simulation starts, nonthermal particles are introduced at the shock front, according to the injection rates derived from the corresponding PIC simulations. The particles are given an absolute velocity of twice the pre-shock velocity of the thermal gas and an isotropic velocity distribution in the post-shock medium. Because we use an isotropic velocity distribution in the post-shock rest frame, we inject twice the number of nonthermal particles. This gives us the desired reflected particle fraction. As will be shown by the PIC results (see Section 3.1), there is a difference in injection rates between simulations with the magnetic field in the plane of the simulation, and those where the magnetic field is out of the plane. This indicates that the 2D results are not fully representative of the 3D reality. Therefore, we choose an injection rate in between the two results two reflect the 3D nature of the problem.

We use a 2D simulation box<fnref n="5" /> that covers 240 × 30 $r_g$, where $r_g$ is the gyro-radius defined by the particle injection velocity and the upstream magnetic field. The grid has a base resolution of 240 × 30 grid cells with the nested OCTREE adaptive mesh allowing for three additional levels of refinement, which gives

us an effective grid of 1920 × 240 cells.<fnref n="6" /> The thermal gas flows in from the outer boundary of the x-axis, passes through the shock halfway along the x-axis, and leaves through the inner x-boundary. Both x-boundaries are transparent to particles, allowing the particles to leave the simulation when they reach them. The y-boundaries are periodic for both thermal gas and particles. The magnetic field in the thermal gas makes an angle $\theta_B$ with the flow in the plane of the simulation. In order to investigate the occurrence of DSA, we need to be able to follow the particles as they move along the field lines. Therefore, the PIC-MHD simulations have to be performed with the magnetic field in the plane of the simulation. We have performed one PIC-MHD simulation with the field out-of-plane for comparison to the PIC simulations. The simulation parameters for the different setup are summarized in Table 2.

Common runs between PIC and PIC-MHD techniques are PIC-MHD01B and PIC07. The PIC-MHD code does not handle relativistic MHD, hence it cannot handle flows with speeds much in excess of $0.1c$; this explains why most of the simulations are done for $V_{sh} = 0.05c$. However, as we discuss below, the injection rates are not strongly dependent on the shock speed at nonrelativistic speeds. A strong dependence does develop at relativistic speeds because of the energy extra requirements for acceleration at such velocity. We have also developed a special relativistic version of the code but it requires much more numerical resources, and the results will be presented elsewhere.

# 3. Results

## 3.1. Results of PIC Simulations

### 3.1.1. Shock Structure

Figure 2 displays a magnetic field map and the phase-space plot of $p_{x\uparrow-x}$ at the shock region of run PIC15. A portion of the upstream ions are accelerated via multiple SDA cycles and reflected back upstream (Caprioli et al. 2015), where they can potentially be further accelerated via DSA. We henceforth denote these ions as injected. The fraction of injected ions is about 2%, which is in rough agreement with Caprioli’s hybrid simulation (Caprioli et al. 2015) for the same shock parameters. The injected ions propagating upstream drive transverse waves with a wavelength of about $\lambda \approx 2r_{gi}$ (Figure 2) and an amplitude of about $\delta B/B_0 \approx 0.1$ at $x = (9 - 10)r_{gi}$, consistent with the expected wavelength of resonant modes that are

<footnote n="5">⁵ van Marle et al. (2019a) performed a 3D study of particle acceleration at shocks.</footnote>

<footnote n="6">⁶ Note that an insufficient length of the simulation box in the direction of flow may change the outcome of the simulation. See Appendix B for further explanation.</footnote>
<table><caption>**Table 1**
PIC Simulation Parameters</caption><thead><tr><th>Name</th><th>θ<sub>B</sub></th><th>In/out-of plane</th><th>M<sub>A</sub></th><th>V<sub>sh</sub>/c</th><th>β<sub>p</sub></th><th>ω<sub>pc</sub>/Ω<sub>e</sub></th><th>N<sub>inj</sub>/N<sub>0</sub></th><th>U<sub>inj</sub>/U<sub>sh</sub></th></tr></thead><tbody><tr><td>PIC01</td><td>65</td><td>In-plane</td><td>20</td><td>0.133</td><td>1</td><td>21.2</td><td>0</td><td>0</td></tr><tr><td>PIC02</td><td>65</td><td>Out-of-plane</td><td>20</td><td>0.133</td><td>1</td><td>21.2</td><td>2.8 × 10<sup>-6</sup></td><td>5.8 × 10<sup>-5</sup></td></tr><tr><td>PIC03</td><td>60</td><td>In-plane</td><td>20</td><td>0.264</td><td>1</td><td>10.6</td><td>2.8 × 10<sup>-6</sup></td><td>6.4 × 10<sup>-5</sup></td></tr><tr><td>PIC04</td><td>60</td><td>Out-of-plane</td><td>20</td><td>0.264</td><td>1</td><td>10.6</td><td>1.3 × 10<sup>-5</sup></td><td>2.5 × 10<sup>-4</sup></td></tr><tr><td>PIC05</td><td>60</td><td>Out-of-plane</td><td>30</td><td>0.264</td><td>1</td><td>15.9</td><td>4.3 × 10<sup>-5</sup></td><td>5.5 × 10<sup>-4</sup></td></tr><tr><td>PIC06</td><td>60</td><td>Out-of-plane</td><td>30</td><td>0.133</td><td>1</td><td>31.9</td><td>9 × 10<sup>-5</sup></td><td>1.1 × 10<sup>-3</sup></td></tr><tr><td>PIC07</td><td>60</td><td>In-plane</td><td>20</td><td>0.133</td><td>1</td><td>21.2</td><td>8.4 × 10<sup>-6</sup></td><td>1.3 × 10<sup>-4</sup></td></tr><tr><td>PIC08</td><td>60</td><td>Out-of-plane</td><td>20</td><td>0.133</td><td>1</td><td>21.2</td><td>6.5 × 10<sup>-5</sup></td><td>8.9 × 10<sup>-4</sup></td></tr><tr><td>PIC09</td><td>60</td><td>Out-of-plane</td><td>20</td><td>0.133</td><td>10</td><td>21.2</td><td>3.7 × 10<sup>-5</sup></td><td>5.7 × 10<sup>-4</sup></td></tr><tr><td>PIC10</td><td>60</td><td>Out-of-plane</td><td>20</td><td>0.133</td><td>0.1</td><td>21.2</td><td>3.3 × 10<sup>-5</sup></td><td>4.5 × 10<sup>-4</sup></td></tr><tr><td>PIC11</td><td>55</td><td>In-plane</td><td>20</td><td>0.133</td><td>1</td><td>21.2</td><td>4 × 10<sup>-4</sup></td><td>4.2 × 10<sup>-3</sup></td></tr><tr><td>PIC12</td><td>55</td><td>Out-of-plane</td><td>20</td><td>0.133</td><td>1</td><td>21.2</td><td>4.1 × 10<sup>-4</sup></td><td>4 × 10<sup>-3</sup></td></tr><tr><td>PIC13</td><td>50</td><td>In-plane</td><td>20</td><td>0.133</td><td>1</td><td>21.2</td><td>2.7 × 10<sup>-3</sup></td><td>0.02</td></tr><tr><td>PIC14</td><td>50</td><td>Out-of-plane</td><td>20</td><td>0.133</td><td>1</td><td>21.2</td><td>2.6 × 10<sup>-3</sup></td><td>0.018</td></tr><tr><td>PIC15</td><td>45</td><td>In-plane</td><td>10</td><td>0.133</td><td>1</td><td>10.6</td><td>0.016</td><td>0.086</td></tr><tr><td>PIC16</td><td>45</td><td>In-plane</td><td>20</td><td>0.133</td><td>1</td><td>21.2</td><td>0.013</td><td>0.068</td></tr><tr><td>PIC17</td><td>45</td><td>Out-of-plane</td><td>20</td><td>0.133</td><td>1</td><td>21.2</td><td>0.022</td><td>0.12</td></tr></tbody></table>

**Note.** Setup of the PIC simulations. Listed are the shock obliquity, in-plane or out-of-plane configuration, the shock speed, the Alfvénic Mach number, $M_A$, and the plasma parameter $\beta_p$. The ion-to-electron mass ratio is equal to 50 for all PIC simulations. The ion injection density and energy density fractions are reported in the last two columns (see Section 3.1.2 for details).

excited by ions streaming along the large-scale magnetic field with $v_\parallel \approx 0.3c$. Caprioli & Spitkovsky (2014b) also found upstream waves of resonant nature, corresponding to a wavelength $\lambda \lesssim 2\pi r_{gi}$. However, we note that conditions at the shock upstream are also suitable for the nonresonant Bell mode (Bell 2004, 2005). According to Equations (16) and (17) from Reville et al. (2006) the growth time and the wavelength of the fastest growing mode are $t_{Bell} \approx 5\Omega_i^{-1}$ and $\lambda_{Bell} \approx 3r_{gi}$ in the asymptotic limit. Unfortunately, the simulation PIC15 is still too short to reach the steady state of the shock evolution and say which mode is dominant in our case.

This upstream turbulence is observed only in run PIC15. This run is characterized by the strongest injected ion flow, and the simulation time is long enough to disturb the upstream magnetic field. The DSA process, however, is still not initialized. Usually that requires a considerably longer simulation time, $T_{sim} > 100\Omega_i^{-1}$, hence the PIC-MHD simulations in our study (see Section 3.2).

### 3.1.2. The Ion Injection Efficiency

To study the ability to drive upstream turbulence, we use PIC simulations to determine the properties of injected ions as a function of the shock parameters. Table 1 lists the density and energy fraction of the injected ions upstream of the shock. The energy fraction is defined as $U_{inj}/U_{sh}$, where $U_{inj}$ is the energy density of injected ions in the upstream reference frame and $U_{sh} = N_0 m_i v_{sh}^2 / 2$. $U_{inj}$ and $N_{inj}$ are calculated over the region $x - x_{sh} = (2-7)r_{gi}$, where $x_{sh}$ is the shock overshoot position, and averaged in time over at least over $5\Omega_i^{-1}$. The fluctuation amplitude for $U_{inj}$ and $N_{inj}$ is at the level of 20%–30%. Note also that the injection efficiency remains stable in the longest run (PIC15), although waves are already evident upstream of the shock.

The fraction of injected ions is mainly defined by the shock obliquity angle, $\theta_B$, and the angle $\phi$. Figure 3 shows the ion energy distributions in the upstream region for all simulations having $v_{sh}/c = 0.133$, $M_A = 20$, and $\beta_p = 1$. The spectra are calculated in the upstream reference frame, and so the upstream

ions are represented by Maxwell distributions predefined by the simulation setup. The injected ions have an energy above unity, on account of SDA. In all simulations, the escaping ions have a similar average speed along the x-axis, which is about $1.6v_{sh}$ in the reference frame of the upstream plasma. To reach a such speed for higher $\theta_B$, the ions must have high velocity along the magnetic field and a high energy. In other words, at highly oblique shocks, the ions must have experienced a large number of SDA cycles, which results in smaller injection rates. Shock simulations performed with out-of-plane magnetic field are usually characterized by larger fractions of injected ions compared to in-plane cases, because the overshoot is more coherent and provides better conditions for SDA.

The ion injection efficiency is different for simulations with $v_{sh}/c = 0.264$ and $v_{sh}/c = 0.133$, the other parameters being the same (runs PIC03-PIC05 and PIC06-PIC08). In both cases, the ions have approximately the same escape speed along the x-axis and it equal $1.6v_{sh}$ if calculated in the upstream reference frame. Therefore, the total velocity is $v_{inj} \approx 1.6v_{sh}/\cos\theta_B = 3.2v_{sh}$ and the average energy of ions propagating back upstream can be estimated as $\varepsilon_{inj} \approx 12\varepsilon_{sh}$ in slower shocks (runs PIC06-PIC08), while for faster shocks (runs PIC03-PIC05), it is $\varepsilon_{inj} \approx 23\varepsilon_{sh}$. This difference is due to relativistic correction, therefore some extra energy (+130% for faster shocks and +20% for slower shocks) compared to the fully nonrelativistic case is needed to reach the same runaway speed. To achieve higher energies, ions should go through a larger number of SDA cycles, and therefore fewer ions are injected. In fact, cases with $v_{sh}/c = 0.133$ are very close to fully nonrelativistic shocks, as the relativistic corrections are about 20% which is commensurate with the variance of the reflected ion fraction in PIC simulations. Therefore, the results can be directly applied for SNR shocks.

Shock simulations which differ only by the upstream plasma beta $\beta_p$ (runs PIC08-PIC10) do not demonstrate a dependence on this parameter. The ion injection efficiency varies at the level just slightly above the random fluctuations while $\beta_p$ is probed over two orders of magnitude.
<table><caption>**Table 2**
PIC-MHD Simulation Parameters</caption><thead><tr><th>Name</th><th>M<sub>A</sub></th><th>θ<sub>B</sub>(deg)</th><th>V<sub>sh</sub>/c</th><th>N<sub>inj</sub>/N<sub>0</sub></th><th>Special physics</th><th>Notes</th></tr></thead><tbody><tr><td>PICMHD01A</td><td>20</td><td>60</td><td>0.05</td><td>5.0 × 10<sup>-5</sup></td><td></td><td>No DSA</td></tr><tr><td>PICMHD01B</td><td>20</td><td>60</td><td>0.1</td><td>2.5 × 10<sup>-5</sup></td><td></td><td>No DSA</td></tr><tr><td>PICMHD02</td><td>20</td><td>60</td><td>0.1</td><td>5 × 10<sup>-5</sup></td><td></td><td>No DSA</td></tr><tr><td>PICMHD03</td><td>20</td><td>60</td><td>0.1</td><td>1 × 10<sup>-4</sup></td><td></td><td>DSA</td></tr><tr><td>PICMHD04</td><td>30</td><td>60</td><td>0.05</td><td>5.0 × 10<sup>-5</sup></td><td></td><td>Start of DSA</td></tr><tr><td>PICMHD05A</td><td>50</td><td>60</td><td>0.05</td><td>5.0 × 10<sup>-5</sup></td><td></td><td>DSA</td></tr><tr><td>PICMHD05B</td><td>50</td><td>60</td><td>0.05</td><td>5.0 × 10<sup>-5</sup></td><td>Out-of-plane</td><td>No DSA</td></tr><tr><td>PICMHD06</td><td>100</td><td>60</td><td>0.05</td><td>5.0 × 10<sup>-5</sup></td><td></td><td>DSA</td></tr><tr><td>PICMHD07</td><td>300</td><td>60</td><td>0.05</td><td>5.0 × 10<sup>-5</sup></td><td></td><td>DSA</td></tr><tr><td>PICMHD08</td><td>50</td><td>60</td><td>0.05</td><td>5.0 × 10<sup>-5</sup></td><td>β = 0.03</td><td>DSA</td></tr><tr><td>PICMHD09</td><td>50</td><td>60</td><td>0.05</td><td>5.0 × 10<sup>-5</sup></td><td>β = 0.1</td><td>DSA</td></tr><tr><td>PICMHD10</td><td>50</td><td>60</td><td>0.05</td><td>5.0 × 10<sup>-5</sup></td><td>β = 10</td><td>DSA</td></tr><tr><td>PICMHD11</td><td>50</td><td>60</td><td>0.05</td><td>5.0 × 10<sup>-5</sup></td><td>β = 30</td><td>DSA</td></tr><tr><td>PICMHD12</td><td>50</td><td>55</td><td>0.05</td><td>4.0 × 10<sup>-4</sup></td><td>β = 30</td><td>DSA</td></tr></tbody></table>

**Note.** List of setups of PIC-MHD simulations. The columns show the Alfvénic Mach number, $M_A$, the magnetic field obliquity, $\theta_B$, the shock speed in units of the speed of light, and the cosmic-ray injection fraction. All simulations have an in-plane magnetic field and $\beta = 1$, unless specified otherwise. N.B., in PIC-MHD, the particles are injected isotropically in the post-shock rest frame. Therefore we inject $2N_{inj}/N_0$ to get the intended number of reflected particles.

PIC simulations featuring shocks with different Alfvénic Mach numbers, $M_A$, demonstrate somewhat higher fluctuations of the ion injection efficiency. The results of simulations PIC06/PIC08 and PIC15/PIC16 are consistent within 30% error bars. The injection efficiency for simulations PIC06 and PIC08 with higher shock velocity ($v_{sh} = 0.264c$) is a factor of three different. At such fast and oblique shocks, the average injection velocity of ions is almost relativistic (~0.85c) and even small fluctuations may result in a factor of a few different injection efficiencies. We also admit that when the ion injection is defined by SDA, the injection efficiency should not depend on $M_A$ (see also Section 3.1.3).

The PIC simulation results suggest that the variable that affects the injection efficiency the most is the shock obliquity angle, $\theta_B$. Indeed, changing the shock obliquity angle from $\theta_B = 45^\circ$ to $\theta_B = 65^\circ$ leads to a sharp drop in the ion injection efficiency by five orders of magnitude, while other parameters change the injection efficiency by a factor of a few at most.

### 3.1.3. An Ion Injection Model

As discussed in Caprioli et al. (2015), ions are accelerated through multiple SDA cycles and with sufficiently high velocity can escape the shock and become injected ions. During an SDA cycle, an ion is reflected back upstream with conservation of the kinetic energy in the shock rest frame. Then gyrating in the upstream region, the ion is accelerated by the motional electric field. The energy increment for cycle $n$ can be estimated as

$$ \Delta\epsilon_{SDA,n} = eE_0\Delta l_n, \qquad (1) $$

where $E_0 = v_{sh}B_0 \sin \theta_B$, $\Delta l_n = \alpha r_{g,n} = \alpha p_{1,n}/(eB_0)$ and $\alpha$ defines the acceleration efficiency of an individual SDA cycle. Therefore the ion kinetic energy after ($n+1$) SDA cycles can be defined as

$$ \epsilon_{SDA,n+1} = \epsilon_{SDA,n} + \alpha c v_{sh} m_i \sqrt{\gamma_n^2 - 1} \sin \theta_B. \quad (2) $$

For each $\theta_B$, we calculate how many SDA cycles are needed to achieve the escape speed of $v_x = 1.6v_{sh}$ in the upstream

reference frame. When the escape speed is fixed the ion injection efficiency is defined by the balance between the escape probability $p$ and the acceleration efficiency $\alpha$. The fraction of ions escaping downstream after each SDA cycle was derived by Caprioli et al. (2015) as $p = 0.75$. Therefore $\alpha$ must be equal to 0.34 for the injection model to match the simulation results. The thermal energy of injected ions is calculated assuming that the x-projection of the thermal velocity is comparable to the velocity difference between the shock and escaping particles, namely, $v_{th} \approx 0.6v_{sh}/\cos \theta_B$. The resulting energy density fractions are marked by blue dots in Figure 4. The dashed line is a fit,

$$ \frac{U_{\text{inj}}}{U_{\text{sh}}} = F(\theta_B) = 10^{a_0 + a_1 \theta_B + a_2 \theta_B^2}, \quad (3) $$

where $a_0 = -7.3$, $a_1 = 0.35$, and $a_2 = -0.0047$. Note that the escape speed in PIC simulations slightly depends on the shock obliquity: it increases from $1.55v_{sh}$ to $1.69v_{sh}$ when the shock obliquity changes from $\theta_B = 45^\circ$ to $\theta_B = 65^\circ$. Hence the escape probability or SDA acceleration efficiency might also slightly depend on the shock obliquity. Taking into account that the final results only weakly depend on these vagaries, here for simplicity we assume that all these parameters are constant.

The red dots in Figure 4 shows results for runs with $v_{sh}/c = 0.133$, $M_A = 20$, and $\beta_p = 1$. Each point represents the average for in-plane and out-of-plane configurations for the $\theta_B$ in question. The error bars include both temporal variations and the difference between in-plane and out-of-plane runs. Figure 4 demonstrates a good match between PIC simulation results and the injection model. The green line represents the critical Mach number, $M_{A,cr} = (F(\theta_B))^{-0.5}$, at which the energy density of injected ions equals that of the magnetic field in the upstream region, and the injected ions are able to drive a considerably stronger upstream turbulence necessary for efficient DSA.

For all shocks with $M_A \gg 20$, we shall assume that the ion injection rate does not depend on the Alfvénic Mach number, because the ion injection process is driven by the shock potential which is the same (when it is normalized to the
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Magnetic field map (top panel) and phase-space distribution of $p_{x_i}-x$ (bottom panel) for run PIC15 at the last simulation step, $t \approx 68 \Omega_i^{-1}$.</figcaption></figure>

upstream bulk energy) for all high Mach number shocks. Also, the energy gain during one SDA cycle does not depend on the Alfvénic Mach number (see Equation 2). At perpendicular high Mach number shocks, the physics of which is very similar to that of oblique shocks, the ion behavior does not depend on $M_A$ for a wide range of Mach numbers, $M_A \approx 20–70$ in Bohdan et al. (2020) and $M_A \approx 130$ in Kato & Takabe (2010).

Summarizing, we can state that as long as SDA works in the usual regime, the ion acceleration efficiency should remain roughly constant and the injection rate is largely defined by $\theta_B$, while the influence of the plasma beta of Alfvénic Mach number is rather minor, if existent at all. Therefore, the best fit to the injection model (the dashed blue line in Figure 4) can be used to represent the ion injection energy fraction in MHD-PIC simulation.

## 3.2. Results of PIC-MHD Simulations

### 3.2.1. Low Mach Numbers and $\beta_p = 1$ Simulations

The results obtained in Section 3.1 show that the injection rate strongly depends on the angle between the magnetic field and the flow, $\theta_B$, and this, in turn, directly influences the evolution of the upstream medium. We have run a series of simulations with an obliquity angle $\theta_B = 60^\circ$ and with Alfvénic Mach number $M_A = 20$ (PICMHD01A-03), based on the PIC results to test at which injection rate DSA starts to occur. The results of these simulations are presented in Figure 5, which shows the spectral energy distribution (SED) of the nonthermal particles at $t = 3000 \ \Omega_{\text{inj}}^{-1}$ with $\Omega_{\text{inj}}$ in gyrofrequency at injection defined by the upstream magnetic field $B_0$ and the injection velocity $v_{\text{inj}}$. For simulations with an injection rate of $5 \times 10^{-5}$ we see no spectral tail with constant slope, and hence no significant DSA. Some of the particles are accelerated, but this is due to the SDA process as was demonstrated in van Marle et al. (2018) and discussed in the previous section. This process creates the “bump” at $\gamma - 1 = 0.04$. The simulation with an injection rate of $2 \times 10^{-4}$ does show an extension toward higher energies. Its most prominent feature is a secondary bump at $\gamma - 1 = 0.1$, indicating that some of the particles that escape from the shock are reflected back toward the shock to go through a second round of shock-drift acceleration. Furthermore, instead of falling off immediately,

this SED shows a nearly constant slope at higher energies, indicating that at least some DSA occurs. The slope has an index that corresponds with $dN/d\gamma \propto (\gamma - 1)^{-3.4}$. This is much steeper than expected for a typical high Mach shock (approximately -1.5). However, the SDA process remains incomplete due to the low amplitude of the upstream instabilities. As a result, many particles escape without being reflected, which leads to a steeper spectrum. However, we find that this required injection rate is higher than what is obtained through PIC simulations. Therefore we can conclude that at $M_A = 20$, DSA does likely not occur in shocks with an obliquity of 60° or more.

### 3.2.2. High Alfvénic Mach Number Runs

In a second series of runs (PICMHD04-07) we inject flows at high Alfvénic Mach numbers ($M_A \ge 30$). This regime is not accessible to PIC simulations (at least for a reasonable computational expense). We do this by changing the magnetic field strength. The other parameters remain the same, including the injection rate, which is fixed at $5 \times 10^{-5}$, based on the results of PIC07 and PIC08, because the PIC simulations showed that there is no strong dependence of the injection rate on the Mach number. This prediction is also born out by (Caprioli & Spitkovsky 2014a, Figure 3), who argue that at high obliquity the injection rate becomes nearly independent of the shock Mach number. PIC models by Ha et al. (2018) show an increase in efficiency for increased Mach numbers, but this was for low Mach number shocks ($M_s < 3.5$) for which the compression ratio increases with the Mach number. This is not the case in our models, which all pertain to high Mach number shocks.

As is shown in Figure 6, the shape of the SED changes strongly between $M_A = 20$ and $M_A = 30$. The spectrum develops a high-energy tail that extends further for higher Mach numbers. This matches the prediction obtained from Figure 4 that, for Alfvénic Mach numbers higher than the critical one (for an obliquity angle of 60° the critical Alfvénic Mach number is around 30), the injection rate would be sufficient to trigger instabilities in the upstream medium.

For comparison, we include a simulation (PICMHD05B) with magnetic field going out of the plane of the simulation, rather than lying in the plane. Figure 7 demonstrates that the simulation with the magnetic field out of the 2D plane shows
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The upstream ion spectra for runs with $v_{sh}/c = 0.133$, $M_A = 20$, and $\beta_p = 1$. The left panel shows in-plane simulations, the right panel presents out-of-plane runs. Spectra are calculated for a region $x - x_{sh} = (2-7)r_{gi}$ ahead of the shocks and are averaged in time over at least $5\Omega_i^{-1}$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** The normalized energy density of injected ions in the upstream region: red dots are PIC simulation results (the error bars include both temporal variations and the difference between in-plane and out-of-plane runs), blue dots represent the ion injection model based on ion acceleration via multiple SDA cycles, and the dashed blue line is the best fit to the injection ion model, $F(\theta_B)$ (Equation 3). The green line gives the critical Alfvénic Mach number, $M_{A,cr} = (F(\theta_B))^{-0.5}$; see right axis.</figcaption></figure>

no sign of DSA. Under these circumstances, the particles cannot trigger an instability since their direction of motion is perpendicular to the plane of the simulation.

### 3.2.3. 55° versus 60°

We compare the results for shocks with obliquity of 55° and 60° (PICMHD05A and PICMHD12) in Figure 8 to demonstrate the transition between shocks that still follow the general model of the quasi-parallel shocks and the truly oblique shocks. The lower obliquity corresponds to a much higher injection rate (see simulations PIC11 and PIC12), for which the issue of DSA is not in doubt. However, the two spectra show a marked difference. The shock with 55° obliquity lacks the “bump” at $(\gamma - 1) = 0.03$, which is clearly visible at 60°. This bump is created by shock-drift acceleration, which is less effective at lower obliquity. Instead, the SED shows only DSA from $(\gamma - 1) = 0.01-0.05$. At these injection rates, the velocity and Alfvénic Mach number of the shock become largely irrelevant because the injected particle density is sufficient to trigger DSA for all nonrelativistic shocks.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** SEDs for simulations PICMHD01B, -02 and -03 with $M_A = 20$ and varying injection rates at time $t = 3000 \Omega_{inj}^{-1}$. At the lower injection rates, we see only the bump at $\gamma - 1 = 0.04$ that is caused by shock-drift acceleration. For the highest injection rate, we see a secondary bump at $\gamma - 1 = 0.1$ populated by particles that experienced a second round of acceleration. At higher energies, the SED shows a constant slope, indicating that at least some particles have been further accelerated. This marks the very beginning of DSA.</figcaption></figure>

### 3.2.4. Investigation of Plasma Beta Effects

The effect of varying plasma beta is demonstrated in Figure 9, which shows the SEDs for simulations PICMHD08-11 compared to the model with the same Alfvénic Mach number and $\beta_p = 1$ (PICMHD05A). Here we vary the plasma beta by changing the thermal pressure and keeping the magnetic field constant. Again, as in Section 3.2.2, we do not vary the injection rate. PIC simulations of shocks with high plasma-$\beta$ show that these shocks have no problem producing nonthermal particles as long as the sonic Mach number remains above 3.5 (Ha et al. 2018). In our models, we do not reach this point and all shocks can be considered to be strong.

Clearly, the variation in $\beta_p$ has little influence on the SED, although there is a slight trend to have less efficient DSA for higher $\beta_p$. This trend can be explained by the fact that at higher plasma beta the thermal pressure tends to counteract variations in the magnetic field, thereby reducing the local distortion. This
<figure><img src="image_8.png" /><figcaption>**Figure 6.** SEDs for simulations with $\beta_p = 1$ and varying Alfvénic Mach numbers at time $t = 3000 \Omega_{\text{inj}}^{-1}$. At $M_A = 20$, we see no significant acceleration through DSA. At $M_A = 30$ a high-energy tail starts to develop, indicating DSA, similar to what is observed for $M_A = 20$ with an exaggerated injection rate. This feature becomes more prominent at higher Alfvénic Mach numbers.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** SEDs for simulations PICMHD05A and -05B at time $t = 3000 \Omega_{\text{inj}}^{-1}$, distinguishing the difference between a magnetic field in and out of the plane of the simulation. Clearly, DSA does not take place when the magnetic field is out-of-plane.</figcaption></figure>

effect was demonstrated for extreme high-$\beta_p$ shocks by van Marle et al. (2019b) and van Marle (2020).

### 3.2.5. Resonant vs. Nonresonant Instability

As already stated in Section 3.1.1, shocks with low Alfvénic Mach number can also trigger the resonant streaming instability over short timescales. This branch has no restriction over the magnetic energy density of the injected ions with respect to the background magnetic energy. Following the discussion in Caprioli & Spitkovsky (2014b: their Equation 9) the ratio of the resonant to the nonresonant growth rate can be written as

$$ \frac{\Gamma_{\text{res}}}{\Gamma_{\text{NRH}}} \approx \frac{1}{\sqrt{N_{\text{inj}}/N_0}} \frac{1}{M_A}. \qquad (4) $$

To derive this expression, we have considered that the injection fraction scales as $1/p$. For numerically manageable injection rates,

<figure><img src="image_10.png" /><figcaption>**Figure 8.** SEDs at time $t = 3000 \Omega_{\text{inj}}^{-1}$ for simulations PICMHD05A and -12 with $\theta_B = 60^{\circ}$ and $55^{\circ}$. Both show clear evidence of DSA, but the $55^{\circ}$ model lacks the “bump” created by shock-drift acceleration.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** SEDs for different plasma beta with $M_A = 50$ at time $t = 3000 \Omega_{\text{inj}}^{-1}$. All SEDs are very close together, indicating that $\beta_p$ has little or no influence on the SED for a given Alfvénic Mach number.</figcaption></figure>

the resonant branch dominates over the nonresonant one only at low Alfvénic Mach numbers, for which the standard formula for the nonresonant growth rate violates the condition $kr_L \gg 1$. In that case, the resonant modes should be seen, if $N_{\text{inj}}/N_0$ is large enough to provide a few growth cycles upstream of the shock. At high Alfvénic Mach numbers the resonant mode is dominant only for very low $N_{\text{inj}}/N_0$ ratios. This does not mean that DSA cannot be supported by the resonant instability in the high-Mach-number regime, but that it would require too long a simulation to be demonstrated. The lower $N_{\text{inj}}/N_0$, the lower are the growth rates and the saturation field strength for both branches of the instability, and the longer is the acceleration timescale (Equation 8 in Caprioli & Spitkovsky 2014b).

Hence, considering in parallel the critical Alfvénic Mach constraint in Section 3.1.3 it is difficult at high shock obliquities to isolate a parameter subspace where the resonant unambiguously dominates the magnetic turbulence spectrum especially close or even below the critical Alfvénic Mach number. This issue deserves a dedicated future study.
## 4. Discussion

PIC-MHD simulations by van Marle et al. (2018) showed that DSA could occur at oblique shocks, provided that the injection rate of nonthermal particles was sufficiently large, and that the energy density of the upstream current exceeded that of the local magnetic field. However, as our PIC simulations show, the injection rate is strongly dependent on the angle between the flow and the magnetic field, with injection rates dropping rapidly for shocks with $\theta_{B_n} > 50$. They also show that the injection is, fundamentally, a 3D process, as demonstrated by the difference in injection rates between in-plane and out-of-plane simulations.

We find that, for $\theta_{B_n} = 60$, the injection rate is insufficient for shocks at $M_A = 20$. However, for shocks with a higher Alfvénic Mach number, the energy density of the injected upstream particles is still sufficient to trigger local instabilities in the magnetic field, thereby producing DSA. In contrast to the strong influence of the Alfvénic Mach number, the plasma beta of the upstream medium proves to be mostly irrelevant, provided that all shocks are still effectively high Mach number shocks.

We can conclude that high Mach number oblique shocks are capable of accelerating particles through DSA, whereas, for slower shocks, acceleration becomes impossible at angles of more than approximately 50°–55°, which matches previous results (e.g., Caprioli & Spitkovsky 2014a, 2014b, 2014c; Haggerty & Caprioli 2019).

Shocks with high Alfvénic Mach numbers are quite common in astrophysical objects. Stellar winds can reach velocities of several thousands kilometers per second (e.g., Krtička & Kubát 2011), which, compared to the typically low magnetic field strength in both the wind bubble and the interstellar medium, can lead to flows with high Alfvénic Mach numbers, as demonstrated by Scherer et al. (2020). Even higher velocities can be found in expanding SNRs. Considering the results obtained in this article, we also expect ions to be injected at large obliquities in SNRs like SN 1006. This object has shock speeds of the order of 3000 km s⁻¹ (see Parizot et al. 2006, and references therein) or even up to 5800 km s⁻¹ (Ressler et al. 2014). Adopting a background magnetic field of 3 µG and a gas density of 0.05 cm⁻³ (Acero et al. 2007), we find that SN 1006 shocks have Alfvénic Mach numbers ~100–200 which would allow ion injection up to obliquities ~65°. These findings raise doubts that the bipolar morphology of SN 1006 can be easily explained as the signature of a relatively homogeneous magnetic field in the environment (Bocchino et al. 2011). Finally, results obtained by Ohira (2016), that shocks propagating into a partially ionized plasma demonstrate an increased injection rate, allow for DSA in quasi-perpendicular shocks. This condition applies to many astrophysical shocks, such as a SNR expanding into cold interstellar medium.

Sato & Hughes (2017) found velocities of up to 7800 km s⁻¹ in Tycho's SNR. But the shock speed seems to be rather position dependent (Morlino & Blasi 2016, and references therein), as in its northeast part the remnant appears to interact with gas denser than typical, 0.2–1 cm⁻³ (Williams et al. 2013). Accounting for these variations and again assuming a background magnetic field strength of 3 µG, we find typical Alfvénic Mach numbers between 200 and 1000, for which obliquities above 65° may allow ion injection.

At the initial supernova break-out, the shock speed can be higher by up to an order of magnitude. Waxman & Katz (2017, page 6 and 27) predict break-out speeds of 10⁹–10¹⁰ cm s⁻¹ for

core-collapse supernovae of red and blue supergiants. Such ejecta would interact with the moderately magnetized free stellar wind, leading to shocks with high Alfvénic Mach numbers up to a few hundreds to one thousand. CR-driven instabilities are expected to grow over day timescales, leading to fast particle acceleration (Marcowith et al. 2018; Inoue et al. 2021). In these configurations, again, highly oblique shocks may allow ion injection.

Finally once turbulence is triggered by ions, electrons having sufficiently high energies can also be accelerated via DSA, giving rise to synchrotron radiation in radio and X-rays. Indeed, electrons at oblique quasi-perpendicular shocks can be accelerated up to relativistic energies (Matsumoto et al. 2017; Xu et al. 2020) injecting them into DSA via the shock internal mechanisms, such as electron shock-surfing acceleration (Matsumoto et al. 2012; Bohdan et al. 2019), magnetic reconnection (Matsumoto et al. 2015; Bohdan et al. 2020), SDA (Xu et al. 2020), and stochastic Fermi-like acceleration (Bohdan et al. 2017; Matsumoto et al. 2017). However, the exact fraction of injected electrons still should be defined with further PIC numerical experiments.

## 5. Conclusions

Combining PIC and PIC-MHD techniques, we investigated the efficiency of ions in the DSA process as a function of the Alfvénic Mach number, $M_A$, the magnetic field obliquity, $\theta_B$, and the plasma $\beta$ parameter. Our main findings are as follows:

1. The ion injection efficiency in PIC simulations of oblique nonrelativistic shocks is largely defined by the shock obliquity angle, $\theta_B$. The influence of the sonic and the Alfvénic Mach number as well as the upstream plasma beta is minor and commensurate with measurement errors. The ion energy injection efficiency can be described with the equation $U_{\text{inj}}/U_{\text{sh}} = F(\theta_B) = 10^{a_0+a_1\theta_B+a_2\theta_B^2}$, where $a_0 = -7.3$, $a_1 = 0.35$, and $a_2 = -0.0047$.

2. Shocks do not permit the triggering of DSA at a fixed obliquity if the shock Alfvénic number is below the empirical critical value $M_{cr} = F(\theta_B)^{-0.5}$, at which the energy density of injected ions in the upstream region is equal to that of the background magnetic field. For example, at oblique shocks with $\theta_B = 60^\circ$ DSA can be already triggered if $M_A \gtrsim 46$.

3. In the case of the $M_A = 50$ shock, the influence of plasma-$\beta$ on both the injection efficiency and the acceleration is weak. However, this will cease to be the case if plasma-$\beta$ is high enough that the sonic Mach number becomes smaller than $M_s \simeq 3.5$, as demonstrated by Ha et al. (2018) and van Marle (2020).

High Alfvénic Mach numbers in excess of 100, which are routinely found in the interstellar medium, can hence trigger efficient DSA up to rather high magnetic obliquity of the order of 60°. However, at higher obliquity, the critical energy criterion quickly becomes impossible to fulfill as it would require Alfvénic Mach numbers of 1000 or more, which are unlikely to occur for most astrophysical shocks. Our study does not preclude that the resonant branch of the streaming instability may control the turbulence development at low Alfvénic Mach numbers close to or even below the critical value.

PIC numerical experiments were conducted on resources provided by the North-German Supercomputing Alliance (HLRN) under the project bbp00033 and by the Prometheus
system at Academic Computer Centre Cyfronet AGH. This
work is supported by the ANR-19-CE31-0014GAMALO
project. PIC-MHD simulations were performed on the OCCI-
GEN machine at CINES under project: A0100412387

# Appendix A
## A Note on Timescales

One of the problems with comparing these simulations is that it is difficult to determine what constitutes “the same moment in time.” There are four different timescales that need to be considered. (1) The dynamic timescale of the shock, which scales with $1/V_{sh}$. (2) The response time of the upstream thermal plasma, which scales with $1/c_s$. (3) The response time of the upstream magnetic field, which scales with $1/V_A$. (4) The kinetic timescale of the nonthermal plasma, which sales with $1/v_p$, $v_p$ being the representative particle velocity. The latter, in particular, is problematic, because the average particle velocity is time and space dependent. Near the shock, the recently injected particles dominate, which have a velocity of $v_{inj}$. However, in the upstream medium, the only particles are those that have been accelerated by the SDA process. As a result, their velocity is typically about 3 $v_{inj}$. The situation becomes even more complex when comparing simulations with different plasma-$\beta_p$, which means that the timescales for gas and magnetic field in the upstream medium are no longer the same. For the sake of simplicity, we define the timescale as $\Omega_{\text{inj}}^{-1} = r_1/v_{\text{inj}}$ with $r_1$ the gyro-radius defined by $v_{\text{inj}}$ and $B_0$. However, we should keep in mind that, depending on the nature of the simulations, this does not guarantee synchronicity for all components of the plasma.

# Appendix B
## The Influence of the Box Size

As particles gain energy, they will travel further into the upstream medium before being reflected. Therefore, the size of the simulation box along the x-axis influences the SED that is obtained. This was demonstrated in van Marle (2020), which showed a much extended spectrum for a longer simulation box. However, in that case, even the small box simulation still showed that DSA took place. In our simulations, we are at the limit where DSA can take place at all and as a result, the box size can also become a critical issue. To demonstrate this, we show the results

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Influence of box length on the resulting SED.</figcaption></figure>

of two simulations: PICMHD04, and a simulation that is identical, but in a box that has only half the length along the x-axis. As shown in Figure 10, the longer box length clearly produces a more extended high-energy tail. Box size limitations can be an issue if one wants to investigate the raise of the maximum particle momentum as function of time.

## ORCID iDs

Allard Jan van Marle https://orcid.org/0000-0002-2387-4515

Artem Bohdan https://orcid.org/0000-0002-5680-0766

Paul J. Morris https://orcid.org/0000-0002-8533-8232

Martin Pohl https://orcid.org/0000-0001-7861-1707

## References

Acero, F., Ballet, J., & Decourchelle, A. 2007, A&A, 475, 883

Bell, A. R. 1978, MNRAS, 182, 147

Bell, A. R. 2004, MNRAS, 353, 550

Bell, A. R. 2005, MNRAS, 358, 181

Blandford, R. D., & Ostriker, J. P. 1978, ApJL, 221, L29

Bocchino, F., Orlando, S., Miceli, M., & Petruk, O. 2011, A&A, 531, A129

Bohdan, A., Niemiec, J., Kobzar, O., & Pohl, M. 2017, ApJ, 847, 71

Bohdan, A., Niemiec, J., Pohl, M., et al. 2019, ApJ, 878, 5

Bohdan, A., Pohl, M., Niemiec, J., et al. 2020, ApJ, 893, 6

Buneman, O. 1993, in Computer Space Plasma Physics: Simulation Techniques and Software, ed. H. Matsumoto & Y. Omura (Tokyo: Terra Scientific), 67

Caprioli, D., Pop, A.-R., & Spitkovsky, A. 2015, ApJL, 798, L28

Caprioli, D., & Spitkovsky, A. 2014a, ApJ, 783, 91

Caprioli, D., & Spitkovsky, A. 2014b, ApJ, 794, 46

Caprioli, D., & Spitkovsky, A. 2014c, ApJ, 794, 47

Dorobisz, A., Kotwica, M., Niemiec, J., et al. 2018, in Parallel Processing and Applied Mathematics, ed. R. Wyrzykowski (Cham: Springer), 156

Drury, L. O. 1983, RPPH, 46, 973

Friedman, A., Ambrosiano, J. J., Boyd, J. K., et al. 1990, Damped time advance methods for particles and EM fields UCRL-JC-106050, Lawrence Livermore National Lab

Ha, J.-H., Ryu, D., Kang, H., & van Marle, A. J. 2018, ApJ, 864, 105

Haggerty, C. C., & Caprioli, D. 2019, ApJ, 887, 165

Hanusch, A., Liseykina, T. V., Malkov, M., & Aharonian, F. 2019, ApJ, 885, 11

Inoue, T., Marcowith, A., Giacinti, G., van Marle, A. J., & Nishino, S. 2021, ApJ, 922, 7

Kato, T. N., & Takabe, H. 2010, ApJ, 721, 828

Krtička, J., & Kubát, J. 2011, A&A, 534, A97

Kumar, N., & Reville, B. 2021, ApJL, 921, L14

Marcowith, A., Dwarkadas, V. V., Renaud, M., Tatischeff, V., & Giacinti, G. 2018, MNRAS, 479, 4470

Marcowith, A., Ferrand, G., Grech, M., et al. 2020, LRCA, 6, 1

Matsumoto, Y., Amano, T., & Hoshino, M. 2012, ApJ, 755, 109

Matsumoto, Y., Amano, T., Kato, T. N., & Hoshino, M. 2015, Sci, 347, 974

Matsumoto, Y., Amano, T., Kato, T. N., & Hoshino, M. 2017, PhRvL, 119, 105101

Morlino, G., & Blasi, P. 2016, A&A, 589, A7

Niemiec, J., Pohl, M., Stroman, T., & Nishikawa, K.-I. 2008, ApJ, 684, 1174

Ohira, Y. 2013, PhRvL, 111, 245002

Ohira, Y. 2016, ApJ, 827, 36

Parizot, E., Marcowith, A., Ballet, J., & Gallant, Y. A. 2006, A&A, 453, 387

Pohl, M., Hoshino, M., & Niemiec, J. 2020, PrPNP, 111, 103751

Ressler, S. M., Katsuda, S., Reynolds, S. P., et al. 2014, ApJ, 790, 85

Reville, B., Kirk, J. G., & Duffy, P. 2006, PPCF, 48, 1741

Sato, T., & Hughes, J. P. 2017, ApJ, 840, 112

Scherer, K., Baalmann, L. R., Fichtner, H., et al. 2020, MNRAS, 493, 4172

van Marle, A. J. 2020, MNRAS, 496, 3198

van Marle, A. J., Casse, F., & Marcowith, A. 2018, MNRAS, 473, 3394

van Marle, A. J., Casse, F., & Marcowith, A. 2019a, MNRAS, 490, 1156

van Marle, A. J., Ryu, D., Kang, H., & Ha, J.-H. 2019b, PFR, 14, 4406119

Vay, J. L. 2008, PhPl, 15, 056701

Waxman, E., & Katz, B. 2017, Shock Breakout Theory, Vol. 967 (Cham: Springer)

Wieland, V., Pohl, M., Niemiec, J., Rafighi, I., & Nishikawa, K.-I. 2016, ApJ, 820, 62

Williams, B. J., Borkowski, K. J., Ghavamian, P., et al. 2013, ApJ, 770, 129

Xu, R., Spitkovsky, A., & Caprioli, D. 2020, ApJL, 897, L41