electron mass ratio $m_i/m_e = 16$. The simulation box was larger $L_x \times L_y \times L_z =$ [200 × 200 × 800]($c/\omega_{pe}$)$^3$ = [10 × 10 × 40]($c/\omega_{pi}$)$^3$ and the simulation time longer $T_{sim} = 480\omega_{pe}^{-1} = 120\omega_{pi}^{-1}$. The main result of this work is that the initial filamentation grows from electron to ion scales. However, the simulation was just long enough to reach the ionic scale, the saturation was just reached and the shock was not completely formed. Hededal and Nishikawa (2005) continue in the same direction by producing 3D simulation with similar parameters but longer time $T_{sim} = 360\omega_{pi}^{-1}$. The authors observed an interesting electron acceleration mechanism when electrons cross the ionic current channels. This produced a non-thermal tail in the electron distribution function $dN/dE \propto E^{-2.7}$. However, here again the shock was not fully formed because downstream ion distribution was still far from isotropy. The same electron energization mechanism was later found by Ardaneh et al. (2015), who studied the jet-ambient medium interaction by means of 3D simulations. These authors suggested that electrons can also be pre-accelerated by SSA mechanism during the shock formation.

Spitkovsky (2005) considered a pure $e^-$-$e^+$ plasma where there is no scale separation. In this way the typical shock formation and evolution time is much shorter than for electron-ion plasma. Different magnetizations were explored $\sigma \in [0, 0.1]$ with the relative Lorentz factor $\gamma_0 = 15$, in 3D3V configuration. The shock was triggered by reflection of incoming flow on a conducting wall. In this way, the interaction of wall-reflected and incoming flows produces the shock [simulation frame = downstream rest frame]. The box size was similar to previous studies [200 × 40 × 40]($c/\omega_{pe}$)$^3$. The shock was formed as the downstream plasma reached the expected Rankine–Hugoniot jump conditions and the distribution function isotropized in the overlap region. As previously, for the unmagnetized case $\sigma = 0$ the shock is mediated by Weibel-filamentation but strongly magnetized shock $\sigma = 0.1$ has a very different structure shaped by the perturbation of the upstream magnetic field: the incoming flow is coherently reflected on the magnetic barrier at the shock front position. The shock is then mediated by magnetic reflection. In all cases the author did not find evidence of non-thermal part in particle distribution that suggested that the acceleration is either slow to setup or not present at all. Very similar conclusions were found by Kato (2007) where 2D simulations of a shock with $\gamma_0 = 2.24$ in pair plasma were performed. This study demonstrated that the small scale magnetic field fluctuations, self consistently produced by Weibel-filamentation, is able to mediate unmagnetized collisionless shocks. The ratio of magnetic energy to the incoming kinetic energy of the upstream flow, $\xi_B = \delta B^2/(4\pi\gamma_{sh}^2\rho c^2)$, peaks at the shock front ($\xi_B = 0.14$) where the incoming mono-directional flow is isotropized and rapidly decreases downstream by phase-mixing.

**Recent multi-dimensional studies** The main difference with the early studies is (i) full formation of pair *and* electron-ion shocks and (ii) the realization of efficient particle acceleration through first-order Fermi process in a self-consistent way by following the evolution of shocks on longer timescale.

Chang et al. (2008) addressed the question of the fate of the magnetic turbulence downstream of relativistic unmagnetized shock in pair plasma. The authors used the

<figure><img src="image_36.png" /></figure>
same code as in Spitkovsky (2005) (TRISTAN-MP) with relatively long simulation time $T_{sim} = 3500 \omega_p^{-1}$, where $\omega_p$ is the total plasma frequency. As expected, the unmagnetized shock is mediated by Weibel-filamentation in the precursor. The filamentary magnetic field in the precursor becomes almost isotropic in the downstream medium once the filaments break at the shock front. The coherence scale of the field just behind the front is small $\ell_c \sim 10 c/\omega_p$ but grows with increasing distance from the shock front downstream. At the same time magnetic field intensity is found to decrease rapidly. The authors compared the simulations with an analytic model where magnetic field decreases by linear response of the plasma. The intensity is predicted to decrease as $\xi_B \propto \delta B^2/8\pi \propto (x_{front} - x)^{-q}$, with $q = 2/3$. The simulations suggested however that $q = 1$ close to the front located at $x_{front}$ and becomes closer to 2/3 far from the shock front. Due to numerical noise in PIC simulations at finite time of the simulations it is still not clear whether the field strength drops to 0 far from the shock, on macroscopic scales.

The work of Spitkovsky (2008b) presented the first self-consistent demonstration of first-order Fermi process in shocks. The same code as in Spitkovsky (2005) was used, but in 2D configuration for an unmagnetized pair plasma with upstream Lorentz factor of $\gamma_0 = 15$. The simulation time extended up to $T_{max} = 10^4 \omega_p^{-1}$ and the box size was $[10^4 \times 400](c/\omega_p)^2$. In the intermediate times the shock structure is identical to Kato (2007) and Chang et al. (2008). At late time, a supplementary population of particles builds up as a small fraction of particles in the bulk downstream plasma is able to scatter back into upstream and participate in a standard DSA (in the relativistic regime). The energy gain per cycle is consistent with the analytic prediction $\Delta E/E \simeq 1$. The non-thermal population carries typically 1% by number and 10% energy fraction of the total incoming plasma (i.e., the ratio of non-thermal electron energy to the total is $\epsilon_e \sim 0.1$).

Spitkovsky (2008a) presented the first study of electron-ion relativistic shocks where the shock is fully formed. Several mass ratios were explored $m_i/m_e = [16, 30, 100, 500, 1000]$ and the upstream plasma was unmagnetized. The most important result of this study is that electrons are brought to sub-equipartition with ions during their crossing of the precursor where they are substantially heated inside the ionic filamentary structures. In the downstream medium one gets $T_e \simeq T_i = (\gamma_0/3)m_ec^2$. This result implies an empiric similarity between shocks in pair plasma and in electron-ion plasma because the relativistic mass of particles in the downstream medium is equal. The simulations where still too short to observe the formation of a non-thermal tail in particle distributions.

Keshet et al. (2009) addressed the long term evolution of $\sigma = 0$ shocks in pair plasma by performing the longest possible simulations allowed by numerical stability. The simulation box was $[63000 \times 1024](c/\omega_{pe})^2$ and the simulation time was $T_{sim} = 12, 600 \omega_{pe}^{-1}$. They demonstrated that, as particles accelerate to larger energies with time, the precursor size increases and the width of the zone filled with magnetic turbulence increases both upstream and downstream. No convergence was reached, which leaves the question of long-term evolution open.

Martins et al. (2009), by means of 2D PIC simulations with the OSIRIS code, demonstrated that DSA works in electron-ion unmagnetized plasma. The mechanism is very
similar to the pair plasma case since electrons are at sub-equipartition with ions in the downstream medium.

The influence of magnetic field orientation with respect to the shock normal in strongly magnetized relativistic shocks ($\sigma = 0.1$ and $\gamma_0 = 15$) was studied by Sironi and Spitkovsky (2009) in pair plasma and by Sironi and Spitkovsky (2011) in electron-ion plasma, using 2D and 3D PIC simulations. These works provide a first survey of parameter space and show in which conditions relativistic shocks are efficient accelerators or not. Even if relativistic shocks are known to be generically quasi-perpendicular, the magnetic field inclination parameter, $\theta_B$, is important for the shock physics. The very special case of quasi-parallel (or subluminal) shocks, even if very rare, is interesting as it shows very different behavior. In all cases, simulations were carried out for long enough time to form the shock and see whether particle acceleration is present or not. For strongly magnetized shocks ($\sigma = 0.1$) the authors demonstrate an important difference between sub-luminal and super-luminal shocks in terms of structure and particle acceleration efficiency. In parallel shocks, the relativistic version of Bell instability is triggered and sustains an efficient DSA process. In oblique, but still sub-luminal, shocks an important contribution from the SDA mechanism was observed in competition with standard DSA. This contribution comes from the fact that the upstream plasma carries a motional electric field that can energize particles when they are reflected on the shock front. Consequently, the power-law slope of the non-thermal particle distribution function, where $dN/dE \propto E^{-\alpha}$, is not equal to the standard prediction but varies between 2.2 and 2.8. In superluminal configuration these authors did not find any particle acceleration. In this case shocks are mediated by the emission of semi-coherent electromagnetic wave from the shock front.

In order to illustrate the output from PIC simulations of shocks, in Fig. 22 we present the structure of perpendicular ($\theta_B = 90^\circ$) relativistic shock in pair plasma for mildly magnetized case $\sigma = 2 \times 10^{-3}$, obtained with the PIC code SMILEI. This structure is similar to the one found by Sironi and Spitkovsky (2009) for superluminal shocks or, more closely, to the mildly magnetized case in Sironi et al. (2013). The magnetization is chosen so that particle acceleration is efficient but maximal energy is limited by the precursor size being of the order of the Larmor radius of incoming particles $R_{L,0} = \gamma_0 m_e c^2 / (eB_0)$: $\gamma_{max} \sim 20\gamma_0$. Panels (a) and (b) present the electron density in the simulation plane and the transversely averaged profile, respectively. The shock front position is delimited by the vertical dashed line and the front propagates from the right to the left side with a velocity $v_{sh|d} \simeq 0.5c$ as measured in the downstream (simulation) frame. Ahead of the shock front oblique filamentary density structures emerge as a result of the interaction between the incoming flow and the cloud of accelerated particles. This region defines the shock precursor. Panels (c) and (d) show the longitudinal phase space $x - u_x$ and transverse phase space $x - u_y$, respectively. The transition at the shock front is clearly seen at the position where the flow becomes isotropic and hot. The cloud of energetic particles ahead of the shock front corresponds to the accelerated population. Finally, panels (e) and (f) present the particle distribution function in energy around the shock front and far downstream, respectively. Particle acceleration operates mainly around the shock front, while far downstream distribution exhibits a Maxwellian part and the start of non-thermal tail at the highest energies.

<figure><img src="image_37.png" /></figure>
<figure><img src="image_38.png" /><figcaption>**Fig. 22** Structure of a relativistic perpendicular shock in a pair plasma, at the simulation time $t\omega_p = 500$, obtained using 2D3V PIC code SMILEI. The Lorentz factor of the upstream incoming flow is $\gamma_0 = 30$ and the upstream plasma magnetization is $\sigma = 2 \times 10^{-3}$. a The electron number density in the simulation plane. b The transversely averaged electron density normalized to the upstream value. c, d The longitudinal phase space $x - u_x$ and transverse phase space $x - u_y$, respectively. e, f Present the particle distribution function in energy around the shock front and far downstream</figcaption></figure>

Haugbølle (2011) explored the differences in the structure of unmagnetized electron-ion shocks between 2D and 3D simulations. While very similar, some quantitative differences emerged in 3D simulations: the cross shock electrostatic field is slightly larger than in 2D, magnetic energy density in the shock transition region is smaller and the index of the power-law tail is closer to 2.2, instead of 2.4 in 2D. The latter is more consistent with analytical expectation (e.g., Achterberg et al. 2001).

Maximal energy of accelerated particles in perpendicular shocks was investigated by Sironi et al. (2013) by means of 2D and 3D long-term simulations. Both pair plasma and electron-ion plasma were explored for a range of magnetizations from unmagnetized case $\sigma = 0$ to strongly magnetized $\sigma = 0.1$ and for different Lorentz factors of the upstream flow ($\gamma_0 = [3, 240]$). The simulation box transverse size was $100 c/\omega_{pi}$ in electron-positron case and $25 c/\omega_{pi}$ in electron-ion case, allowing to capture at least several filaments when Weibel-filamentation mediates the shock. It was found that the
<figure><img src="image_39.png" /><figcaption>**Fig. 23** Dependence of the structure of relativistic perpendicular shocks in a pair plasma on the flow magnetization $\sigma$, at the simulation time $t\omega_p = 900$. Five representative cases are shown, from top to bottom, $\sigma = 8 \times 10^{-6}, 6 \times 10^{-5}, 4 \times 10^{-4}, 3 \times 10^{-3}$ and $2 \times 10^{-3}$. The left column shows the absolute value of the transverse magnetic field increment and the right column shows the longitudinal $x - p_x$ phase space distribution in the shock transition region. Images adapted from Plotnikov et al. (2018)</figcaption></figure>

maximum particle energy increases in time as $E_{max} \propto t^{1/2}$ for both electron-positron and electron-ion shocks. This result emerges from small-angle scattering regime of the accelerated particles in the self-excited micro-turbulence, where one expects the spatial diffusion coefficient to scale as $D \propto E^2$. The other important result of Sironi et al. (2013) study is evidencing the critical magnetization above which relativistic perpendicular shocks are not accelerating particles. For electron-positron composition the critical magnetization value is $\sigma_{crit} \approx 3 \times 10^{-3}$ and for electron-ion composition it is $\sigma_{crit} \approx 3 \times 10^{-5}$. Weakly magnetized shocks with $\sigma < \sigma_{crit}$ were found to be mediated by Weibel-filamentation that generates strong small-scale magnetic field in the vicinity of the front. In this regime DSA is efficient, with the maximum particle energy scaling as $E_{max} \propto \sigma^{-1/4}$, and a fraction of energy transmitted to the supra-thermal particles $\xi_{CR} \sim 10\%$. On the other side, for $\sigma > \sigma_{crit}$ DSA is inhibited as the shock structure is no longer dominated by the filamentation instability.

Bret et al. (2014) studied the shock formation mechanism in the unmagnetized $\sigma = 0$ case for pair plasma. An analytical model was developed, based on the growth and saturation time of the Weibel-filamentation instability. The formation time is estimated as a multiple of the instability $e$-folding time, $3\tau_{sat}$. At saturation, the density in the overlap region is 2, then the phase of density accumulation up to Rankine-Hugoniot conditions is expected to be linear in time as the incoming plasma supplies
the downstream region. The analytic model is then compared with 2D PIC simulation obtaining a reasonably good agreement.

Using 2D simulations, Plotnikov et al. (2018) provided a more systematic investigation over $\sigma$ than previously done, from unmagnetized to strongly magnetized shocks. Two different PIC codes were used, Finite Difference Time Domain (FDTD) and pseudo-spectral. Shock formation time, jump conditions, shock structure transition from low-$\sigma$ to high-$\sigma$ were investigated. The shock structure evolution for five different values of $\sigma$ is presented in Fig. 23. It shows the gradual transition from filamentation mediated shocks ($\sigma < 10^{-3}$), where particle acceleration is efficient, to magnetic reflection-shaped shocks ($\sigma > 10^{-2}$) where particle acceleration is inhibited, confirming the findings of Sironi et al. (2013). The shock formation time was found to be significantly longer than predicted by Bret et al. (2014). This points out the importance of other physical process than only saturation of Weibel-filamentation instability. For example, the studies of Vanthieghem et al. (2018), Ruyer and Fiuza (2018) demonstrate a dominant role of the drift-kink instability in the non-linear phase during which the shock front really forms. The particular focus of the study of Plotnikov et al. (2018) was on particle scattering properties, directly extracted by following self-consistent particle dynamics. The results demonstrated that the particle diffusion coefficient scales as $D = (\Delta x^2)/2\Delta t \propto E^2$ in weakly magnetized shocks, which justifies the increase in particle maximum energy of accelerated particles as $\gamma_{max} \propto \sqrt{t}$, evidenced by Sironi et al. (2013). In moderately magnetized shocks, the diffusion coefficient is modified by the presence of the ordered component that imposes a saturation of the maximum particle energy once particles get advected downstream under the effect of regular gyration.

### 4.1.2 Mildly relativistic shocks

The current consensus is that ultra-relativistic shocks are not very efficient particle accelerators for particle energies above PeV energies, mainly because of the quadratic dependence of the spatial diffusion coefficient on the particle energy, $D \propto E^2$. In the non-relativistic case, supernova remnants are considered to accelerate protons up to several hundreds of TeV or a few PeV at best and iron nuclei at energies 26 times higher. The question is then: how do the particles get accelerated to $10^{20}$ eV, maximal energy of CRs as measured at Earth? One of the promising scenarios considers mildly relativistic shocks as viable candidates (for example, trans-relativistic phases of supernova explosions or internal shocks in GRBs and jets of the AGNs). The reason is that the energy gain approaches the relativistic limit ($\Delta E/E \simeq 1$ per cycle), while a number of intrinsic limitations of the ultra-relativistic regime are alleviated, such as the generic superluminal configuration imposed by strong contraction of the pre-shock magnetic field by the shock front. Also, recent non-linear Monte-Carlo simulations demonstrated the efficiency of trans-relativistic shocks to accelerate particles to very high energies (Ellison et al. 2013).

The mildly relativistic regime is still poorly explored with kinetic simulations as of now. In the review by Marcowith et al. (2016, section 4.3) was devoted to the discussion of mildly relativistic shocks. The studies discussed there concerned mainly the plasma
physics of shock formation, but not the long term evolution. Here we provide a short
update in light of the most recent studies.

**Electron-positron plasma** Using the PIC code EPOCH Dieckmann and Bret (2017), Dieckmann and Bret (2018) investigate the generation of instabilities in a 2D configuration in the case of two interpenetrating pair plasma clouds. One of the beams is produced by the reflection at a wall of the incoming beam. The simulations focus on the generation of micro-instabilities and the shock formation process, but are not long enough to investigate non-thermal particle production. Dieckmann and Bret (2018) consider a pair plasma moving at a speed of c/2 and perform three simulations, one in 1D with a resolution of 67500c/ω<sub>pe</sub> and 3.4 × 10<sup>7</sup> macro particles, and two in 2D with the best resolution at [67500 × 1500](c/ω<sub>pe</sub>)<sup>2</sup> and using 1 billion macro particles. The two-stream and Weibel instabilities are found to rule the wave growth at the shock transition layer in this regime and to take over the filamentation instability, which nevertheless may develop upstream. We note that the Debye length has to be resolved in order to accurately capture the two-stream instability, which limits the spatial and temporal extension of the simulations.

Another interesting configuration can be found in the study of the expansion of a mildly relativistic pair plasma in a background electron–proton plasma. This setup approaches the scientific case studied in the so-called two-flow model developed to investigate gamma-ray emission in blazar jets (Sol et al. 1989). Two setups have been considered either in an unmagnetized plasma (Dieckmann et al. 2018a, b) or with a guiding magnetic field oriented along the pair plasma drift direction (Dieckmann et al. 2019). While in the former work the pair were hot with a mildly relativistic temperature of 1 MeV, the two latter works have a similar setup: the simulations are 1D with a cold pair plasma with a temperature of 400 keV moving at 0.9c, the background plasma has a realistic proton/electron mass ratio of 1836. The study follows the formation of the pair jet and the interaction between the two plasmas. It results from the free expansion of the pair beam the production of an electromagnetic piston that expels and compresses ambient electrons. The excess of negative current decelerates further the electrons but accelerates the positrons than can drift ahead the jet’s head, and reach kinetic energies of ~ MeV. In the meantime in both configurations (unmagnetized and magnetized) the pair beam and the background plasma interact through the filamentation instability which builds up a turbulent electro-magnetic field and contributes to accelerate the ambient protons also to MeV energies.

**Electron-ion plasma** Early studies of mildly relativistic shocks in electron-ion plasma (e.g., Dieckmann et al. 2008) presented important insights on the shock formation process but unfortunately their simulations were not long enough to follow-up on particle acceleration efficiency. The longest (and largest in transverse dimension) 2D PIC simulations to date were performed by Crumley et al. (2019) allowing the full formation and mid-term evolution of the shock. These authors studied the regime where the shock front velocity (in the pre-shock frame) is $\beta_{sh} \approx 0.83c$, Lorentz factor $\gamma_{sh} \approx 1.8$, and the Aflévnic Mach number of $M_A = 15$. Two different magnetic field inclinations to the shock-normal where investigated: $\theta_{Bn} = 15^\circ$ (sub-luminal) and $\theta_{Bn} = 55^\circ$ (super-luminal). The main finding is that sub-luminal (quasi-parallel) shocks are effi-
cient particle accelerators (for both electrons and ions) but not super-luminal shocks.<fnref n="42" /> When particle acceleration is efficient, the energy fraction transferred from the shock to supra-thermal ions was found to be $\epsilon_p \simeq 0.1$ (same as in non-relativistic and ultra-relativistic cases) and the energy fraction in accelerated electrons $\epsilon_e \simeq 5 \times 10^{-4}$ was found to be higher than in the non-relativistic shocks ($\epsilon_e \sim 10^{-4}$; see next subsection) but still much smaller than in ultra-relativistic shocks where electrons are in equipartition with ions, hence $\epsilon_e \sim 0.1$.

Some details of plasma physics underlying the particle acceleration efficiency were also addressed by Crumley et al. (2019). The presence of whistler waves was found in the simulation of the quasi-parallel shock, confirming the finding of Dieckmann et al. (2008), but their role in electron acceleration or injection was found to be sub-dominant, i.e., with increasing mass ratio $m_i/m_e$ from 64 to 160 whistler dynamics were expected to play more important role in electron acceleration efficiency but this effect was not observed. The maximal energy of the accelerated particles was found to increase linearly in time similarly to non-relativistic shocks, implying that the diffusion coefficient scales as $D \propto E$, resulting from efficient excitation of Bell instability in the shock precursor. On numerical side, these authors also evidenced the importance of large transverse size of the simulation box, showing that too narrow box suppresses the electron acceleration efficiency. The reason is that too narrow box suppresses Bell modes, which dominate the non-linear physics of the shock precursor.

### 4.1.3 Non-relativistic shocks

If compared to ultra-relativistic (UR) or mildly-relativistic (MR) shocks, the difficulty of capturing the full development of non-thermal tail in non-relativistic (NR) shocks comes from the fact that the energy gain per Fermi cycle (upstream → downstream → upstream) is much smaller than in UR and MR cases, as $\Delta E/E \simeq u_{sh}/c$. As a consequence, the duration of simulations must be long enough to capture at least a dozen of cycles in order to get a well-developed power-law tail while in UR shocks only a couple of cycles provides a distinguishable tail.

Despite significant efforts, only the initiation and very early stages of particle acceleration process were conveniently addressed using PIC and hybrid-PIC simulations. DSA is the most accepted model for particle acceleration at shocks. As already discussed in Sect. 2.2.3 one major difficulty for DSA to operate is the process requires particle to have a Larmor radius larger than the shock width, typically of the order of a few thermal ion Larmor radii. This concern is particularly stringent for electrons which at sub-relativistic energies have very small Larmor radii. We here discuss recent PIC simulations which address the problem of injection of electrons and ions, while the acceleration performances on dynamical timescale of the shock will be discussed in the following sections.

Table 2 presents a (non-exhaustive) list of PIC numerical experiments applied to NR shock studies. All these works are discussed in the text below.

<footnote n="42">⁴² We note that some caution has to be taken when assimilating super-luminal shocks with quasi-perpendicular shocks. In non-relativistic cases, quasi-perpendicular—but still sub-luminal—shocks can be efficient electron accelerators.</footnote>
<table><caption>**Table 2** The different PIC experiments discussed in the text referenced by author's names and publication dates</caption><thead><tr><th>References</th><th>Dim.</th><th>m<sub>i</sub>/m<sub>e</sub></th><th>θ<sup>o</sup><sub>B</sub></th><th>M<sub>A</sub></th><th>M<sub>s</sub></th><th>u<sub>sh</sub>/c</th></tr></thead><tbody><tr><td>Shimada and Hoshino (2000)</td><td>1D</td><td>20</td><td>90</td><td>3.4–10.5</td><td>1.3–4</td><td>–</td></tr><tr><td>Hoshino and Shimada (2002)</td><td>1D</td><td>20</td><td>90</td><td>32</td><td>3.3</td><td>0.375</td></tr><tr><td>Amano and Hoshino (2007)</td><td>1D</td><td>100</td><td>80</td><td>15</td><td>188</td><td>0.075</td></tr><tr><td>Amano and Hoshino (2009)</td><td>2D</td><td>25</td><td>90</td><td>14</td><td>28</td><td>0.28</td></tr><tr><td>Lembège et al. (2009)</td><td>2D</td><td>400</td><td>90</td><td>5</td><td>–</td><td>0.42</td></tr><tr><td>Riquelme and Spitkovsky (2011)</td><td>2D, 3D</td><td>40–1600</td><td>45–90</td><td>3.5–14</td><td>–</td><td>0.042–0.14</td></tr><tr><td>Matsumoto et al. (2012)</td><td>2D</td><td>100</td><td>90</td><td>14–30</td><td>28–60</td><td>0.3</td></tr><tr><td>Guo et al. (2014a)</td><td>2D, 3D</td><td>100,400</td><td>63</td><td>≤ 5</td><td>3</td><td>0.225</td></tr><tr><td>Guo et al. (2014b)</td><td>2D, 3V</td><td>100</td><td>13–80</td><td>≤ 5</td><td>3</td><td>0.225</td></tr><tr><td>Park et al. (2015)</td><td>1D</td><td>100</td><td>30</td><td>20</td><td>40</td><td>0.1</td></tr><tr><td>Kato (2015)</td><td>1D</td><td>30</td><td>30</td><td>23</td><td>–</td><td>0.37</td></tr><tr><td>Wieland et al. (2016)</td><td>2D</td><td>50</td><td>45</td><td>27.6</td><td>755</td><td>0.39</td></tr><tr><td>Bohdan et al. (2017)</td><td>2D</td><td>100</td><td>90</td><td>31–35</td><td>50, 1550</td><td>0.26–0.3</td></tr><tr><td>Matsumoto et al. (2017)</td><td>3D</td><td>64</td><td>75</td><td>21</td><td>23</td><td>–</td></tr></tbody></table>

The different columns display simulation dimensionality in space (in velocity space all simulations are '3V', tracking all three components), the ion to electron mass ratio, the magnetic field obliquity angle, the Alfvénic and sonic Mach numbers and/or the ratio of the shock velocity to the speed of light when they are available

<figure><img src="image_40.png" /></figure>
**Electron injection at non-relativistic shocks** Non-thermal electrons are responsible for the radio synchrotron emission from SNRs (Vink 2012) and are observed at interplanetary shocks (Masters et al. 2016). Several PIC simulations have explored the injection of electrons for different shock regimes and magnetic field obliquity.

At quasi-perpendicular shocks electrons can be accelerated by SSA if large amplitude electrostatic waves can develop at the shock front. These waves can be produced in the non-linear regime of the Buneman instability triggered by the relative streaming of reflected ions and incoming electrons (Wu et al. 1984; Shimada and Hoshino 2000). Hoshino and Shimada (2002) and Amano and Hoshino (2007) using 1D PIC simulations investigate the acceleration of electrons in a perpendicular fast shock. Electrons are heated by the interplay of the Buneman instability and trapped in these electrostatic waves and reflected back to the shock front by the motional upstream electric field. The maximum energy is expected to be at best $E_{max} = m_i c^2 (u_{sh}/c)$. These simulations have been generalized to multi-dimensions (2D in configuration space 3D in velocity space) by Amano and Hoshino (2009), Matsumoto et al. (2012), Dieckmann et al. (2012), Wieland et al. (2016), Bohdan et al. (2017). These works first set the criteria for SSA to occur<sup>43</sup>: (1) the thermal speed of electrons has to be smaller than the drift speed between ions and electrons, hence the shock Mach number must satisfy $M_s \ge (1+\alpha)/\sqrt{2\sqrt{m_i T_e/m_e T_i}}$ where $\alpha$ is the density ratio of reflected to incoming ions and $T_{i,e}$ are the background ion and electron temperatures, and (2) Buneman modes have to be destabilized, this requires the shock Alfvénic Mach number to satisfy $M_A \ge (1+\alpha)(m_i/m_e)^{2/3}$. These conditions depend on the ion to electron mass ratio adopted in the PIC simulations. The development of the Buneman instability and the intensity of electrostatic waves vary considerably with the number of reflected ions in the shock reformation process. Acceleration of electrons to non-thermal energies is confirmed (but see the discussion in Dieckmann et al. 2012), but the efficiency of the process depends on the background magnetic orientation with respect to the simulation plane (recall that simulations are 2D). Acceleration is the most efficient when Buneman instability-generated waves have the highest intensities, which happens when the magnetic field is out of the plane of the simulation (Bohdan et al. 2017). Electron acceleration also depends on shock non-stationarity associated to its reformation (Lembège et al. 2009). Matsumoto et al. (2017) perform 3D perpendicular shock PIC simulations of an oblique shock ($\theta_B \simeq 75^\circ$) with a high Alfvén Mach number ($M_A \simeq 21$). They find a two-step electron acceleration: first electrons gain energy via SSA in the electrostatic waves driven by the Buneman instability as above, but then they further gain energy by interacting with turbulent fields produced by the Weibel ion-ion instability triggered by the interaction of reflected and background ions. The downstream electron distribution shows the formation of a power-law energy spectrum with an index $\sim -3.5$.

Riquelme and Spitkovsky (2011) perform an extensive survey of shock conditions to investigate electron acceleration. They study the effect of variations of the shock speed, ambient medium magnetization, electron to ion mass ratio and magnetic field obliquity over non-thermal electron injection at shocks. However their simulations are restricted to rather modest Alfvénic Mach numbers $M_A < 14$. One important issue

<footnote n="43"><sup>43</sup> Bohdan et al. (2017) provide an update of these conditions for the orientation of the magnetic field with respect to the simulation plane.</footnote>
raised by the authors is that a small ion to electron mass ratio suppresses the prop-
agation of oblique whistler waves (Scholer and Matsukiyo 2004),<sup>44</sup> whereas these
waves can become over-dominant to heat/energize electrons in the foot. A criterion
for whistler wave to grow is M<sub>A</sub>/(m<sub>i</sub>/m<sub>e</sub>)<sup>1/2</sup> &lt; 1 (Matsukiyo and Scholer 2003). The
acceleration mechanism relies on the property of oblique whistler waves to have an
electric field component parallel to the magnetic field. Particles are then first acceler-
ated by this electric field before the complementary action of the convective electric
field. Electrons are preferentially accelerated at high obliquity θ<sub>B</sub> ~ 70<sup>o</sup> (at M<sub>A</sub> = 7),
where the downstream energy index is ~ 3.6. At smaller obliquities particles are not
sufficiently confined at the shock front whereas for quasi-perpendicular shocks par-
ticles can not propagate in the foot. Electron acceleration efficiency depends mostly
on the Alfvén Mach number first through the condition on whistler wave production
recalled above. The electron distribution is the hardest for Alfvénic Mach numbers in
the range 3–7. The energy index changes from 2.6 to 4 as M<sub>A</sub> changes from 3.5 to
14. A complementary study of electron acceleration in low Mach number (M<sub>a</sub> ≤ 5)
shocks was performed by Guo et al. (2014b, c) using 2D PIC simulations in order
to get better understanding of electron acceleration in galaxy cluster shocks. These
authors found that a measurable fraction of incoming upstream electron (up to 15%)
bounces back upstream and forms a non-thermal tail in the distribution function with
power-law index in energy p ≃ 2.4. These particles scatter back to the shock front
on self-generated waves via firehose instability and participate in the SDA process.
This acceleration process was found to be efficient if upstream plasma is high beta
(β ≥ 20) for nearly any magnetic field obliquity.

High Alfvénic Mach number, quasi-parallel shocks could also allow electron injec-
tion. These shocks are likely good proton injectors (see below). In turn protons (ions)
can trigger magnetic perturbations, as the magnetic field grows in the pre-shock
medium then lowering M_A and its transverse component can be compressed at the
front. The conditions then resemble the case of highly-oblique moderate Alfvén Mach
number shocks discussed by Riquelme and Spitkovsky (2011) [see also Caprioli and
Spitkovsky (2014b)]. Park et al. (2015) perform long term 1D PIC simulations of
high M_A quasi-parallel shocks (see also Kato 2015). Protons destabilize non-resonant
streaming (Bell) modes and electrons are accelerated by a combination of SDA and
Fermi processes as they are scattered by the non-resonant streaming modes. Interest-
ingly, non-thermal electrons entering in the relativistic regime show a E⁻² energy spec-
trum and a non-thermal electron to proton ratio ~ 10⁻³ roughly proportional to uₛₜ/c.

**Ion injection at non-relativistic shocks**
The major drawback of full-PIC simulations to address the ion injection into DSA is the need to resolve both electron- and ion-scale physics, that bakes typical simulation not longer than a ~ 10ω<sub>ci</sub><sup>-1</sup>. This difficulty is partly bypassed using hybrid-PIC simulations where ions are still treated kinetically but electrons are treated as massless fluid (see,e.g., Lipatov 2002; Gargaté et al. 2007; Kunz et al. 2014). In this approach all ion-scale kinetic physics are preserved while the global numerical cost is about two orders of magnitude lower than in full-PIC simulations.

<footnote n="44"><sup>44</sup> Whistler waves are likely excited due to the cross-field drift of background electrons with respect to either reflected or background ions, this is the so-called modified two-stream instability or MTSI.</footnote>
The injection of thermal ions into the acceleration process was investigated by Guo and Giacalone (2013) and Caprioli et al. (2015) using multi-dimensional hybrid-PIC simulations. The first main finding of these studies is that protons are not injected by 'thermal leakage' of downstream thermalized distribution into pre-shock medium but by specular reflection on time-varying shock barrier. For quasi-parallel shocks ($\theta_{Bn} \le 45^\circ$) and high $M_a > 5$, the injection efficiency is larger than 10%. As evidenced by Caprioli et al. (2015), protons gain energy through SDA in consecutive reflections on the shock front and inject into DSA when their energy is large enough to escape upstream. They also propose a quantitative model that accounts for the drop in injection efficiency of quasi-perpendicular shocks with $\theta_{Bn} \ge 45^\circ$, as more than 4 SDA cycles are required for injection into DSA, while at each SDA cycle a large fraction of ions ($\sim 75\%$) is lost downstream. This effect explains the rapid drop in injection efficiency of $\theta_{Bn} \ge 45^\circ$ shocks (see, however, Ohira 2016).

Other studies addressed the thermalization of heavy ions in post-shock medium and chemical enhancement in shock accelerated particles. It was found that each species acquires downstream temperature proportional to its mass, $T_d \propto A_i$, where $A_i$ is the atomic number (Kropotina et al. 2016; Caprioli et al. 2017). The efficiency of ion injection into DSA increases with $A/Z$ ratio, where $Z$ is the charge. Caprioli et al. (2017) show that there is preferential acceleration of ions with large $A/Z$ in quasi-parallel shocks. For $M_a > 10$ these authors find that the fraction of DSA-accelerated ions scale as $(A/Z)^2$, in quantitative agreement with abundance ratios in Galactic Cosmic Rays. The injection mechanism of heavy ions is different from proton injection, since they do not have any dynamical impact on the shock structure. Instead of reflecting specularly on the shock front, heavy ions directly thermalize in the post-shock medium. If the downstream isotropization time is shorter than advection time, a small fraction can back stream into pre-shock medium and participate in DSA. We note that Hanusch et al. (2019), using 2D hybrid-PIC simulations as well, confirm the preferential injection of heavy ions up to $A/Z \sim 10$ but find that there is saturation for higher $A/Z$ values, in contrast with Caprioli et al. (2017) findings. For quasi-perpendicular shocks with $\theta_{Bn} > 50^\circ$, similarly to pure electron-proton shocks (Caprioli and Spitkovsky 2014a), there is no injection into DSA of any ion species, since advection time becomes shorter than isotropization time.

### 4.1.4 Discussion

Let us summarise and briefly discuss the micro-physical studies of particle acceleration at collisionless shocks. Several key points emerge in recent studies:

- The self-consistent shock structure produces non-thermal particle distributions for a broad range of parameter space. The main requirement for efficient particle acceleration is the ability of particle impinging the front to escape back into the pre-shock medium and trigger wave growth through kinetic instabilities. This condition is generally met in high-Mach number sub-luminal shocks (or weakly magnetized shocks).

- Unmagnetized relativistic shocks are efficient particle accelerators in both electron-positron and electron-ion plasma. Here, it was demonstrated that the
Weibel-filamentation instability mediates the shock transition region. It leads to significant magnetic field amplification and concomitant particle acceleration. Typically, about 1% by particle number fraction and 10% by the shock kinetic energy fraction is channelled to non-thermal particles. However, the acceleration rate is slow as it scales quadratically with particle energy: $t_{acc} \propto E^2$. Thus, one infers the maximum energy of protons achievable in ultra-relativistic shocks of GRBs as $E_{max} \sim 10^{16}$ eV.

- When the relativistic shock front propagates into highly magnetised medium ($\sigma > 0.01$) one has to distinguish between sub-luminal and superluminal configurations. In the former case, particle acceleration was found to be efficient, sometimes even more than in unmagnetised case. Yet, the sub-luminal configuration is statistically disfavoured in relativistic case. When the configuration is superluminal, shock-processed particles are advected downstream and are unable to undergo Fermi process.

- At intermediate magnetization, $\sigma_{crit} < \sigma < 0.01$, limited particle acceleration occurs and the maximum energy scales as $E_{max} \propto \sigma^{-1/4}$.

- The interesting case of mildly-relativistic shocks, e.g., $\gamma_{sh} \ge 1$, where energy gain per Fermi cycle is large and the shock can easily be subluminal is poorly studied. Recent study by Crumley et al. (2019) found that shock physics in quasi-parallel case is similar to non-relativistic shocks.

- Non-relativistic shocks are the most common and studied in literature. The most representative case is the external shock of the Supernova Remnants. In this regime, when particles are efficiently reflected on the front, several instabilities can be in competition in the precursor region (e.g., Buneman, firehose, Whistler, Weibel, gyroresonant, Bell), depending on the Mach number, magnetic field strength and obliquity. Hence, the phenomenology is more complex then in ultra-relativistic shocks. For example, quasi-parallel shocks are common, contrary to the ultra-relativistic regime, and lead to efficient magnetic field amplification through resonant and Bell instability in the shock precursor.

There are several open questions under active investigation or to be addressed in the near future.

1. Are ultra-relativistic shocks always locked in the slow acceleration rate, i.e., $t_{acc} \propto E^2$, or an additional source of magnetic turbulence can produce faster acceleration? How is the long term evolution of unmagnetized shocks where the shock transition is governed by self-excited microturbulence? No steady state was reached with current simulations. For recent progress in this field, see Lemoine et al. (2019a).

2. The question of which configuration, quasi-parallel or quasi-perpendicular, in non-relativistic shocks is more efficient for ion/electron acceleration? The regime of quasi-perpendicular but still sub-luminal shocks is a particular case that requires clarification.

3. How promising is the mildly relativistic regime?

4. On numerical side, important efforts are undertaken to push PIC simulations to the largest scales and longest time benefiting from modern computational resources. While largely needed, gaining one order of magnitude in system size and simulation time becomes rapidly prohibitive even with the largest available supercomputers.
In this respect, hybrid approaches such as MHD-PIC are promising when one
is mainly interested in dynamics of supra-thermal particles (see Sect. 5.5). Yet,
this requires to robustly prescribe how some part of thermal particles (simulated
using the fluid approach) are promoted to non-thermal status (simulated using PIC
or Vlasov approach). For example, in the shock problem one prescribes a fixed
fraction of shock-processed particles to be injected into the non-thermal pool (Bai
et al. 2015; van Marle et al. 2018) but more accurate parametrisation is required
when the shock structure becomes modified by non-thermal particles.

In conclusion, PIC simulations provide detailed non-linear solutions of the shock problem. Therefore, they are an efficient tool to probe the efficiency of particle acceleration for a given parameter set. They are, however, of limited duration (a couple of ion gyro-periods for full-PIC and ~ 10³ ion gyro-periods for hybrid-PIC) and box sizes are typically less than thousands of ion skin depths. Even with the most powerful current computational facilities, simulating the global astrophysical systems is unachievable with this approach. The goal and common approach is to provide robust scalings which can be included in fluid simulations as sub-grid prescriptions.

## 4.2 Kinetics of magnetic reconnection

There is a vast literature on magnetic reconnection, both for the collisional case based on resistive MHD and the non-collisional case based on the Maxwell-Vlasov equations. Nevertheless, many questions are still far from understood, including ‘what triggers reconnection events in real astrophysical objects’, ‘what are the physical processes which accelerate particles to super-high Lorentz-factors’, ‘do associated energy spectra always show power-law slope and what is the spectral index’, ‘can ions be accelerated to equal energies as electrons’, ‘is there an upper limit for the Lorentz-factors that can be achieved and which process sets this upper limit’? A deeper understanding of these questions will definitely help to answer relevant astrophysical questions: (1) To what degree is magnetic reconnection important for the dynamics of large scale flows like the launching of jets from compact objects or driven shock waves? (2) To what degree is REC responsible for the production of thermal and non-thermal high-energy photons observed from the Sun to AGNs? (3) To what degree can REC accelerate ions to relativistic speeds and can thus contribute to the cosmic ray flux and the hadronic channel of emissivity of photons and energetic neutrinos?

As the literature is vast there is no chance to refer to all papers. In a hopefully not
too biased view, basic ideas are thus presented on the basis of selected papers in 4
subsections. (1) What kinetic simulations can achieve and why we decisively need
them, (2) Some key results based on kinetic simulations, (3) The most prominently
discussed physical processes able to accelerate particles, (4) A critical discussion and
outlook. For further important points, which are not discussed due to lack of space,
we refer to the reviews given at the beginning of Sect. 2.5.

4.2.1 What kinetic simulations can achieve and why we decisively need them

Microphysical studies have the great advantage that they rely only on fundamental physics. Difficult questions—like which equations of the MHD family best
model reconnection events and which values of the transport coefficients (resistivity, viscosity, Hall-parameters, or even higher moments) are most appropriate—can be omitted. All this comes self-consistently from the kinetic physics solved, described by Vlasov-Maxwell equations given in Sect. 3.1. The numerical method most widely used in astrophysics is the PIC method described in Sect. 3.3, though solvers for the Vlasov equations in the 6+1 dimensional phase space start to appear.

The price to pay is that the computational costs to solve kinetic equations are much higher than those to solve any MHD model—even complex ones like the popular 10-moment closure model, see for instance Wang et al. (2015) or Lautenbach and Grauer (2018). With kinetic simulations, even huge ones, only local aspects can be addressed, on spatial scales which measure at most some thousand electron inertial lengths and thus up to about hundred proton inertial lengths. This is sufficient to study a single X-point or even a small current sheet breaking up into a plasmoid chain. It is, however, not sufficient to answer questions like how these reconnection sites are embedded in the large scale environment and how they have been formed. One way out, which is included in the discussion below, is to apply some splitting between MHD models and localized particle models or to post-process MHD solutions by propagation of test-particles. In the future, more numerical codes will probably be available which dynamically couple MHD and particle dynamics. Some references to such codes are given at the end of Sect. 2.5.

Another limitation of kinetic simulations is that they do not include photons and particle physics. Such aspects are ultimately relevant for high-energy plasma processes. Power law slopes and cutoffs may change when considering energy losses by the emission of synchrotron radiation or inverse Compton scattering of electrons on colder photons originating from the reconnection site itself or from other, external processes, possibly far away from the reconnection site. In relativistic REC, where the energy involved exceeds the rest masses of electrons (and maybe even of protons), the building of electron-positron-pairs is very likely to take place, which significantly back-react on the reconnection dynamics. The same may be true if accelerated protons can create pions and higher hadronic resonances and neutrinos. Attempts to account for radiative losses have recently been made, but are only at the very beginning.

Despite these limitations, kinetic simulations have brought, in about the last 10 to 15 years, an immense progress in our understanding of both, magnetic reconnection and shock waves, as well as of associated particle acceleration. This progress could not have been achieved on the basis of pure MHD simulations.

**Simulation setup** Most kinetic simulations of magnetic reconnection have been performed in two spatial dimensions or in simple prolongations into the third dimension (but see the paragraph on dimensionality) of the setup sketched below. In two space-dimensions, there are known two analytical, stationary, though unstable configuration of current sheets for the Vlasov-Maxwell equations. These solutions are typically used to setup a kinetic simulation to study magnetic reconnection.

(1) *Harris-sheet*: the reconnecting magnetic field $B_0$ is oppositely oriented along one axis (say the x-axis in Fig. 24). Normal to it (along the y-axis), the field strength varies as $\mathbf{B}_{x,\text{rec}}(y) = \hat{x}B_0 \tanh(y/\delta)$. The thickness of the sheet, $\delta$, typically is a few electron inertial lengths and thus, for realistic mass ratios between
<figure><img src="image_41.png" /><figcaption>**Fig. 24** Acceleration sites in a Harris-sheet broken into a chain of plasmons and X-points (from an unpublished PIC simulation by D. Folini and R. Walder). Shown are accelerated electrons on top of the electron number density (top panel), magnetic field orientation (middle panel), and magnetic field strength (bottom panel). Possible acceleration sites are indicated by different labels: (A) Reconnecting electric field. (B) First order Fermi-mechanism between the converging inflows to the reconnection site. (C) Contracting plasmons. (D) Merging plasmons. (E) Drifts in inhomogeneous magnetic fields. (F) Turbulence</figcaption></figure>

electrons and ions, less than a proton inertial length. The plasma density varies as $n(y) = n(0)/\cosh^2(y/\delta)$ along the y-axis. Together with an appropriate temperature, the thermal pressure within the sheet balances the magnetic pressure from outside the sheet. The induced current then points into the normal direction and writes $J(x) = cB_0/4\pi\delta \sech^2(x/\delta)\hat{z}$.

(2) *Force free equilibrium*: a basically similar equilibrium can be obtained using the force-free equations (see, e.g., Guo et al. (2015); Wilson et al. (2017) for details).

The parameters of both configurations are symmetric with respect to the center of the sheet, except that the direction of the magnetic field is reversed. Some authors use for their setup many such current sheets aligned in parallel, (see Biskamp 2000; Drake et al. 2010; Kowal et al. 2011; Kagan et al. 2013; Werner et al. 2018). This can be quite natural, for instance within the frame of striped wind models of the solar corona (Li et al. 2014) or at the termination shock of the solar wind at the heliopause (Drake et al. 2010) or in neutron star nebulae (Pétri and Lyubarsky 2007; Kirk et al. 2009; Sironi and Spitkovsky 2012; Uzdensky and Spitkovsky 2014; Cerutti et al. 2016; Cerutti and Philippov 2017).

An extension, and a step towards more realistic models, is to use asymmetric initial configurations (Cassak and Shay 2007; Belmont et al. 2012; Hesse et al. 2013; Aunai et al. 2013; Pritchett 2013; Eastwood et al. 2013). Within this review we do not further
elaborate on this situation. We concentrate on the 2D Harris-sheet instead, what can be learned from it and where the limitations of this toy setup lie.

If the thickness of such current sheets is less than a typical field diffusion length, they are linearly unstable to tearing modes. Eventually, X-points will develop where REC will start. In simulations, one often uses overpressured inflows to accelerate the development of the instability. As described below, each X-point will develop its exhausts. Exhausts of different X-points may collide to form magnetic islands (also called plasmoids or O-points), see Fig. 24. Waves are generated and with them turbulence develops within the sheet and outside of it, on both its sides in the diffusion regions. These ingredients turn out to be important drivers of non-thermal particle acceleration.

Some authors trigger one single X-point in their simulations, from which huge exhausts develop. Within the exhausts, secondary X-points and islands can develop. Many features described below are valid for both, the triggered and the un-triggered approach. Both approaches develop a power-law spectrum of accelerated non-thermal particles (Figs. 25, 26, 27). Differences between the two approaches and whether one approach is more close to a reconnection event in astrophysics are not yet worked out in detail (but see some points discussed below).

Different boundaries can be used and these will be discussed when presenting simulations and associated results below.

**The parameters governing the physics** It turns out that the evolution of the sheet and the acceleration of particles depend on the strength of the reconnecting magnetic field, the particle density, and the field temperatures of electrons and ions. Often, the magnetization is given as the ratio $\sigma_s$ of magnetic energy density to enthalpy density, where s stands for the particle species, either electrons or ions, s = e, i. To account for temperature effects, one may define $\beta_s$ as the ratio of particle thermal pressure and magnetic pressure. Other, equivalent, characterizations used are $\sigma_{\text{hot},s}$, the ratio of the energy flux in the reconnecting magnetic field to that in the particles for each species (thermal and bulk, and particle rest mass) and $\sigma_{\text{cold},s}$, which does not consider temperature and bulk flow.

A particular parameter is the background field. The term refers to the field component which is not going to reconnect, i.e., a component which encloses a certain angle, $\theta_B$, with the inverted field components that drive the reconnection process. If a background field is present, the reconnection physics is going to alter drastically.

As discussed in Sect. 2.5, the physics of collisionless REC is different for pair plasmas on the one hand and for an electron-ion plasma on the other hand. The much larger masses of ions as compared to those of electrons results in their much earlier de-magnetization when the plasma flows into the reconnection region (consult Fig. 10), with a series of consequences (Melzani et al. 2014b). Relativistic settings also result in peculiarities (Melzani et al. 2014a, b). In the ultrarelativistic limit (i.e., when for both, electrons and ions, the energy largely exceeds the rest-mass), plasma time- and lengths-scales (cyclotron frequencies, skin depths, Larmor radii) become independent of the particle rest mass and depend only on the particle energy. Guo et al. (2016) showed that, in this ultra-relativistic limit, a pair and ion/electron plasma behave essentially similarly.

<figure><img src="image_42.png" /></figure>
<figure><img src="image_43.png" /><figcaption>**Fig. 25** Upper row: final electron spectrum for various simulations, without (left panel) and with a guide field (right panel). Lower row, left panel: Lorentz factor distributions of electrons of a simulation by Melzani et al. (2014b), using $\omega_{ce}/\omega_{pe} = 3$ and a magnetization in the background plasma $\sigma_{hot,i,e} = 3.6, 83$. Red curves indicate particles originally found in the setup current sheet, green curves indicate particles from the plasma which flows into the sheet during the reconnection process. For the green curves, times are ordered as dark to light green, with values 0, 750, 1500, 2250, 3000, 3750 $\omega_{ce}^{-1}$, i.e., one curve every 750 $\omega_{ce}^{-1} = 250 \omega_{pe}^{-1} = 50 \omega_{pi}^{-1} = 30 \omega_{ci}^{-1}$. The blue dashed line indicates the final power-law slope of the background accelerated particles. Lower right panel: the time-evolution of the maximum Lorentz factor of the background particles for various simulations with $m_i/m_e = 25$ or 1. Solid lines are for electrons, dashed lines for ions and represent $m_i/m_e \gamma_{1,max} \cdot \tilde{\omega}_{c,e} = \omega_{c,e}/\omega_{p,e}$. The index s refers to the power law index $t^s$. Images adapted from Melzani et al. (2014b)</figcaption></figure>

Regarding PIC simulations themselves, care must be taken that the number of super-particles per discretization cell is sufficiently high to ensure that collisionless kinetic processes remain faster than collisional effects, e.g., for thermalization. For this, the PIC-plasma parameter, $\Lambda^{PIC}$, the number of superparticles per Debye sphere, has to be sufficiently large, see Melzani et al. (2013) for a more thorough discussion. In addition, when using an explicit PIC scheme on a Yee grid (in particular of low order), care must be taken that energy conservation is sufficiently well guaranteed and that artificial Cherenkov radiation (Greenwood et al. 2004) does not dominate the scene. The problem appears in particular for particles traveling close to the speed of light, where the high-frequency waves from artificial Cherenkov radiation actively interact with the particles and thus may influence the acceleration process in a non-physical way (see section 3.3.3).
<figure><img src="image_44.png" /><figcaption>**Fig. 26** Upper row: final electron (blue-solid) and proton (green-dotted) spectrum of non-thermal particles for a low ($\sigma_i = 0.03$, left panel) and a high magnetization ($\sigma_i = 30$, right panel). Lower left panel: fit for the electron power law index $p$ (Eq. 46). Lower right panel: final energy partition between background electrons and ions at different ionizations $\sigma_i$. As expected, there is no difference between electrons and ions in the ultra-relativistic regime. The fit relates to Eq. (45). Images adapted from Werner et al. (2018)</figcaption></figure>

### 4.2.2 Some key results based on kinetic simulations

In this section, selected results are presented regarding the energy spectra of accelerated particles, the energy partitioning between ions and electrons, and the overall efficiency of the energy transfer from the reconnecting magnetic field to the accelerated particles. These results set the stage for a more detailed view on the physical mechanisms behind particle acceleration, to be presented in the next section.

We concentrate on the case of a relativistic electron-ion plasma in an idealized 2D Harris-sheet set up. The non-relativistic case is important to understand processes on the Sun, space-weather, but also technical devices, from plasma thrusters to tokomaks. The case of relativistic pairs is interesting in pulsar winds and perhaps in certain regions of a black hole corona—in every environment where pair-cascades could develop. However, these two cases will not be discussed due to lack of space and because many results have been presented in other reviews (references given at beginning of Sect. 2.5).

We call REC semi-relativistic when the energy of the magnetic field exceeds the electron rest mass energy but is smaller than the ion rest mass energy, otherwise we

<figure><img src="image_45.png" /></figure>
<figure><img src="image_46.png" /><figcaption>**Fig. 27** Upper left panel: the $\beta_i$-dependence of the final ($t/t_A = 2$) electron spectra at low magnetization. Upper right panel: time-evolution of the electron spectrum for $\sigma_i = 1$ and $\beta_i = 0.16$. The spectrum starts to develop an additional component after about one $t_A$. Lower panels: electron spectra for different initial and boundary conditions. For details, see text. Image adapted from Ball et al. (2018)</figcaption></figure>

call it relativistic. Ultra-relativistic REC terms the situation where the magnetic energy largely exceeds the rest mass energy of both species, electrons and ions.

There are only a few papers that have addressed the relativistic electron-ion regime though it is decisive for our understanding astrophysical high-energy objects: the dynamics and emission of accretion disks around and jets from either black holes or neutron stars (see Table 3).

Melzani et al. (2014a) and Melzani et al. (2014b) present a first study of collisionless relativistic magnetic reconnection of an electron–proton–plasma in two space dimensions. They use reduced mass ratios $m_p/m_e = 25$ and $50$ respectively, and started from a Harris-equilibrium. Periodic boundary conditions normal to the current sheet and reflecting conditions parallel to the sheet are used. Tearing instability and subsequent REC develop from numerical noise. The magnetization varies as $10 \le \sigma_e \le 260$ and as $0.4 \le \sigma_i \le 14$. Except for one simulation $2.8 \cdot 10^{-5} \le \beta_{p,i} = \beta_{p,e} \le 2 \cdot 10^{-3}$. They present cases without a guide field and with a guide field of $B_G = 0.5B_0$ and $B_G = B_0$, where $B_0$ is the reconnecting component of the field, the component (anti-)parallel to the current sheet.

Melzani et al. (2014a) describe the general structure of the reconnection process for both, the cases with and without guide field. The bulk inertia is identified as the main
<table><caption>**Table 3** Order of magnitude for physical parameters in astrophysical environments. Adapted from Melzani et al. (2014b)</caption><thead><tr><th>Objects with ion-electron plasmas (with also pairs)</th><th>B (G)</th><th>n<sub>e</sub> (cm<sup>-3</sup>)</th><th>T<sub>e</sub> (K)</th><th>σ<sub>e</sub><sup>cold</sup></th><th>U<sub>A,in</sub>/c</th></tr></thead><tbody><tr><td>Microquasar coronae, X-ray emitting region<sup>a</sup></td><td>10<sup>5</sup>–10<sup>7</sup></td><td>10<sup>13</sup>–10<sup>16</sup></td><td>10<sup>9</sup></td><td>10<sup>-1</sup>–10<sup>5</sup></td><td>0.003–1</td></tr><tr><td>AGN coronae, X-ray emitting region<sup>b</sup></td><td></td><td></td><td>10<sup>9</sup></td><td>1.7–180</td><td>0.03–0.3</td></tr><tr><td>Giant radio galaxy lobes<sup>c</sup></td><td>10<sup>-6</sup>–10<sup>-5</sup></td><td>3 × 10<sup>-6</sup></td><td>10<sup>6</sup></td><td>0.8–80</td><td>0.02–0.2</td></tr><tr><td>Extragalactic jet, γ-ray emitting region (< 0.05 pc)<sup>d</sup></td><td>12</td><td>80</td><td></td><td>2 × 10<sup>5</sup></td><td>~ 1</td></tr><tr><td>Extragalactic jet, radio emitting region (kpc scales)<sup>e</sup></td><td>1–3 × 10<sup>-5</sup></td><td>0.8–5 × 10<sup>-8</sup></td><td></td><td>500–2500</td><td>~ 1</td></tr><tr><td>GRB jet, at radius of fast reconnection<sup>f</sup></td><td>7 × 10<sup>8</sup></td><td>10<sup>10</sup></td><td>10<sup>8</sup></td><td>5 × 10<sup>12</sup></td><td>0.9</td></tr><tr><td>Objects with pair plasmas</td><td>B (G)</td><td>n<sub>e</sub> (cm<sup>-3</sup>)</td><td>T<sub>e</sub> (K)</td><td>σ<sub>e</sub><sup>cold</sup></td><td>V<sub>A,in</sub>/c</td></tr><tr><td>At the termination shock of pulsar winds<sup>g</sup></td><td>10<sup>4</sup></td><td>0.1–10</td><td></td><td>10<sup>13</sup></td><td>~ 1</td></tr><tr><td>In a pulsar wind nebulae<sup>h</sup></td><td>5 × 10<sup>-3</sup></td><td>5 to 10<sup>3</sup></td><td>γ ~ 10<sup>-9</sup></td><td>&lt; 0.5</td><td>0.6</td></tr></tbody></table>

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

<footnote><sup>a</sup> Analytical disk and corona models, de Gouveia dal Pino and Lazarian (2005), Di Matteo (1998), Merloni and Fabian (2001), Reis and Miller (2013); matching observed spectra with radiation models Del Santo et al. (2013), Romero et al. (2014)</footnote>

<footnote><sup>b</sup> Analytical disk and corona models, Di Matteo (1998), Merloni and Fabian (2001), Reis and Miller (2013)</footnote>

<footnote><sup>c</sup> Observations Kronberg et al. (2004)</footnote>

<footnote><sup>d</sup> Analytical model assuming σ<sub>cold</sub> = 100, Giannios et al. (2009), See also Giroletti et al. (2004) for magnetic field measurements (0.2G, but on larger scales)</footnote>

<footnote><sup>e</sup> Observations, Schwartz et al. (2006). See also Romanova and Lovelace (1992)</footnote>

<footnote><sup>f</sup> Analytical model, McKinney and Uzdensky (2012). Pairs are also present, with n<sub>pair</sub> ~ 10n<sub>e</sub></footnote>

<footnote><sup>g</sup> Analytical model and observations, Bucciantini et al. (2011), Sironi and Spitkovsky (2011)</footnote>

<footnote><sup>h</sup> Analytical model and observations, Bucciantini et al. (2011), Sironi and Spitkovsky (2011). The plasma distribution function is a broken power-law with Lorentz factors γ in the indicated range. We note that Cerutti et al. (2013) considers only the high-energy electron population, and hence has larger magnetizations</footnote>
non-ideal process, which de-magnetizes both, electrons and ions. In the energy flux of the exhausts the thermal component dominates over the bulk component. Protons are generally hotter than electrons in the exhausts. The numerical results correspond to analytical estimates given. They identified a good measure for the relativistic reconnection rate,

$$E^* = E_{\text{Rec}}/B_0 U_{\text{in}}^{\text{R}}, \qquad (44)$$

with $U_{\text{in}}^{\text{R}}$ the relativistic Alfvén speed, $B_0$ the reconnecting field, and $E_{\text{Rec}}$ the reconnection electric field. This rate varies between 0.14 and 0.25, showing higher values with lower background densities. Generally, the rate is higher than in the non-relativistic case (for which $0.07 \le E^* \le 0.14$), which is in line with the findings for pair plasmas, e.g., by Zenitani and Hoshino (2007) ($E^* = 0.2$), by Cerutti et al. (2012b) ($E^* = 0.17$), or by Bessho and Bhattacharjee (2012) ($E^* = 0.19$ and $0.36$). It follows directly that the reconnection electric field is very large, $E_{\text{Rec}}/B_0 \sim 0.2U_{\text{A,in}}^{\text{R}} \sim 0.2c$, which turns out to be important for the acceleration process of particles.

Melzani et al. (2014b) further elaborate on associated results, notably on the acceleration of non-thermal particles to very high Lorentz-factors. The results are summarized in Fig. 25.

REC produces a population of non-thermal particles. This population shows a power-law spectrum, which gets harder for increasing magnetization (upper left panel). To first order, the slope of the power law is independent of the presence of a background field (upper right panel). There is a clear difference in the spectrum of accelerated particles (lower left panel), depending on whether the particles have been present in the initial current sheet (population CS) or whether particles have floated into the sheet after REC has started (population BG). Population BG forms a power-law, population CS is substantially heated but its spectrum is still Maxwellian. This points to a different acceleration mechanism of the two populations, see next section.

The lower right panel of Fig. 25 shows the temporal evolution of the maximum Lorentz-factor $\gamma_{\text{max}}(t) \sim t^s$, which is close to an exponential growth. Higher magnetization $\sigma$ in a larger coefficient $s$ and thus a faster growth of the spectrum to high Lorentz-factors. While, as said above, non-thermal accelerated particles have the same spectrum, independently whether a background field is present or not, the growth time of the power-law cutoff $\gamma_{\text{max}}$ is clearly slower in the presence of a background field. Also, to first order, $\gamma_{\text{max}}$ for protons grow slower than electrons by a factor $m_i/m_e$.

Finally, while in the thermal exhausts there is much more energy in the ions than in the electrons, the situation is slightly different for the accelerated particles. Without a background field, there is more energy in the ions. But the situation reverses when a background field is present, when more energy is carried by the electrons. We add a note of caution to this result as the mass-ratio between electrons and ions used is either 25 or 50. Simulations with a realistic mass ratio should be undertaken and their results be confronted with the findings of Melzani et al. (2014b).

Guo et al. (2016) perform PIC simulations of relativistic electron-ion reconnection (magnetically dominated in their terms) without a guide field, starting from a force-free equilibrium of the current sheet. They use different mass-ratios between electrons and ions, between 1 and 1836 and use different domain sizes and different inflow temperatures. They explore magnetizations reaching from the relativistic to the
ultra-relativistic case, but cover not really the semi-relativistic case, in contrast to the other work discussed in this section. This should be kept in mind when looking at the following results. They find that for low mass ratios, ions gain slightly (1.1 times) more energy than electrons while for a real mass ratio, the ions gain 1.5 times more energy than the electrons. All power-law slopes are hard – when put into the form $f \propto (\gamma - 1)^{-s}$, s is between 1 and 2 and very close to 1 for high magnetizations. For their high magnetizations, the electron power-law slope does not depend on the mass ratio and is the same for a pure pair plasma, for a mass ratio of 100 and for a realistic mass ratio. In energy space, the slopes for electrons and ions slightly differ, e.g., 1.35 for electrons and 1.2 for ions for a simulation with $\sigma_0 = B_0^2/(4\pi(m_i+m_e)nc^2 = 100$ ($n$ is the particle number density) and $m_i/m_e = 100$. However, the momentum distribution shows $s_p = 1.35$ for both species. They argue that this can only be achieved by a Fermi-like acceleration mechanism. They find a slight dependence of the power laws on the size of the current sheet, not so much for the all over slope, but for secondary variations of the slope. Smaller domain sizes show more variations, with a significant change of the slope for different energies ($\gamma_e - 1$). Larger domain sizes show a much more smooth, unique slope over the entire energy range. Finally, these authors emphasize that all their power-laws show an exponential cutoff.

Werner et al. (2018) present an extensive 2D study of collisionless REC in the semi-relativistic and relativistic regime, varying the ion magnetization, $\sigma_i$, from $10^{-3}$ to $10^4$. The plasma-beta value of ions was 0.01 for all simulations. A realistic mass ratio between ions and electrons is used. For their setup, they use a doubled Harris sheet with two field reversals and they use period boundary conditions in both directions. The authors analyze the plasma flows in the 'thermal' regime and describe the Hall signatures due to the different sizes of the electron and ion diffusion region (see Fig. 10). They find a reconnection rate of about 0.1 of the Alfvénic rate (their slightly different definition of the reconnection rate) across all regimes, slightly below (0.08) in the semi-relativistic regime and slightly above (0.12) in the ultra-relativistic regime. In the ultra-relativistic limit, the release of magnetic energy during REC is distributed equally between electrons and protons, but protons gain more in the semi-relativistic regime, up to 75 % for the weakest $\sigma_i$ (see Fig. 26, lower right panel). The authors present a formula for the fractional energy gain of the electrons, $q_e$,

$$q_e = \frac{1}{4} \left( 1 + \sqrt{\frac{\sigma_i/5}{1 + \sigma_i/5}} \right). \qquad (45)$$

Integration is done till 2000 Larmor time-scales. Particles are accelerated to a non-thermal regime. Power-laws and energy-cut off seem to saturate till the end of the simulation. By fitting their data to a power-law, they find for electrons the time-saturated relation

$$f(\epsilon) \sim \epsilon^{-p}; \qquad p(\sigma_i) \approx 1.9 + 0.7\sigma_i^{-1/2}, \qquad (46)$$

where $\epsilon$ denotes the kinetic energy of the electrons without their rest-mass (see Fig. 26, lower left panel). They emphasize that this index can be understood on the basis of the
bouncing of electrons between approaching islands (see next section). The normalized cutoff energy $\epsilon_c/\sigma_e m_e c^2$ rises slowly with $\sigma_i$ in the semi-relativistic regime, from around 2.5 to 4 or 4.5 as $\sigma_i$ goes from 0.1 to 10. The authors emphasize that it is not yet clear, whether the computed power-law indices and cut-off energies are truly independent of the simulation length $L$ of the current sheet. In the ultra-relativistic regime, $\sigma_i > 10$, the ion spectra show a power law which closely matches that of the electrons. For lower $\sigma_i$, the situation is somehow puzzling (see Fig. 26, upper row). A possible power law only appears at high ion energies while at lower energies the spectrum is rather flat and much harder, with a slope of about 1. This part is neither a clear power law nor, as it is much broader, a Maxwellian. This flat region in the spectrum turns downward significantly below the electron cut-off energy. However, at these energies, there are always more ions than electrons. The authors find no explanation for this behavior but speculate that indeed the power laws may show a break at energies where protons become trans-relativistic ($\omega \approx m_i c^2 \approx 103$ MeV). We add that yet another possibility is different dominant acceleration process for protons depending on energy.

Ball et al. (2018) present another study of collisionless relativistic REC of an electron-ion plasma, adding new facets to the picture (see Fig. 27). Firstly, the dependence on the $\beta_i$ parameter is systematically explored, in addition to the dependence on $\sigma_i$. They explore $\sigma_i = 0.1, 0.3, 1$ and $3$, while $\beta_{p,i}$ varies from $\beta_{p,i} = \beta_{p,e} = 10^{-4}$ up to the maximum possible value of $\beta_{max} \approx 1/4\sigma_i$. A realistic mass ratio is used. The computational domains are larger by at least a factor 5 than those of previous studies. The initial current sheet width is $\delta = 80c/\omega_{p,e}$. Periodic boundaries are used normal to the sheet. A moving injector and a dynamically-enlarging box is used. A description of the implementation of this boundary type is given in Sironi and Spitkovsky (2009). In this way, the magnetic flux is not limited to the one present at initial times. Finally, they trigger the start of REC by removing the pressure at one point within the sheet's center. Two exhausts with powerful depolarization fronts develop. As the boundaries are periodic normal to the sheet, the fronts meet to form a large plasmoid. When the fronts first meet, a secondary current sheet is created normal to the first sheet. In addition, they perform a study with different initial and boundary conditions, see below and Fig. 27, lower panels. All runs of this study are integrated till $2 t_A$, independently of the box size and value of other relevant parameters, with $t_A$ is given by the ratio between the length of the sheet and the Alfvén speed of the inflow, $t_A = L_x/U_A$.

The initial current sheet fragments to a certain degree into secondary plasmoids and secondary X-points. For fixed $\sigma_i$, the fragmentation is less pronounced for higher $\beta_{p,i}$ and, for fixed $\beta_{p,i}$, is more pronounced for higher $\sigma_i$. Sironi et al. (2016) found the same dependences for the ultra-relativistic case.

Ball et al. (2018) find that the electron spectrum in the reconnection region is non-thermal and can be modeled as a power-law with slope $p$ which depends on $\sigma_i$ and $\beta_{p,i}$, as

$$p(\sigma_i, \beta_{p,i}) = A_p + B_p \tanh(C_p \beta_{p,i}) \\ A_p = 1.8 + 0.7/\sqrt{\sigma}, \quad B_p = 3.7\sigma_i^{-0.19}; \quad C_p = 23.4\sigma_i^{0.26}. \qquad (47)$$
Thus, at low $\beta_{p,i}$, the slope is (nearly) independent of $\beta_{p,i}$ and hardens with increasing $\sigma_i$, having a (nearly) equal form as found by Werner et al. (2018). At higher values of $\beta_{p,i}$, the electron power law steepens and the electron spectrum eventually approaches a Maxwellian distribution for all values of $\sigma$, see Fig. 27, upper left panel. At values of $\beta_{p,i}$ near $\beta_{p,i,\max} \approx 1/4\sigma_i$, when both electrons and protons are relativistically hot prior to REC, the spectra of both species display an additional component at high energies, containing a few percent of particles, see Fig. 27, upper right panel.

Importantly, when using the same $\sigma_i$ and $\beta_{p,i}$ as Werner et al. (2018), the non-thermal particle spectra found by Ball et al. (2018) are systematically softer than the one by Werner et al. (2018). It is shown in Ball et al. (2018) that this discrepancy may be caused by two reasons: firstly, the authors find numerically that a larger box size makes the spectra generally softer—though there are weak indications for a certain saturation of the slopes at the two biggest box-lengths they have used ($L_x c/\omega_{p,e} = 5'440$ and 10'880). Secondly, also found numerically, runs with a (single) triggered initial X-point show a softer spectrum than runs in which the tearing instability produces spontaneously many X-points (as used by Melzani et al. 2014b and Werner et al. 2018). Indeed, Ball et al. (2018) reproduce exactly the same slopes as Werner et al. (2018) when using exactly the same setup and boundary conditions. The lower row of Fig. 27 summarizes the influence of the initial and boundary condition for the electron spectrum. Another finding is that protons have up to a magnitude larger mean energy than electrons, though protons show a steeper slope in their spectrum than electrons.

An important result is that for all low $\beta_{p,i}$, the time-evolution of non-thermal acceleration is different for electrons and protons. While electrons immediately develop a non-thermal tail in their spectrum, protons develop a non-thermal tail only after $t \approx 0.8t_A$, corresponding approximately to the time when the two reconnection fronts interact across the periodic boundaries. This may indicate another acceleration mechanism for the two species.

**Result summary** The last few years have brought progress in our understanding of the population of the non-thermal high-energy particles accelerated by relativistic REC. However, we emphasize again that all these results have been achieved on the basis of only one particular setting, 2D Harris or force free sheets. One should always keep in mind that this setting is a very particular one out of many other, probably more realistic settings. Under these conditions, the reconnection rate is given by Eq. (44). This rate is about 0.2 and thus higher than in the non-relativistic case. Two parameters are responsible for the formation of the power-law slope of the distribution of non-thermal particles accelerated in the reconnection event. The ‘cold’ magnetization $\sigma_i$ and the magnetization $\beta_{p,i}$, which includes thermal aspects of the plasma in which the REC takes place. A pretty good expression for this power-law slope is given by Eq. (47). Note that for cold plasma (low $\beta_{p,i}$), this expression becomes fairly independent of $\beta_{p,i}$ and matches the simpler expression given by Eq. (46). Without a background field—unless in the highly relativistic limit—the ions gain more energy, up to 75% of the released magnetic energy in the semi-relativistic regime (see Eq. (45)).

There are, however, some aspects which disturb the picture and point to the need of additional work. Firstly, the exact shape of the power-law seems to depend, to a certain degree on the initial and boundary conditions (Fig. 27, lower row) and on the
length of the Harris sheet. Secondly, if electrons and ions are hot before reconnection, there seems to exist an additional power-law slope in the electron spectra, at least at late times (Fig. 27, upper right panel). This may indicate that at least two different acceleration mechanisms are dominantly at work.

While the questions just addressed can be judged as minor, there remain two more fundamental largely open questions. One is the role of a background field (and such a background field is always present in a real environment). Melzani et al. (2014b) found indications that a background field (at least up to $B_G = B_0$) does not affect the power-law slope of the electron distribution, but the spectrum evolves more slowly than without a background field (Fig. 25, right column panels). In addition, this study found that in the presence of a background field, electrons may gain more energy than protons, reversing the ratio found in simulations without a background field. However, *one study is no study*, the more as the low mass ratio used for these simulations may have spoiled the result. The second question is the exact shape of the proton population distribution. At low $\sigma_i$, these spectra seem to consists of a very flat part at intermediate energies, before a clear power-law is established at higher energies (Fig. 26, upper left panel). This may again point to two different acceleration processes at work. The slope of the high-energy power-law for protons may be close to the one of electrons (Melzani et al. 2014b; Werner et al. 2018) or may be steeper (Ball et al. 2018):

### 4.2.3 Selected physical processes which accelerate particles

So far there is no comprehensive picture on which physical process is responsible for the acceleration of non-thermal particles in magnetic reconnection—though several reasonable ideas exist and some processes could be identified to be active, but possibly together with other processes. Therefore, at this place, we present those acceleration channels most discussed in the literature and the arguments of the authors who advocate them. We attempt a provisionally ranking in the next section.

**Thermal exhausts** As discussed in Sect. 2.5, REC at single X-points or on larger current-sheets produces exhausts where the particles leave the reconnection region (consult Fig. 10). Strictly speaking, the flow in the exhaust is not completely thermal. For instance, temperature is non-isotropic and electrons and ions do not have the same temperature, see e.g., Fig. 10 of Melzani et al. (2014a). However, the magnitude of the speed in the exhausts, $u_{out} \sim \sqrt{2} u_{A,in}$ (Eq. 15), and the temperature in the exhausts can be understood on the basis of flow conservation laws between the inflow and the outflow.

Note that in a highly magnetized environment such as relativistic REC (but not exclusively), the Lorentz-factor of the exhausts may already be of order of a few, because the Alfvén speed of the inflow may already be very close to the speed of light.

Below, the dynamics of chains of X-points and current sheets are discussed. There, exhausts of the different sites collide and form plasmons (Fig. 24). But in nature there are also reconnection sites found with just one X-point. There, the exhausts can freely expand to form jets. This mechanism is used to explain jets on all scales, from the solar atmosphere and corona (Zharkova et al. 2011), but also from compact objects (de Gouveia dal Pino and Lazarian 2005). The more, small reconnection events and
associated exhaust within large scale jets may be at the origin of fast TeV variability in blazars (Giannios et al. 2009; Khiali et al. 2015).

**Dynamics of plasmons and chains of plasmons** We now look at a longer current sheet that break apart and, consequently, allows a variety of potential particle acceleration mechanisms—to be further detailed below—to act on the plasma. If the sheet is long enough, many X-points will develop, with associated exhausts. Exhausts of neighboring X-points collide and form plasmons, plasma regions which are bound by a strong circularly closed magnetic field (consult Fig. 24). The formation of plasmons (→ Colliding plasmons) has the potential to accelerate particles as the border of two associated exhausts, called dipolarization fronts, form approaching magnetic mirrors, see e.g. Lapenta et al. (2015).

The field within the plasmoid tends to zero. Thermal plasma within the plasmons cannot escape as the strong encircling field deflects particles immediately back to the interior. As REC goes on and more of the inflowing plasma is being processed, islands potentially grow in size and their encircling fields grow in strength. Plasmons will also merge and grow in size. After merging, plasmons will contract and give rise to an important acceleration mechanism (→ Contracting plasmoid). The process of breaking a current sheet will also generate turbulence, another important source of particle acceleration (→ Turbulence).

At this place, the stability of current sheets, the formation and dynamics of plas-
moids cannot be reviewed in detail. This process sets the stage for the acceleration of
ultra-fast particles, but not directly causes it. For a deeper understanding of the process,
the reader may consult the vast literature on the subject, e.g. Loureiro et al. (2007),
Lapenta (2008), Samtaney et al. (2009), Lapenta and Lazarian (2012), Kowal et al.
(2012), Markidis et al. (2013), Kagan et al. (2013), Loureiro and Uzdensky (2016),
Sironi et al. (2016), Kowal et al. (2017).

**Contracting plasmons** PIC simulations show that contracting plasmons (see structure B in Fig. 24) can efficiently accelerate electrons (Drake et al. 2006) and ions (Drake et al. 2010) to super-Alfvénic speeds by a Fermi-like process. Such particles are injected into magnetic islands from the exhausts of X-points of a reconnecting current sheet. Particles approaching the strong magnetic fields encircling magnetic islands will be turned around by these fields and are captured in this way within the island, bouncing in the island backwards and forwards or move along the magnetic field lines. When islands are freshly formed, the encircling field also lets contract the island, up to a point where a quasi-stationary equilibrium is reached between the pressure of the enclosed particles and the tension force of the magnetic structure. In this way, the particle will gain energy at each bounce. When more energized, the particles will start to diffuse out of the island and either escape or, in a more complex situation, may be kept by another island held together with a stronger field.

Counting on this mechanism, Drake et al. (2010) can explain the anomalous cosmic ray (ACR) energy spectrum observed by both voyager missions in the region between the solar wind termination shock and the heliopause (note that the same idea was also developed in Lazarian and Opher 2009). Drake et al. (2010) firstly perform large scale MHD simulations of the solar wind which show stripes of inversed field directions
due to the non-alignment of the magnetic dipole of the Sun with its rotation axis. Current sheets develop in the region where the field polarity changes. (This is very similar to the striped neutron star winds (Kirk 2004; Uzdensky et al. 2011; Cerutti et al. 2012a, 2014a; Bühler and Blandford 2014).) The MHD simulations show that stripes and current sheets are compressed in the passage of the solar wind through the solar wind termination shock, setting the stage for REC in the now unstable sheets. The post-shock situation can be well approximated by multi-layered Harris-sheets.

Subsequent PIC simulations using the plasma-parameters found in the MHD simulations demonstrate that particles can be accelerated to several tens Alfvénic speeds, forming a power law with a spectral index of about 1.5. The multi-layered Harris-sheets break off, forming a network of plasmoids which collide and collapse. By a thorough analysis of the data, Drake et al. (2010) show the outstanding role contracting islands play in the acceleration process.

The model assumes that interstellar ions are picked up by the solar wind and advected subsequently through the solar wind termination shock. Thus the model can also explain similarities in the spectra of different ion species observed by the voyager missions. The result achieved demonstrates the power of combining large-scale MHD simulations with detailed kinetic simulations.

Subsequently, Kowal et al. (2011) performed MHD simulations in two and three space dimensions of the same multi-layered Harris sheets. A dynamic network of magnetic islands develops, similar to the one observed in the PIC simulations of Drake et al. (2010). Test-particles are injected into this configuration and integrated with the sixth-order implicit Runge-Kutta-Gauss method (Sanz-Serna and Calvo 1994) which conserves particle energy and momentum. The accelerated test-particles show also a super-Alfvénic distribution though there are differences in the distribution as compared to Drake et al. (2010). Contracting island again play a dominant role in the acceleration process. But the authors report also on drift acceleration along the magnetic field discontinuity and the Fermi-process between the converging inflows (see the corresponding paragraphs within this section). The reasons behind these differences remain open. The paper extends the multi-layered Harris sheet configuration to three space dimensions and repeats the analysis. They find substantial differences in the acceleration process and stress the need for more comparative studies between 2D and 3D settings.

Montag et al. (2017) develop a general framework which includes compressibility and non-uniform fields to analyze the process of contracting plasmoids. They derive an expression for the power-law scaling of the distribution function and for the rate at which the power-law develops in time. In analogy to the case of acceleration by the converging inflows, the spectrum gets harder if the compression increases, which is generally true for Fermi-like processes. The authors also find that a guide field of order unity suppresses the development of power-law distributions.

**Colliding plasmoids** As described, plasmoids eventually collide. Subsequently, they first grow and then contract, thereby accelerating particles.

There are two other interesting points in colliding plasmoids. Firstly, at the contact interface between the coalescing plasmoids secondary current-sheets with secondary REC develop, mostly normal to the primary sheet. This secondary REC may support or
suppress particle acceleration, its effect seems, however, to be small, see the discussion in Bessho and Bhattacharjee (2012). Secondly, even before collision, in the phase they approach each other, a Fermi-like acceleration process may work as particles are reflected on the plasmons and travel back and forth between neighboring approaching plasmons (Lapenta et al. 2015, for instance).

Werner et al. (2018) use this process to explain the expression for the power-law index they found (Eq. 45). The acceleration by approaching plasmons is a second order Fermi process as the movement of the plasmons can be regarded as stochastic. The energy gain per bounce is thus $\Delta\epsilon \simeq (u_{pl}/c)^2\epsilon$, where $\epsilon$ is the particle energy without the rest-mass energy and $u_{pl}$ the speed of the plasmoid. The typical bounce time can be approximated by $t_b \sim \lambda_{pl}/c$. Thus, the acceleration time-scale is $t_{acc} \equiv \epsilon\Delta t_b/\Delta\epsilon = \text{const} \cdot c\lambda_{pl}/u_{pl}^2$. The process ends when the two plasmons collide and the particle escapes, at $t_{esc} \simeq \lambda_{pl}/u_{pl}$, where $\lambda_{pl}$ denotes a typical distance between plasmons. By this, the Fermi II power law index, $p$, can be written as $p = 1 + t_{acc}/t_{esc} = 1 + \text{const } c/u_{pl}$. Werner et al. (2018) assume that $u_{pl} \simeq U_A$. In the regime where $\sigma_i \ll 1$ (but still $\sigma_e \gg 1$), $U_A \approx c\sigma_i^{1/2}$, and thus $p = 1 + C\sigma_i^{-1/2}$, which has the form of Eq. (45). We note that one can repeat the exercise for any other Fermi II process taking place within the sheet, in particular also for acceleration due to turbulence. Ball et al. (2018) account for a variant of this process to explain the additional component in the power-law at high energies when the inflow is hot ($\beta_{p,i}$ close to $\beta_{p,i,max}$).

**Turbulence** As discussed in Sect. 2.5, the break apart of current sheets by tearing mode instabilities (or a combination of tearing and kink mode in 3 spatial dimensions) generates turbulence in the region where REC takes place. This turbulence, in combination with turbulence in the inflows, is responsible that collisional REC becomes as fast as observed (Kowal et al. 2009; Lapenta and Lazarian 2012). In the frame of test particles, it was shown that MHD fast and slow compressible modes accelerate energetic particles through the second order Fermi acceleration. Density fluctuations in converging flows can enable first order Fermi acceleration of particles (Kowal et al. 2012). Kowal et al. (2017) provide a thorough analysis of the statistics of the reconnection-driven MHD turbulence. Lazarian et al. (2012) and Lazarian et al. (2016) review the development and the character of (self-generated) turbulence related to magnetic reconnection, in the classic as well as in the relativistic regime.

Corresponding studies for the collisionless case only recently start to appear. Zhdankin et al. (2018) study forced turbulence using PIC simulations. They show that particles can indeed be accelerated to very high energy also in the case of collisionless turbulence. Figure 28, top panel, shows the time-evolution of the spectrum of accelerated particles, converging to a power law of the form $\gamma^{-3}$. The lower panel shows the dependence of the domain resolution, from 512³ to 1536³ cubed.

In a preprint, Zhdankin et al. (2019) present corresponding results for an ion-electron plasma and come to the conclusion that non-thermal acceleration is efficient for both species in the fully relativistic regime. In the semi-relativistic regime, while still efficient for cosmic ray production, the mechanism is less efficient for electrons. For the production of hard non-thermal radiative signatures a (very) low $\beta_{p,i}$ is necessary, or nearly relativistic ion-temperatures. The authors emphasize that the result is still somehow preliminary as it is not yet established how this result scales to a large
<figure><img src="image_48.png" /><figcaption>**Fig. 28** Top panel: evolution of the particle Lorentz factor (energy) distribution $f(\gamma)$ for the $1536^3$ simulation of forced collisionless turbulence. The dotted line indicates the slope expected for a second order Fermi-process. Bottom panel: compensated distribution $f(\gamma)\gamma^3$, at fixed time $tU_{A0}/L = 7.0$, for varying system sizes. Power laws with pre-compensated index $-3.0$ (black dashed) and $-2.7$ (black dash-dotted) are also shown, along with the mean energy $\langle\gamma\rangle$ (green dash-dotted) and system-size cutoff $\gamma_{max}$ (green dotted) for the $1536^3$ case. Images reproduced with permission from Zhdankin et al. (2018), copyright by AAS</figcaption></figure>

system size. Related to these simulations, the authors have statistically analyzed collisionless turbulence (Zhdankin et al. 2017) and make a connection to the Fokker-Planck framework (Wong et al. 2019).

A word of caution, equally raised by the authors, must be added. The forcing of the turbulence implies a steady energy input into the system. As no energy sinks are present, energy may pile up to give artificial heating and non-thermal particle acceleration. More studies are clearly needed. Also, non-global compression was present which would rise the value for the power-law index (Bell 1978a; Drury 2012).

Nevertheless, it can be firmly stated that turbulence is present in reconnection sites, both, collisional and non-collisional, and that this turbulence can contribute to particle acceleration. We note these results are also very promising to advance our understanding of particle acceleration in shocks.
**The reconnecting electric field** This field is an important source of particle energization in the relativistic case.

As particles are de-magnetized in the current sheet, this field can directly accelerate particles without perturbation by other structures. However, the width of the region where the electric field is strong is not too large, i.e., it is given by $E > cB$, or with a guide field, by $E \cdot B \neq 0$.

In non-relativistic REC, this process cannot be too effective as the electric field is too weak as compared to the width of the acceleration zone and the particles escape soon (Drake et al. 2010; Kowal et al. 2011; Drury 2012).

In the relativistic case, however, the process becomes important and the accelerated particles will form a power law. Particles enter the region with a large electric field at all distances from the central X-points. The particles close to the X-point will only slowly be turned, by $B_x$, towards the exhausts, in contrast to those particles which enter the sheet at some distance from the X-point. The longer the particle is held in the central region, the longer it can be accelerated by the reconnection electric field. As emphasized by Zenitani and Hoshino (2001) the gyro-radius of relativistic particles grows with their Lorentz-factors. Already accelerated particles thus stay even longer within the region of a large electric field and hard tails of the distribution can grow. These authors predict, for a pair plasma, a power-law distribution of non-thermal particles created by the reconnection electric field of $p \propto cB_x/E_{Rec}$ Bessho and Bhattacharjee (2012) refined the analysis in two aspects, also for a pair plasma. Firstly, they looked closer at the creation of the power law distribution but then included also the process when the particles finally have to leave the region immediately around the X-point. The latter produces an exponential cut-off such that the particle distribution function can be written in the form

$$f(\gamma) \sim \gamma^{-1/4} \exp -a\sqrt{\gamma}, \qquad (48)$$

with a a constant of order $cB_x/E_{Rec}$. Note that this is a very hard distribution. The authors emphasize that particles accelerated by the reconnection electric field may be further accelerated by other processes, in the field when they swirl around contracting islands or when they encounter other X-points. We add that any other of the processes described in this section could be tapped. PIC simulations of relativistic collisionless REC confirm that this process significantly contributes to the high-energy spectrum, for both cases, of a electron/positron pair plasma (Bessho and Bhattacharjee 2012; Sironi and Spitkovsky 2014) and an ion-electron (Melzani et al. 2014b; Ball et al. 2018).

**Fermi process within the inflow region** Recall that the inverse field components float, due to $E \times B$-drift, towards the region where they will reconnect (Fig. 1, right panel).

de Gouveia dal Pino and Lazarian (2005) point out that this may, in principle, accelerate particles by a first order Fermi-process (see Sect. 2.2.1), if particles are allowed to cross the reconnection region and to be ping-ponged between the lower to the upper inflow into the sheet. Repeating the argumentation given by Bell (1978b) for
the shock case, de Gouveia dal Pino and Lazarian (2005) derive an emerging spectrum of accelerated particles of $N(E) \propto E^{-5/2}$ for this acceleration mechanism.

An important point is raised by Giannios (2010) who states that this mechanism does not depend on the not well known reconnection mechanism but works solely on the converging inflows to the reconnection region. Be $\delta$ the width of the dissipation region and $R_G$ the gyro-radius associated to the undisturbed magnetic field and $\theta$ the angle under which the particle enters the reconnection layer. Be this angle moderately small (as compared to the length of the sheet, $L$) and consider the case where the particle can freely pass the reconnection region. Then, the particle will gyrate around the field drifting into the reconnection site and will cross back the sheet to its other side, where it gyrates back and forth. Of course, turbulent field fluctuations or collisions may scatter the particle as well, but their presence is not really necessary. In addition, the particle is gaining energy (though less) at each passage in the reconnection electric field. Giannios (2010) refers to this process as the betatron effect.

The energy amplification factor, $A$, for one cycle can then be derived to:

$$A(\theta) = \Gamma_r^2 (1 + \beta_r \cos \theta)^2, \qquad (49)$$

where $\theta$ denotes the angle at which the particle enters the reconnection layer and $\Gamma_r$, $\beta_r$ the Lorentz-factor and normalized speed of the inflow. We note that with $\theta = 0$ and $\beta_r \ll 1$, the formula of de Gouveia dal Pino and Lazarian (2005) is recovered. The time-scale of the acceleration process is

$$t_{Acc} = \frac{2\pi\gamma mc^2}{(1 - 1/A)eBc} \qquad (50)$$

and is thus of the order of the gyration period of the largest $\gamma_e$ factor. Giannios (2010) shows that for an isotropic particle distribution (and thus $\theta$), the amplification amplitude, $A$, is only about one quarter smaller as for the case $\theta = 0$.

The process is limited by two factors, the size of the current sheet and by radiative losses. These are dominated by synchrotron losses while photo-pion production turn out to be much less important (Waxman 1995). Assuming jets carrying field-reversed components of length of the jet diameter, Giannios (2010) estimates that protons can reach energies up to $10^{20}$ eV in GRBs or luminous AGN jets. For iron, an important species in very high energy cosmic rays, the same limit applies. For electrons, however, the acceleration process stops much earlier due to synchrotron losses. The author emphasizes that an important pre-requisite for the process to work is that it can only be initiated if the particles are sufficiently pre-accelerated, because, otherwise, the gyro-radius of the particles is not large enough to let the particle stream over the reconnection region.

Drury (2012) points out that the spectrum derived by de Gouveia dal Pino and Lazarian (2005) is probably much too weak (i.e., the power-law slope should be much smaller than 5/2). This, because the compression of the plasma within the sheet is not considered. The paper provides arguments that, in the case of REC, the compression is probably larger than in the shock case, where mostly a compression ratio of 4 is assumed, which corresponds to an adiabatic shock. For REC, spectra like $f(p) \propto p^{-2}$
or $N(E) \propto E^{-1}$ have to be expected. This was the first indication that spectra of particles accelerated in magnetic reconnection events may indeed be very hard, much harder than for particles accelerated by the Fermi process in shocks.

Bosch-Ramon (2012), considering compressible effects, realistic cross B-field diffusion coefficients, and accounting for synchrotron cooling, shows that for protons a maximum energy of

$$E_{\max} \approx 60 (\chi/k)^{1/2} (v/c)^{1/2} \frac{1}{\sqrt{B_0}} \text{TeV.} \qquad (51)$$

can be reached. k and $\chi$ depend on the concrete microphysics, with good arguments for $\chi/k \approx 1/10$. The spectrum of the accelerated particles is harder than standard Fermi I.

Provornikova et al. (2016) analyze possible compression rates of REC in the solar corona. They look at configurations such as a Harris current sheet, a force-free current sheet, and two merging flux ropes. Plasma parameters are taken to be characteristic of the solar corona. It is found that plasma compression is expected to be strongest in low-beta plasma $\beta_p \sim 0.01-0.07$ at reconnection magnetic nulls and can be as high as a factor 10. Provornikova et al. (2018) derive magnetic wave solutions in compressible reconnection sites and claim that such waves may act as scattering centers for the Drury-mechanism.

In conclusion one may say that this is a promising process whenever there are pre-accelerated particles present which can initiate the back- and forth-bouncing.

**Drift acceleration** If the particles enter the diffusion region or if they leave the current sheet the magnetic field changes substantially. In this situation, the particles encounter magnetic drifts and magnetic drift acceleration (see Sect. 2.3 for the description of the equivalent mechanism in shocks). Kowal et al. (2009) and Guo et al. (2016) report that this mechanism contributes to the acceleration of particles, both in the collisional (Kowal et al. 2009) and non-collisional (Guo et al. 2016) case.

**Associated shocks** de Gouveia dal Pino and Lazarian (2005) point to another interesting possibility: analogous to the sun, eruptive reconnection may push matter away and produce a shock-wave in which particles can be accelerated by the 'ordinary' Fermi-mechanism related to shocks, leading to a power-law slope of 2.

**Dimensionality** The role of dimensionality in the particle acceleration process is largely unexplored.

2D spatial setups can easily be extended to 3D spatial setups by just expanding the 2D configuration to the third dimension in a planar way. Such configurations will still be subject of the tearing mode, but, in addition, the 3D extension will be subject to the kink instability. Cerutti et al. (2014b), using PIC simulations with radiative feedback, compare the cases where either the kink or the tearing instability grows faster and discuss application to the Crab nebula. In their setup they found the kink mode to dominate, leading to a disruption of the current sheet and associated turbulence unless the kink-mode is stabilized by a background magnetic field. The same result is also
found by Oishi et al. (2015) in 3D MHD simulations. Other authors, using different configurations, found that rather the kink mode develops slower than the tearing mode even without a guide field. Here, REC develops essentially similarly as in two spatial dimension. Oblique modes, a combination of tearing and kink modes, are possible and may even grow fastest (Daughton et al. 2011).

The results of different 3D kinetic studies of current-sheet reconnection do not yet converge to a unique picture. Some authors claim that there are small differences in the reconnection rate between 2D and 3D, e.g., (Liu et al. 2012; Daughton et al. 2014, for non-relativistic REC) and (Guo et al. 2014a, for the relativistic case). On the other hand Sironi and Spitkovsky (2014) found a four times lower rate for 3D REC of a relativistic pair plasma as compared to 2D REC.

In a series of papers it was advocated that the drift kink instability can modify the electric and magnetic field structures in an anti-parallel reconnection layer and prohibit non-thermal acceleration (Zenitani and Hoshino 2005a, b, 2007, 2008). But other authors found that the kink mode cannot suppress the acceleration of particles (Liu et al. 2011; Sironi and Spitkovsky 2014).

There are very few generic 3D configurations which include 3D nulls. Based on observations, Baumann and Nordlund (2012) reconstruct a field configuration within the solar corona and use it as initial condition for 3D PIC simulations. Olshevsky et al. (2013) simulate REC, starting from a cluster of eight null points. Much more work will be necessary in future to get a more comprehensive picture of 3D reconnection events.

### 4.2.4 A critical discussion and outlook

The first point to state is that there has been a tremendous progress in our understanding of magnetic reconnection in the last few years. Given a correct environment, both collisional and non-collisional magnetic reconnection proof to be fast. For the collisional case, both, self-generated and/or external turbulence is the key-ingredient to make the process fast. The collisionless case turns out to be always fast. With this progress, we now can understand qualitatively, and even to a good degree quantitatively the ‘thermal aspects’ of REC.

This review has for subject the non-thermal, ultra-energetic particles. Mostly PIC-simulations have shown that such particles can originate from reconnection sites. Many questions remain, however. We want to address three of them, A–C, in the following paragraphs.

(A) Open questions concerning small scale kinetic simulations As summarized at the end of Sect. 4.2.2, we have, on the basis of 2D current sheet simulations some definitive results. However, and as emphasized also there, many details remain unanswered also for this case. There remains, in our view, three large questions which need to be answered before we can state that our understanding even for this simple case is sufficiently secured.

1. There is only one study which includes a guide field (at least for ion/electrons). This study has brought interesting results, the independence of the spectral slope on
the presence and the strength of a guide field and that the energy partition between electrons and ions may critically depend on the presence of a guide field. Guide fields are necessarily present in a current sheet like reconnection event in space. Corresponding large parameter studies, which derive the relations for ion-electron energy partition, the dependence of the power-law slope on the magnetization and the thermal state of the inflow when a guide field is present are urgently needed.

2. The situation in three spatial dimensions remains unclear, even for the most simple case, the extension of a 2D current sheet towards the third dimension and a comprehensive study is lacking.

Beyond this simple case, there are much more complicated reconnection topologies present in 3D than the simple extension of a 2D current sheet to 3D. An overview of such topologies can be found in Birn and Priest (2007, Chapter 2) and Pontin (2011). None of these topologies have been systematically addressed by kinetic simulations.

3. Finally, REC is always accompanied with radiative emission. In a magnetized environment, synchrotron radiation is always present and cool the particles which emit. In addition, these photons undergo inverse Compton scattering with the hot particles, cooling them again. Many reconnection sites are embedded in a strong external cold radiative source, e.g., companion stars in X-ray binaries, leading to additional cooling. Radiative effects introduce new parameters into REC not discussed here so far. For instance, REC in a micro-quasar corona close to the hole and in the γ-ray emitting region of an extragalactic jet takes place with the same magnetizations (see Table 3). Without considering radiative effects, REC in both environments features the same reconnection rate, particle spectra, or energy distribution. However, the field strength differs by six orders of magnitude and thus the effect of synchrotron radiation is largely different.

There are first attempts to account for radiative emission in PIC-simulations. Werner et al. (2018) perform a relativistic 2D Harris sheet study of a pair plasma and describe effects of external inverse Compton cooling on the basic dynamics, the non-thermal particle acceleration, and radiative signatures. They find the reconnection rate and the overall dynamics basically unchanged. Important differences are found for the particle spectra. They still show a hard power law (index $\ge -2$) as in nonradiative REC, but transition to a steeper power law that extends to a cooling-dependent cutoff. The steep power-law index fluctuates in time between roughly -3 and -5. Some other studies which include radiative losses address mostly REC in pulsar winds (Cerutti et al. 2014b, a, 2016). Also the community of the laser-plasma facilities starts to study radiative effects with PIC codes as newer, more powerful lasers establishes a regime where radiative losses essentially co-determine the dynamics, see for instance Wallin et al. (2015), Gonoskov and Marklund (2018), Blackburn and Marklund (2018).

**(B) How are particles accelerated?** So far, it is not possible to rank the relative importance of the different acceleration processes discussed in Sect. 4.2.3. It seems likely that different mechanisms are at work even for a simple setting, but the more if one also considers the prevailing large scale physical conditions, e.g., hot non-relativistic pair plasma or cold ultra-relativistic electron proton plasma (see also part C below).

<figure><img src="image_49.png" /></figure>
Most papers mentioned above come to this conclusion. In the same direction points the finding that the distribution functions of the non-thermal particles, electrons and protons, often show broken power-law slopes, with up to three different slopes and an exponential cut-off. Even if some studies suggest that the shape of the distribution functions for protons and electrons are close to identical, other studies shows differences between the shapes of electron and proton distribution functions. In relativistic REC, the direct acceleration by the reconnection electric field is unambiguously identified as one important ingredient—in contrast to the non-relativistic case. However, all studies showed that it is accompanied by some stochastic Fermi-process. Whether this process operates between magnetic islands or other magnetic structures or is just a consequence of kinetic turbulence in the diffusion region and in the current sheet, is not yet clear.

First order Fermi processes can also be present. Such a process certainly works in contracting islands and seems to accelerate the population of particles which is originally located in the current sheet. If particles can cross island-boundaries, this process may contribute to the acceleration of other particles as well. Whether, or under what physical conditions, the first order Fermi process between the converging inflows can work is an open question. In none of the relativistic simulations it has been observed. The reason may be that a Harris sheet is a too symmetric constellation where it is very hard to initiate particle motions normal to the direction of the sheets which are sufficiently fast to move the particle out of the diffusion region. However, if one would be able to find a mechanism to pre-accelerate particles, this mechanism could become operational. Corrugated sheets, spine-fan or other topologies may support the initialization of this mechanism. This question remains open and is part of the next point of discussion.

**(C) Magnetic reconnection in a large scale environment and in real objects** Under-standing microphysical processes is decisively important but, on the other hand, only one part of the game to understand REC and to what degree REC contributes to the emission of real objects via both, the leptonic and hadronic channel. And it does not answer the question whether REC contributes to the cosmic ray flux. To get the answers to such questions, we definitely need to combine microphysics with large scale flow.

It is important to understand what triggers a reconnection event and the nature of the event. Are, in real astrophysical objects, quasi-stationary configurations present, like Harris-sheets or other magnetic nulls which eventually get unstable? Or are reconnection events driven by large-scale flows more important. Related, does REC start at one point or will a network of reconnection points establish? All this is largely unexplored.

We therefore advocate simulations of large scale MHD flows, entire accretion disks and large portions of jets in high-energy objects. Such MHD simulations must control magnetic diffusion, to be able to identify reconnection sites and topologies. Ideally, such large scale simulations are also resistive. Only on the basis of such MHD simulations we will obtain a good estimate of the spatial and temporal distribution of reconnection events in such objects. Subsequent microphysical studies which adapt the configurations found by the MHD solutions will then allow to model photon and possible neutrino emission. MHD simulations of magnetic reconnection will be reviewed in a future version of the review.
# 5 Macro scale numerical particle acceleration studies

At scales comparable with the system size dynamics are often treated in the magnetohydrodynamics (MHD) approximation. In this approximation different methods have been developed to handle the acceleration and propagation of energetic supra-thermal particles, which are listed and described below. Before treating this aspect we present in Sects. 5.1 and 5.2 the main MHD solvers used in most of modern codes. A discussion concerning relativistic MHD is also included. Detailed monographs and reviews on the subject can be found, e.g., in Leveque (1998), Martí and Müller (2015) and references therein. The next sections treat the way CRs can be coupled with MHD. Section 5.3 discusses the multi-fluid approach where energetic particles are treated as a fluid. Section 5.4 describes the procedure to combine kinetic and HD/MHD methods, in particular the way to treat energetic particles back-reaction over fluid solutions (see Sect. 3.6). Section 5.5 presents the P(MHD)IC method in some details, this method combines simulation techniques exposed in Sect. 3.3 to investigate microscale physics but uses the electromagnetic field derived from the MHD equations. Section 5.7 discusses semi-analytical calculations of the problem of DSA, introduced in Sect. 2.2.

## 5.1 The equations of magnetohydrodynamics

The use of a magnetohydrodynamic (MHD) approach is crucial for the description of the large scale astrophysical phenomena involving collisional plasma, i.e. if the timescale associated to collision is shorter than the system dynamical time. This circumstance can occur for instance in supernova remnants, stellar bubbles and accretion discs. However, most of astrophysical shocks (including SNR shocks) are non-collisional as previously explained and their structure can not be dynamically described using a MHD model. However, if collective interactions are sufficiently frequent to keep the system isotropic, if electroneutrality can be assumed (which is the case for scales larger than the Debye length), in the cold plasma approximation the MHD equations reduce to the one-fluid system of equation derived in Eq. 52 below. This system describes long wavelength and low frequency perturbations.

### 5.1.1 Classical magnetohydrodynamics

MHD equations couple fluid mechanics equations and Maxwell's equations. They are obtained by averaging the moments of the Boltzmann equation over the velocity space. Particle density conservation is deduced by taking the zeroth order moment of the Boltzmann equation, the momentum conservation equation is obtained by taking the first order moment, and the energy conservation equation is obtained by taking the second order moment. Each of these moment equations introduces a new unknown function: the continuity equation introduces the velocity, the moment conservation equation introduces the pressure and the energy conservation equation introduces the internal energy. Additional assumptions are then required to close the system. To that aim, an equation of state is usually introduced which expresses the internal energy as function of density and pressure, but other assumptions can be used concerning energy fluxes.
The full MHD equations for a single fluid are given in a conservative form as follows:

$$ 
\begin{aligned}
\partial_t \rho + \nabla \cdot (\rho \mathbf{u}) &= 0, \\
\partial_t (\rho \mathbf{u}) + \nabla \cdot \left( \rho \mathbf{u}\mathbf{u} - \mathbf{B}\mathbf{B} + \left(p + \frac{B^2}{2}\right) \mathbf{I} \right) &= \mathbf{0}, \\
\partial_t (e) + \nabla \cdot \left( \left(e + p + \frac{B^2}{2}\right) \mathbf{u} - \left(\mathbf{B} \cdot \mathbf{u}\right) \mathbf{B} \right) &= 0, \\
\partial_t (\mathbf{B}) - \nabla \times (\mathbf{u} \times \mathbf{B}) &= \mathbf{0},
\end{aligned}
\quad (52) $$

where $\rho$ is proper rest mass density, $p$ is the thermal pressure, $\mathbf{u}$ is the Eulerian fluid velocity vector, $\mathbf{B}$ is the magnetic field, and $e = \frac{p}{\gamma_{ad}-1} + \frac{\rho u^2}{2} + \frac{B^2}{2}$ is the total energy density, where $\gamma_{ad}$ is the gas adiabatic index. The ideal Ohm's law (perfect conductivity) is retained: $\mathbf{E} = -\mathbf{u}/c \times \mathbf{B}$.

These equations are justified for a plasma where the relevant time scales are long in comparison with microscopic particle motion time scale and spatial scales are large in comparison with the thermal ion gyroradius and Debye length.

### 5.1.2 Equation of state (EOS)

The set of partial differential Eq. (52) governing compressible fluid dynamics is incomplete. There are more unknowns than equations and then an additional closure equation is required. Usually, this equation involves the internal energy, the thermal pressure, and the density, and in some cases the temperature as well. This closure equation describes the thermodynamical processes in the fluid. In a real plasma the different particle populations can have different thermodynamical behaviour. For simplicity, the most used closure equation is the polytropic equation of state,

$$ \frac{d}{dt} \left( \frac{p}{\rho \gamma_{ad}} \right) = 0, \quad (53) $$

in an adiabatic gas with one degree of freedom $\gamma_{ad} = 5/3$ and in an isothermal gas $\gamma_{ad} = 1$.

## 5.2 Numerical solutions

The finite volume method (FVM) is widely applied to solve problems described by hyperbolic partial differential equations (PDEs) as is the case for MHD equations. In this method the set of Eq. (52) is written in a shortened conservative form:

$$ \partial_t \mathbf{U} + \partial_i \mathbf{F}^i(\mathbf{U}) = \mathbf{0}, \quad (54) $$

where $i = (x, y, z)$ are the space variables, the conserved variables are given by
<figure><img src="image_50.png" /><figcaption>Fig. 29 Finite volume scheme with variables set at cell center and flux computed at face center</figcaption></figure>

$$ \mathbf{U} = \begin{bmatrix} \rho \\ \rho \mathbf{u} \\ e \\ \mathbf{B} \end{bmatrix}, \qquad (55) $$

and the associate fluxes are

$$ \mathbf{F} = \begin{bmatrix} \rho \mathbf{u} \\ \rho \mathbf{u} \mathbf{u} - \mathbf{B}\mathbf{B} + \left(p + \frac{B^2}{2}\right) \mathbf{I} \\ \left(e + p + \frac{B^2}{2}\right) \mathbf{u} - \left(\mathbf{B} \cdot \mathbf{u}\right) \mathbf{B} \\ \mathbf{u}\mathbf{B} - \mathbf{B}\mathbf{u} \end{bmatrix} \qquad (56) $$

Numerically Eq. (54) are discretized in space and time as

$$ \mathbf{U}^{n+1} = \mathbf{U}^n + \Delta t^n \left[ \sum_{idim=1}^{idim=ndim} \frac{\mathbf{F}_{idim}^{i+\frac{1}{2}} - \mathbf{F}_{idim}^{i-\frac{1}{2}}}{\Delta x_{idim}} \right], \qquad (57) $$

where $\mathbf{U}^n$ and $\mathbf{U}^{n+1}$ are the conserved variables respectively at time $t^n$ and $t^{n+1} = t^n + \Delta t^n$ where $\Delta t^n$ is the time step. The time-averaged fluxes $\mathbf{F}_{idim}^{i+\frac{1}{2}}$ in the time interval $[t^n, t^{n+1}]$ at the interface between a cell with indices i and its neighbour i + 1 in the direction idim are calculated from the solution of Riemann problems.<fnref n="45" /> This is illustrated in Fig. 29.

In MHD, the Riemann problem is described by a 7-waves pattern. These seven eigenvalues correspond to the left and right going Alfvèn waves and four magnetosonic waves (two fast and two slow), and between these two propagating waves there is the entropy wave. These waves are defined in a direction idim as

<footnote n="45">⁴⁵ A Riemann problem is an initial value problem of a conservative equation in fluid dynamics which involves a discontinuous distribution of the conserved variables. It leads to the derivation of characteristics or the eigenvalues of the problem corresponding to different wave solutions.</footnote>
$$ 
\begin{aligned}
\lambda_{2,6}^{\text{idim}} &= u^{\text{idim}} \mp U_{a}^{\text{idim}} \\
\lambda_{1,7}^{\text{idim}} &= u^{\text{idim}} \mp U_{f}^{\text{idim}} \\
\lambda_{3,5}^{\text{idim}} &= u^{\text{idim}} \mp U_{s}^{\text{idim}} \\
\lambda_{4}^{\text{idim}} &= u^{\text{idim}}
\end{aligned}
\qquad (58) $$

where the Alfvèn speed $U_a^{\text{idim}}$, the fast $U_f^{\text{idim}}$ and the slow $U_s^{\text{idim}}$ speed in direction idim in ideal MHD case are defined as (Ryu et al. 1995)

$$ 
\begin{aligned}
U_a^{\text{idim}} &= \frac{|B^{\text{idim}}|}{\sqrt{\rho}}, \\
\left(U_{f,s}^{\text{idim}}\right)^2 &= \frac{c_s^2 + U_a^2 \pm \sqrt{(c_s^2 + U_a^2)^2 - 4 c_s^2 U_a^{\text{idim}^2}}}{2}.
\end{aligned}
\qquad (59) $$

Here $c_s = \sqrt{\left(\frac{\partial p}{\partial \rho}\right)}$ is the local sound speed and $U_a = \frac{|B|}{\sqrt{\rho}}$ is the local Alfvèn speed. We can note here that the characteristic MHD waves are direction dependent. These waves are depicted in Fig. 30. The seven eigenvalues satisfy the inequalities

$$ \lambda_1 \le \lambda_2 \le \lambda_3 \le \lambda_4 \le \lambda_5 \le \lambda_6 \le \lambda_7. \qquad (60) $$

However, some eigenvalues can coincide depending on the direction and strength of the magnetic field. Therefore, the MHD equations form a non-strictly hyperbolic system (Brio and Wu 1988).

### 5.2.1 The Courant–Friedrichs–Lewy condition

Courant, Friedrichs and Lewy (1928) showed that the stability of numerical schemes requires the use of all the information contained in the initial state that will influence the solution in a given spatial cell. To satisfy this condition, the ratio between the spatial discretization $\Delta x$ and the time step $\Delta t$ should be smaller than the largest velocity of the signal solution of the PDEs; i.e. $\max(|\lambda_1|, |\lambda_7|)$ (maximum speed propagating to the left and to the right). This inequality is called the CFL condition:

$$ \Delta t \le \frac{\Delta x}{\max(|\lambda_1|, |\lambda_7|)}. \qquad (61) $$

Satisfying this condition is necessary for the convergence of explicit difference schemes.

### 5.2.2 Riemann solvers

The Riemann solvers are a fundamental tool in the development of FVM. They are based on a simplification of the hyperbolic equations. Moreover, at each step they simplify the physical state to constant piecewise values with jump discontinuities at some
<figure><img src="image_51.png" /><figcaption>Fig. 30 The Riemann fan. See text for the nomenclature of the waves.</figcaption></figure>

or all eigenvectors (and associated eigenvalues) that characterize the hyperbolic equations. In the case of MHD equations, the solution of the Riemann problem is controlled by seven waves, either discontinuities or rarefaction fans. Each wave is associated with one eigenvalue [characteristic velocities (Eq. 58)], where $\lambda_{1,3,5,7}$ are associated with shock or rarefaction waves, $\lambda_{2,6}$ are associated with rotational discontinuity and $\lambda_4$ is associated with a contact discontinuity (see Fig. 30).

There are different types of Riemann problem solvers (see Leveque 1998; Martí and Müller 2015):

- An exact Riemann solver requires an iterative method and thus is impractical for a MHD code (e.g., Takahashi et al. 2014; Torrilhon 2003). However, the results of the exact Riemann problem are used as reference solution to check the numerical precision and performance of approximate Riemann solvers.

- A linearized Riemann solver such as the Roe solver (Roe 1981) requires a decomposition of the left and right eigenvectors into characteristics. This solver is complex and time-consuming in MHD and HD cases. Moreover, the Roe solver can lead to negative solutions (Einfeldt et al. 1991).

- A guess-based Riemann solver such as the HLL solver (Harten–Lax–van Leer) introduces an estimate of the wave speed and the solution is averaged over the Riemann fan.

In this review, we focus on the most used Riemann solvers for MHD, all in the guess-based category. There are three main such solvers:

- The TVDLF (Total Variation Diminishing Lax–Friedrich Rusanov) Riemann solver (e.g., Bouchut 2004) is constructed by assuming a mean state which is given by the fastest wave. Indeed, the flux at cells interface $i + 1/2$ is given by,

$$F^{i+1/2} = \frac{1}{2} [F(U_{L,i+1/2}) + F(U_{R,i+1/2}) - S_{\max}(U_{R,i+1/2} - U_{L,i+1/2})], \quad (62)$$
where the fastest propagating speed at cell interface is:

$$S_{\max} = \max(S_{L,\max}, S_{R,\max}),$$

with $S_{L,\max}$, $S_{R,\max}$ are the eigenvalues associated with the fastest wave propagating respectively at the left and at the right of cell interface $i + 1/2$. They are defined as,

$$S_{L,\max} = \max_{(1 \le k \le 7)} (|\lambda_{k,L}|), \quad S_{R,\max} = \max_{(1 \le k \le 7)} (|\lambda_{k,R}|). \qquad (63)$$

- The HLL solver proposed by Harten et al. (1983) is constructed by assuming an average intermediate state between the fastest and slowest waves.

$$F^{i+1/2} = \begin{cases} F(U_{L,i+1/2}) & S_{L,i+1/2} > 0 \\ F(\text{HLL}, i + 1/2) & S_{L,i+1/2} \le 0 \le S_{R,i+1/2} \\ F(U_{R,i+1/2}) & S_{R,i+1/2} < 0 \end{cases} \qquad (64)$$

where

$$F(\text{HLL}, i + 1/2) = \frac{S_{R,i+1/2} U_{R,i+1/2} - S_{L,i+1/2} U_{L,i+1/2} + F(U_{L,i+1/2}) - F(U_{R,i+1/2})}{S_{R,i+1/2} - S_{L,i+1/2}} \qquad (65)$$

The outermost wave speed $S_{L,i+1/2}$ and $S_{R,i+1/2}$ are estimated using the left and right states,

$$S_L = \min_{(1 \le k \le 7)} (\lambda_{k,L}, \lambda_{k,R}) \\ S_R = \max_{(1 \le k \le 7)} (\lambda_{k,L}, \lambda_{k,R}) \qquad (66)$$

- The HLLC solver proposed by Toro et al. (1994) is a two-state HLL Riemann solver. It introduces sub-structures associated with a contact discontinuity into the sub-slow state of the HLL Riemann solver.

- The HLLD solver (Miyoshi and Kusano 2005) is a four-state HLL Riemann solver. The HLLD Riemann solver introduces sub-structures associated with the two rotational discontinuities $\lambda_{a,L}$, $\lambda_{a,R}$ separated by the contact discontinuity.

In the Riemann solvers presented above, the numerical estimation of the fluxes at interface $F^{i+1/2}$ requires the values of the conserved variables to the left $U_{L,i+1/2}$ and to the right $U_{R,i+1/2}$ of the cell interface located at $x_{i+1/2}$. These are reconstructed from cell centred values $U$. These spatial reconstructions can be performed by using a slope limited scheme to keep the reconstruction monotonic. There are various slope limited schemes such as minmod, super-bee and monotonized central difference limiter (MCD) (Toro 2009), ppm (Colella and Woodward 1984), KOREN (Koren and van der
Maarel 1993), van Leer (van Leer 1979). The minmod limiter is the most stable in the presence of strong discontinuities and is very efficient in decreasing numerical instabilities. MCD and Super-bee limiters are more efficient in the vicinity of smooth flows because they permit to retrieve a centered slope. PPM and KOREN are higher order limiters and thus with lower numerical dissipation, they can be used on large classes of problems with smooth flow, intermediate discontinuities and some strong shocks.

Numerical resolution of the partial differential Eq. (52) also requires high order temporal accuracy. This is realized by using a second order predictor-corrector scheme or a higher order scheme such as the strong stability preserving Runge-Kutta scheme.

### 5.2.3 Semi-implicit and implicit schemes

The various complex physical processes associated with astrophysical plasma phenomena act on widely different time scales. The precise treatment of these phenomena with an explicit MHD approach relies on the accuracy to which one can capture the dynamics as imposed by the CFL condition. In some cases, the resulting short time-resolution can make the MHD simulation computationally intractable with the explicit method. To make these simulations computationally feasible, it is necessary to integrate the MHD equations with larger time steps. This is possible when we are not interested in tracking all fast waves in the system and we can step over some unimportant propagating waves. This can be realized by the use of implicit schemes, which are preferable for complex astrophysical MHD simulations. However, parallel computation with implicit schemes becomes less efficient than with explicit schemes since implicit schemes use iterative algorithms that request more communication between the different processes.

In the implicit method the original system of differential Eq. (54) is rewritten as follows:

$$ \partial_t U = R(U), \qquad (67) $$

where $R(U) = -\partial_t \mathbf{F}^i(U) + S(U)$ and $S(U)$ represents source terms. $R(U)$ is a non-linear function of $U$. Eq. (67) is usually discretized in time using a third-level Backward Differentiation Formula method

$$ U^{n+1} = U^n + \Delta t \left[ \beta \mathcal{R}(U^{n+1}) + (1-\beta) \frac{U^n - U^{n+1}}{\Delta t^{n-1}} \right], \qquad (68) $$

where $\beta = (\Delta t^n + \Delta t^{n-1}) / (2 \Delta t^n + \Delta t^{n-1})$. In the case of constant time steps $\beta = 2/3$ and Eq. (68) is simplified to a second-order Backward Differentiation formula; with $\beta = 1$ Eq. (68) corresponds to a backward-Euler scheme. In the case $\beta = 0.5$ the scheme is fully implicit.

In the implicit method, at each time step, the multidimensional nonlinear system

$$ G(U^{n+1}) = U^{n+1} - \left\{ U^n \Delta t \left[ \beta \mathcal{R}(U^{n+1}) + (1-\beta) \frac{U^n - U^{n+1}}{\Delta t^{n-1}} \right] \right\} = 0, \qquad (69) $$

must be solved to determine the time-updated solution for the state $U^{n+1}$.
In some specific simple cases, such as 2D isothermal hydrodynamics, Eq. (69) can be preconditioned analytically to obtain a diagonal matrix, leading to a system of linear equations that can be solved by relaxation techniques (Harlow and Amsden 1968). The linearization (analytical preconditioning) of Eq. (69) can be done either about the initial equilibrium $G(U^n)$ or about the current state in order to construct the implicit operator $G(U^{n+1}) = 0$ needed to advance to the next time step. These classes of implicit schemes are used for specific problems and geometry (e.g., Harned and Schnack 1986). Another class of implicit schemes uses Newton-Krylov techniques for the resolution of nonlinear systems like Eq. (69). Newton's method consists of the local linearization of Eq. (69) for each state $U^n(m)$ at iteration $m$ according to

$$G(U^n(m)) = (U^n(m-1)) + \frac{\partial G(U)}{\partial U} (U^n(m) - U^n(m-1)) + \mathcal{O}(\Delta t^2), \quad (70)$$

which can be substituted into Eq. (69) giving the equation to be solved,

$$\left(I - \Delta t \beta \frac{\partial G(U)}{\partial U}\right) \delta U = \Delta t \left( \beta U^{n(m-1)} + (1 - \beta) \frac{U^{n(m-1)} - U^{n(m-2)}}{\Delta t^{n(m-2)}} \right), \quad (71)$$

and after the values of $U$ is updated as

$$U^{n+1} = U^n + \delta U. \quad (72)$$

However, this nonlinear resolution method benefits tremendously from accurate initial guesses, $U^{n(0)}$. In many cases an explicit predictor is used to provide an initial guess to the implicit scheme (e.g., Reynolds et al. 2006). The iterative resolution of the linear Eq. (71) is performed in general by using preconditioned Krylov (sub)space solvers (e.g., Tóth et al. 2006) with the non-restarted generalized minimum residual method (GMRES) iterative solver (Brown and Hindmarsh 1989). These linear solvers are very efficient for large-scale problems since they do not require storage of the matrix.

### 5.2.4 Magnetic divergence-free algorithms

The resolution of Euler and Maxwell equations using the standard Godunov schemes does not work by default in maintaining the divergence-free property of the magnetic field $\nabla \cdot \mathbf{B} = 0$. The resulting error accumulated during the simulation may grow to the point that it produces unphysical forces (Tóth 2000). Several strategies have been undertaken to handle the magnetic field evolution in numerical MHD. They are classified into two main categories.

We first find the divergence-cleaning schemes, where the evolution of magnetic field components is treated as any other MHD variables and only in a second step a divergence-cleaning procedure is applied. In these schemes the magnetic field components are defined at cells center as others variables. Then, the MHD equations are solved by adding source terms function of $\nabla \cdot \mathbf{B}$. There are various methods for divergence-cleaning schemes, such as the Generalized Lagrange Multiplier (GLM)
(Dedner et al. 2002). In the GLM method, a new transport variables $\Psi$ and its governing equation is introduced into the MHD equations system, which plays the role of advection and dissipation of the local divergence error. The divergence-cleaning can be treated as well by the eight-wave formulation approach Powell et al. (1999). There is also the projection method (Brackbill and Barnes 1980). In this scheme at each iteration the Poisson equation $\nabla^2\Phi = -\nabla \cdot \mathbf{B}$ is solved and in the second step the magnetic divergence part is removed from the magnetic field $\mathbf{B} = \mathbf{B} - \nabla\Phi$. Let us also mention the vector divergence-cleaning scheme (Balsara 1998) which is an extension of the projection method. Finally, another divergence-cleaning scheme based on the use of an artificial diffusivity added at each time step, following the completion of the TVD Lax-Friedrich scheme (van der Holst and Keppens 2007), where a term $\eta\nabla \cdot \mathbf{B}$ is added to the magnetic field to diffuse the magnetic divergence.

A second category of schemes is based on constrained transport methods, originally introduced by Evans and Hawley (1988). Such schemes use a staggered mesh formulation which is inherently divergence-free. In this method the magnetic field is defined at face centers and the remaining fluid variables are defined at cell centers. In this approach, the electric field is set along the cell edges. This method sustains a specified discretization of the magnetic field divergence around machine round off error.

The constrained transport method is attractive from a physical point of view, however, it requires specific treatment for magnetic field variables different from others variables, which is inconvenient for implementation specifically in the AMR. The diffusive method reduces the numerical error of $\nabla \cdot \mathbf{B}$ by adding a source term in the induction equation and the energy equation. The projection method involves an additional Poisson equation which significantly increases the computational cost. The GLM method is based on the use of central cell magnetic field and thus it can easily be applied on general grids.

### 5.2.5 Adaptive mesh refinement techniques

The numerical resolution of the PDEs (Eq. 57) uses a discrete domain. Therefore its precision depends on the mesh resolution (spacing), determined according to the scales of the phenomenon under study. In fluid mechanics a broad variety of spatial perturbations exist and can interact with each other. The complexity of these interactions requires the resolution of the problem at all scales. With uniform meshes, if high resolution is required throughout the computational domain the simulation can become computationally extremely costly.

Adaptive Mesh Refinement (AMR) addresses the problem of resolving this wide range of scales by increasing the spatial resolution around small scale structures. It is achieved by increasing locally the mesh resolution and then adjusting the computational effort locally to maintain a uniform level of accuracy throughout the computational domain. This type of AMR approach is called the *h*-type refinement. It consists in the splitting of existing elements into smaller ones. The development and use of the AMR starts with Berger and Oliger (1984), the transition from serial to parallel computing occurred after Griebel and Zumbusch (1999). Various AMR approaches exist, depending on the cells shape and the logical grouping of the cells

<figure><img src="image_52.png" /></figure>
on the mesh: gathering the cells according to their size $h$, inducing a particularly strict ansatz in the hierarchy, $h = n^{-L}$, corresponding to some refinement level $L$ (L is an integer). The cells volume at level $l$ are $h^{-ND}$ (ND = number of dimensions) smaller than the cells at coarse level $l = 1$. The most used hierarchy is $n = 2$ where the coarse cell is divided by $2^{ND}$ when it is refined.

In the block structured methods the cells are arranged in blocks according to their levels only. This method does not have any constraints on the size or shape of these blocks. Tree-base methods impose constraints on the block size and shape. With this method, the blocks are organized hierarchically as a quadtree in 2D (octree in 3D). In this distribution, the blocks in use at level $l$ represent the tree leaves since they have no children and they have an associated parent (ascendant) blocks. The parent blocks that are at the lowest refinement level $l = 1$ represent the tree roots. In some simulations, there is a need to use multiple tree roots arranged in an unstructured AMR, giving rise to a forest of trees.

### 5.2.6 Errors estimator

The AMR consists in the use of a coarse grid over the entire computational domain and refined grids only in some specific regions where the local truncation errors are judged to be too large to maintain a given numerical accuracy. These errors can be computed using the Richardson-estimator which compares the evolution of the variables at two successive levels. This estimator is accurate however it requires a lot of memory and is time consuming. Another estimator is the Löhner-estimator (Lohner 1987), a modified central second derivative normalized by the sum of first-order forward and backward gradients. It has the advantage of using mostly local calculations of any variables of the simulation and their combinations.

### 5.2.7 Load balance

In the simulation box, the distribution of blocks across the processors requires the use of a space filling curve. The most used space filling curves (SFCs) are Hilbert, Peano, and Sierpinski curves (Bader 2013). The Hilbert and Peano curves use recursive algorithms. The Morton order is also widely used due to its ease of implementation in the space filling. However the curves it generates are not continuous and thus do not fit into the family of finite SFCs.

With this curve the blocks are organized over a forest of trees and thus they are distributed over all the processors. There are two approaches for this forest partition. In the first approach, each tree and its leaves are associated to one owner process. In the second approach, the tree leaves can belong to multiple processes. The first approach is simple to implement however it does not provide a right load balance since the number of blocks per process may differ. The second approach, even if it presents issues with shared tree between processes, provides a perfect load balance, the difference in block distribution over processes is at most one. This last approach provides the best scalability (Keppens et al. 2012).
<figure><img src="image_53.png" /><figcaption>**Fig. 31** Sod shock tube test in MHD. The density, pressure, velocity components, transverse component of the magnetic field, and specific internal energy [scaled by $(\gamma_{ad} - 1)$] for the Brio & Wu (1988) shock tube problem are plotted at $t = 0.08$, computed with 400 grid points, second-order spatial reconstruction, and Roe fluxes. The solid line is a reference solution computed with 104 grid points. Image reproduced with permission from Stone et al. (2008), copyright by AAS</figcaption></figure>

### 5.2.8 Standard numerical tests

The development of numerical tools requests extended tests for all implemented physics and algorithms. Many comprehensive and well documented sets of tests have been presented in the literature. Tests are set for all physics, all dimensions, and implemented geometries and AMR schemes. For many 1D tests an exact solution exists and it is possible to compute the deviation of the numerical result from it using some error norm (e.g., norm error $L_1 = \sum_i ||q_i - q_i^0||/N$, where N is number of point, $q_i^0$ is the exact solution at cell i, $q_i^0$ is the numerical result at cell i).

One of the standard one-dimensional tests is the Sod MHD shock tube (Brio and Wu 1988). The shock tube test consists of two constant states, one on the left and the second on the right, separated by a discontinuity. This test allows to check the ability of a numerical code to treat correctly the Riemann problem in 1D and resolve the Riemann fan evolution, see Fig. 31.

For two-dimensional MHD tests, one of the standard tests in the vortex of Orszag & Tang (1979), see Fig. 32. This test consists in a doubly periodic fluid configuration leading to 2D supersonic MHD turbulence. The density and pressure are set to constant values, while the velocity and magnetic field are set as

<figure><img src="image_54.png" /></figure>
<figure><img src="image_55.png" /><figcaption>Fig. 32 Orszag–Tang vortex test. Contours of selected variables at $t_f = 1/2$, computed using a grid of 192 × 192 cells, third-order reconstruction, and Roe fluxes. Thirty equally spaced contours between the minimum and maximum are used for each plot. Image reproduced with permission from Stone et al. (2008), copyright by AAS</figcaption></figure>

$$ 
\begin{aligned}
\mathbf{u} &= (-\sin y, \sin(x), 0), \\
\mathbf{B} &= (-\sin y, \sin(2x), 0).
\end{aligned} 
\qquad (73) 
$$

This test does not have an analytical solution and the results have to be compared between different codes.

### 5.2.9 Relativistic magnetohydrodynamics

The MHD description introduced in the previous paragraphs is relevant for non-relativistic plasma velocities and energies. It is adequate for plasmas in the interstellar medium and in the vicinity of stars. However, some astrophysical phenomena involve relativistic plasma flows with energies of the order of the mass energy. In order to model these plasmas in the framework of the fluid model, the MHD formulation has to be revised using the special relativity (SR) framework. In flat space-time, the plasma is described by the following SR-MHD equations (see Martí and Müller 2015)
$$ \partial_t (\gamma \rho) + \nabla \cdot (\rho \Gamma \mathbf{u}) = 0, $$

$$ \partial_t (\gamma^2 h \mathbf{u} + \mathbf{E} \times \mathbf{B}) + \nabla \cdot (\gamma^2 h \mathbf{u} \mathbf{u} - \mathbf{E}\mathbf{E} - \mathbf{B}\mathbf{B} + p\mathbf{I}) = \mathbf{0}, $$

$$ \partial_t \left( h\gamma^2 - p - \gamma\rho + \frac{E^2 + B^2}{2} \right) + \nabla \cdot \left( \left( h\gamma^2 - p - \gamma\rho \right) \mathbf{u} + \mathbf{E} \times \mathbf{B} \right) = 0, $$

$$ \partial_t \mathbf{B} + \nabla \cdot (\mathbf{u}\mathbf{B} - \mathbf{B}\mathbf{u}) = \mathbf{0}, \qquad (74) $$

where the closure of this system of equations is provided by the equation of state expressed with the enthalpy $h = h(\rho, p)$. The total pressure is $p_t = p + \frac{E^2+B^2}{2}$ where the electric field is given by Ohm's law. In the ideal case $\mathbf{E} = -\mathbf{u}/c \times \mathbf{B}$, finally $\Gamma = (1 - (u/c)^2)^{-1/2}$ is the Lorentz factor of the flow.

The numerical resolution of SR-MHD equations exploits the same type of algorithms as presented before, by solving Eq. (74) in conservative form

$$ \partial_t U + \partial_t F^i(U) = 0, \qquad (75) $$

where the conserved variables are

$$ U = \begin{pmatrix} \gamma \rho \\ h\gamma^2 \mathbf{u} + \mathbf{E} \times \mathbf{B} \\ h\gamma^2 - p - \gamma \rho + \frac{E^2 + B^2}{2} \end{pmatrix}, \qquad (76) $$

and the corresponding flux tensor

$$ F(U) = \begin{pmatrix} \gamma \rho \mathbf{u} \\ h\gamma^2 \mathbf{u}\mathbf{u} - \mathbf{E}\mathbf{E} - \mathbf{B}\mathbf{B} + p\mathbf{I} \\ (h\gamma^2 - p - \gamma\rho)\mathbf{u} + \mathbf{E} \times \mathbf{B} \\ \mathbf{u}\mathbf{B} - \mathbf{B}\mathbf{u} \end{pmatrix}. \qquad (77) $$

Most of the complications in relativistic MHD comes from the non-linear relation between the primitive variables $(\rho, \mathbf{u}, p)$ and conserved variables $U$:

$$ \xi - p - \tau - D + B^2 - \frac{1}{2} \left( \left( \frac{B}{\gamma} \right)^2 + \left( \frac{\mathbf{S} \cdot \mathbf{B}}{\xi} \right)^2 \right) = 0, $$ $$ \mathbf{u} \cdot \mathbf{B} = \frac{\mathbf{S} \cdot \mathbf{B}}{\xi}, $$

$$ \frac{1}{\gamma^2} = 1 - \frac{(\mathbf{S} + (\mathbf{u} \cdot \mathbf{B})\mathbf{B})^2}{(\xi + B^2)^2}, $$ $$ \mathbf{u} = \frac{\mathbf{S} + (\mathbf{u} \cdot \mathbf{B})\mathbf{B}}{\xi + B^2}, \qquad (78) $$

where $D = \gamma\rho$ is laboratory frame density, $\xi = \gamma^2\rho h$, $\mathbf{S} = h\gamma^2\mathbf{u} + \mathbf{E} \times \mathbf{B}$ the momentum and $\tau = h\gamma^2 - p - \gamma\rho + \frac{E^2+B^2}{2}$. Eq. (78) can be handled only numerically,

<figure><img src="image_56.png" /></figure>
the most used scheme is the Newton–Raphson method (Noble et al. 2006). The iteration is performed on the pressure, on the enthalpy or the velocity. In SR hydrodynamics the iterative method can be avoided and instead a quartic equation can be solved (Schneider et al. 1993).

As in classical MHD, the SRMHD schemes exploits characteristic speed of plasma normal modes, but these are limited by the light speed, i.e.

$$-1 \le \lambda_1 \le \lambda_2 \le \lambda_3 \le \lambda_4 \le \lambda_5 \le \lambda_6 \le \lambda_7 \le 1. \qquad (79)$$

For a specific direction idim, the entropy wave as in classical MHD travels with speed $\lambda_4 = v^{idim}$, the Alfvén wave has speed

$$\lambda_{2,6} = u^{idim} \mp \frac{B^{idim}}{\sqrt{(\rho h + B^2) \mp (v \cdot B)}}, \qquad (80)$$

and the magneto-acoustic speeds are found from the quartic equation

$$\rho h (1 - c_s) \gamma^4 (\lambda - u^{idim})^4 - (1 - \lambda^2) \times \left( \gamma^2 (\rho h c_s^2 - B^2) (\lambda - u^{idim})^2 - c_s^2 \left( \Gamma (u \cdot B) (\lambda - u^{idim}) - \frac{B^{idim}}{\gamma} \right)^2 \right) = 0 \qquad (81)$$

where $c_s$ is the sound speed. In the hydrodynamics case, Eq. (81) becomes a single quadratic expression, and then the use of exact analytic formulae for the root evaluation is straightforward. In SRMHD Eq. (81) is quartic, although an analytical solution exists (Del Zanna et al. 2003), it is more easily obtained by numerical iteration, for which a Laguerre method can be used.

**Standard numerical tests** In special relativistic hydrodynamics and magnetohydrodynamics several numerical experiments are used as benchmark tests. Like in classical HD and MHD there is the Sod shock tube and the equivalent of the Orzag–Tang vortex (see above). Here we show the rotor test (Del Zanna et al. 2003), in Fig. 33. It consists in a disk of radius 0.1 with higher density, $\rho = 10$, positioned at the center of the computational domain $[0, 1] \times [0, 1]$, rotating at high relativistic speed, $\Omega = 9.95$, thus the Lorentz factor at the disk edge $\Gamma_{max} \simeq 10.0$, the rotor is embedded in a static background with $\rho = 1$, $p = 1$, and uniform magnetic field $B_x = 1$ with polytropic index $\gamma = 5/3$.

## 5.3 Multi-fluid methods

The previous section presented the case of a single fluid. However, in astrophysics the thermal plasma pressure is usually in competition with the pressure in other components: cosmic rays, radiation fields, dust, neutral species, or requires to include an
<figure><img src="image_57.png" /><figcaption>**Fig. 33** The relativistic rotor problem. The top panel shows the density structure (both in linear scale as well as using a Schlieren plot), the bottom panel the magnetic field components, at time $t = 0.4$. Image reproduced with permission from Keppens et al. (2012), copyright by Elsevier</figcaption></figure>

electron fluid. In bi- or multi-fluid models it is in principle necessary to add as many fluid equations as the number of species to be included in the simulations, and to account for the friction forces induced by collisions among the species in the momentum conservation equation (see, e.g., O'Sullivan and Downes 2007).

As in this review we are mostly interested in non-thermal particle acceleration and transport we restrict our discussion to the case of the interaction of a thermal fluid and a non-thermal (or cosmic ray) component (Drury and Voelk 1981; Rasera and Chandran 2008).

### 5.3.1 Model equations

CRs impact fluid dynamics through the effect of the gradient of their pressure $p_{CR}$ (Eq. 9). The momentum equation is modified as

<figure><img src="image_58.png" /></figure>
$$ \partial_t (\rho \mathbf{u}) + \nabla \cdot \left( \rho \mathbf{u}\mathbf{u} - \mathbf{B}\mathbf{B} + \left( p + p_{\text{CR}} + \frac{B^2}{2} \right) \mathbf{I} \right) = 0, \quad (82) $$

and the energy equation now includes the work of this force plus another component produced by the CR diffusion

$$ \partial_t (e) + \nabla \cdot \left( \left( e + p + p_{\text{CR}} + \frac{B^2}{2} \right) \mathbf{u} - (\mathbf{B} \cdot \mathbf{u}) \mathbf{B} + \mathbf{F}_{\text{CR}} \right) = 0, \quad (83) $$

where the CR flux is $\mathbf{F}_{\text{CR}} = -\bar{\mathbf{D}} \cdot \nabla e_{\text{CR}}$ and $e_{\text{CR}}$ is the CR energy density. The CR diffusion coefficient can be decomposed into parallel and perpendicular components with respect to the background magnetic field, namely

$$ D_{ij} = D_{\perp} \delta_{ij} + (D_{\parallel} - D_{\perp}) b_i b_j \quad (84) $$

where $\mathbf{b} = \mathbf{B}/B$ (Hanasz and Lesch 2003). An energy equation for the CR energy density $e_{\text{CR}}$ is now required. It reads

$$ \partial_t (e_{\text{CR}}) + \nabla \cdot (e_{\text{CR}} \mathbf{u} + \mathbf{F}_{\text{CR}}) = -p_{\text{CR}} \nabla \cdot \mathbf{u}. \quad (85) $$

CR pressure and energy density are linked by $p_{\text{CR}} = (\gamma_{\text{CR}} - 1)e_{\text{CR}}$, where a CR gas adiabatic index $\gamma_{\text{CR}}$ is introduced. The above energy equation does not have a flux-conservative form. Kudoh and Hanawa (2016) (see also Pfrommer et al. 2006) propose an alternative approach leading to a full set of flux-conservative equations for the CR-MHD system. To proceed the authors introduce the CR mass density $\rho_{\text{CR}} = p_{\text{CR}}^{1/\gamma_{\text{CR}}}$, thus approximating CRs as a polytropic gas. With this assumption the CR energy equation can be recast into a continuity equation for the CR gas

$$ \partial_t (\rho_{\text{CR}}) + \nabla \cdot (\rho_{\text{CR}} \mathbf{u}) = 0. \quad (86) $$

The previous equation implies that that $\rho_{\text{CR}}/\rho$ is conserved along a streamline.

### 5.3.2 Specific numerical schemes of CR-fluid systems

The new set of equations can be solved using FVM as for standard MHD equations. The solver now has to account for the CR pressure which modifies the local sound speed $c_s = \sqrt{\gamma_{\text{ad}} P/\rho + \gamma_{\text{CR}} P_{\text{CR}}/\rho}$. The CR pressure is dominated by the relativistic part of the CR distribution hence $\gamma_{\text{CR}} \simeq 4/3$. The standard solvers detailed in Sect. 5.2.2 can be used to treat the above CR-MHD system. Working in the framework of the fluid-kinetic approach (Sect. 5.4), Miniati (2007) develop a modified Glimm-Godunov solver where the CR mediation is included in the Riemann problem. Kudoh and Hanawa (2016) propose a CR+MHD solver, second order accurate in space and time for a bi-fluid system based on a Roe solver.
**Semi-implicit and implicit methods** An important difficulty of explicit schemes for CR-HD or CR-MHD systems comes from the CFL stability criterion for the CR diffusion, which imposes a timestep $\Delta t \le X_{CFL} \Delta x^2/(2D)$, where $X_{CFL}$ is the CFL number. This criterion, because the time step scales non-linearly with the grid resolution, imposes severe slowing down limitations, especially in multi-scale problems where AMR is active. Several semi-implicit or implicit approaches have been proposed to cure this issue, usually in the context of thermal conduction studies (Balsara et al. 2008). The CFL constraint from CR diffusion is alleviated by calculating the diffusion operator using an implicit method (see Sect. 5.2.3 and next). Another way to reduce the computation time is the so-called super-time stepping technique (O'Sullivan and Downes 2007; Balsara et al. 2008) where the CFL condition is imposed over a large time interval composed of multiple elementary substeps over which the stability condition can be relaxed. To increment the guessed solutions at every substep from the previous guess a Runge-Kutta method is applied, using a polynomial recursion relation (using either Chebyshev or Legendre polynomials).

Here we discuss more specifically the implicit scheme proposed by Dubois and Commerçon (2016) which is well-adapted to simulations with AMR. The diffusion operator in the energy Eq (83) can be discretized as (we write it here in 2D for a Cartesian grid but it can easily be generalized to 3D)

$$e_{i,j}^{n+1} = e_{i,j}^{n} - \frac{\Delta t}{\Delta x} \left( F_{i+1/2,j}^{n+1} + F_{i,j+1/2}^{n+1} - F_{i-1/2,j}^{n+1} - F_{i,j-1/2}^{n+1} \right) \quad (87)$$

where the energy density is calculated at the cell center (i,j) and the fluxes are obtained at cell interfaces. The quantities are expressed at the final time step $t^{n+1}$, forming a linear system that can be solved by matrix inversion. In this scheme the anisotropic part of the fluxes at the cell interfaces are calculated from the fluxes at the cell corners, i.e., $F_{i+1/2,j}^{n+1} = 0.5(F_{i+1/2,j+1/2}^{n+1} + F_{i+1/2,j-1/2}^{n+1})$. The anisotropic part of the flux corresponds to the diffusion along the background magnetic field $D_{\|}bb \cdot \nabla e$.

**Cosmic-ray streaming regularization** CRs are not exactly advected with the background fluid but with the scattering centers (MHD waves) carried by the background fluid (Skilling 1975a). Considering slab-type waves that CRs can self-generate as they stream along the background magnetic field, we find a streaming velocity

$$\mathbf{u}_{\text{st}} = \mathbf{u} + \left\langle \frac{3}{2}(1 - \mu^2) \frac{\nu_{+} - \nu_{-}}{\nu_{+} + \nu_{-}} \right\rangle \mathbf{u}_{\text{A}} \quad (88)$$

where $\mathbf{u}$ and $\mathbf{u}_A$ are resp. the fluid and waves speed, $\mu$ is the CR pitch-angle cosine, and $\nu_{\pm}$ are the angular scattering frequency produced by the forward/backward (+/-) propagating waves along the background magnetic field. This calculation assumes that the quasi-linear theory of CR transport applies (Schlickeiser 2002). A detailed calculation of this velocity requires to know the scattering frequencies, which is only possible by adding two other energy equations for each type of propagating wave. This is the purpose of the next paragraph. However, self-generated waves are preferentially produced when the local CR pressure is in excess with respect to the gas and magnetic
pressure (as is likely the case close to CR sources), along a CR gradient. In that case, waves are preferentially triggered in one direction and one can write $\mathbf{u}_{\text{st}} = \mathbf{u} - \mathbf{u}_a\mathbf{B} \cdot \nabla p_{\text{CR}} / |\mathbf{B} \cdot \nabla p_{\text{CR}}|$ (Pfrommer et al. 2017). The CR energy Eq. (85) needs to be modified according to

$$ \partial_t (e_{\text{CR}}) + \nabla \cdot (e_{\text{CR}} \mathbf{u}_{\text{st}} + \mathbf{F}_{\text{CR}}) = -p_{\text{CR}} \nabla \cdot \mathbf{u}_{\text{st}}. \quad (89) $$

This equation is non-linear as $\mathbf{u}_{\text{st}}$ depends on $\nabla e_{\text{cr}}$. The LHS of Eq. (89) resembles an advection equation but with a speed dependent on the sign of the gradient of the CR energy density, which introduces some spurious oscillations at extrema where the gradient changes its sign (Sharma et al. 2010). These authors propose a regularization of the energy equation by replacing the sign of $\mathbf{B} \cdot \nabla p_{\text{CR}}$ in Eq. (89) with a smooth function. The drift speed is rewritten as

$$ u_{\text{drift}} = u_A \tanh \left( \frac{X}{\epsilon} \right), \quad (90) $$

where $X = \mathbf{B} \cdot \nabla p_{\text{CR}} / |\mathbf{B} \cdot \nabla p_{\text{CR}}|$, and $\epsilon$ is a small parameter to be adjusted. The energy equation becomes diffusive at streaming speed extrema with a diffusion coefficient dependent on $\epsilon$. Sharma et al. (2010) show that an implicit non-linear integration scheme (see their Eq. 3.9) produces a rapid convergence, but requires a sufficiently small time step to be adjusted with the value of $\epsilon$. Jiang and Oh (2018) propose an alternative approach where Eq. (89) is replaced by a system of two equations

$$ \begin{aligned} \partial_t (e_{\text{CR}}) + \nabla \cdot (\mathbf{\Phi}_{\text{CR}}) &= -p_{\text{CR}} \nabla \cdot \mathbf{u}_{\text{st}} \\ \frac{1}{u_m^2} \partial_t (\mathbf{\Phi}_{\text{CR}}) + \nabla P_{\text{CR}} &= -\frac{1}{D} \mathbf{\Phi}_{\text{CR}}. \end{aligned} \quad (91) $$

where $\mathbf{\Phi}_{\text{CR}} = \mathbf{F}_{\text{CR}} + e_{\text{CR}}\mathbf{u}_{\text{st}}$ is the total CR flux and $U_m$ is a speed in practice taken larger than the maximum natural mode speed of the plasma. In this approach the problematic term sgn($X$) is replaced by $\nabla \cdot \mathbf{\Phi}_{\text{CR}}$. The system is closed by choosing the diffusion coefficient as $D = D_0 - \mathbf{u}_{\text{st}}(e_{\text{CR}} + p_{\text{CR}}) / (\mathbf{b} \cdot \nabla p_{\text{CR}})$, where $D_0$ is a background coefficient produced by large-scale injected turbulence. Thomas and Pfrommer (2019) criticize the two previous methods. They present CR transport mediated by self-generated waves including an accurate description of CR pitch-angle scattering up to the second order in accuracy in $(u_a/c)$. The latter is possible because of the addition of two supplementary energy equations respectively for forward and backward propagating waves (see the next paragraph). The authors employed a system of equations similar to Eq. (91) but with supplementary terms derived from the effect of CR scattering off self-generated waves.

**A four-fluid approach** A more complete description of the CR-fluid system includes the description of the forward and backward propagating CR self-generated waves (Skilling 1975a). It includes two supplementary fluid energy equations (Ko 1992). Now the total pressure includes the contribution of forward (backward) waves $P_{w,+}$ ($P_{w,-}$), and the total energy density includes wave terms as well $e_{w,+}$ ($e_{w,-}$), see Ko
(1992), Thomas and Pfrommer (2019) for a complete derivation of the new system of equations. This system is restricted to the quasi-linear theory framework. Hence, coupled with a MHD code the amplitude of self-generated waves has to be small with respect to the amplitude of the background magnetic field. The four-fluid system has been numerically solved by Thomas and Pfrommer (2019) using an explicit scheme separating the CR and wave fluid equations.

The interest of this approach resides in the more accurate calculation of the CR streaming speed given by Eq. (88) as well as the CR diffusion coefficient in space and energy, both dependent on the amplitude of scattering frequencies $v_\pm$ off forward and backward waves. In particular, it is possible to account for a fluid description of a stochastic acceleration term in the CR fluid energy density equation. This term has the form (Thomas and Pfrommer 2019)

$$4 \frac{\bar{v}_+ + \bar{v}_-}{\bar{v}_+ + \bar{v}_-} \frac{v_a^2}{c^2} (e_{CR} + P_{CR}), \qquad (92)$$

where $\bar{v}$ is the momentum-averaged scattering frequency.

**Other numerical strategies** Pfrommer et al. (2006) and Enßlin et al. (2007) develop a method based on a smooth particle hydrodynamics (SPH) approach. In SPH fluid dynamics is treated using Lagrangian particles. To each particle is attached relevant fluid properties (e.g. density, pressure ...) calculated using a SPH kernel dependent on a smoothing length (Monaghan 1992). The CR distribution follows a power-law distribution whose normalization and power-law index vary under the effect of adiabatic gas variation (compression/expansion) as a function of the gas density (itself calculated in the SPH code).

### 5.3.3 Numerical tests

Aside the standard numerical tests for HD and MHD codes described in Sect. 5.2.8, we detail here some specific setups aiming at testing the transport (either passive or active) of CRs.

Hanasz and Lesch (2003) propose a series of tests of the CR flux term in Eq. (83) using the PIERNIK code (Hanasz et al. 2010). A first, straightforward test is 1D diffusion of CRs along the background magnetic field directed in one direction of the Cartesian grid, setting the perpendicular diffusion to zero, and no CR backreaction (turning off the CR pressure gradient), in a static medium. Diffusion can also be tested along an inclined magnetic field with contributions of the different directions to the CR flux. Figure 34 shows the profile of the CR energy density along and perpendicular to the ellipsoid solution of the propagation of an initial spheroidal distribution in 3D (see Fig. 2 in Hanasz and Lesch 2003).

The authors propose also the same test but now turning on the effect of CR pressure gradient. Figure 35 shows two plots: CR energy density and magnetic field lines (left), gas density and velocity field lines (right). It can be seen that the gas accelerates up to a few km/s under the effect of a strong CR gradient preferentially along the background magnetic field line. This gas motion leads to a drift of CRs along the magnetic field.
<figure><img src="image_59.png" /><figcaption>**Fig. 34** CR propagation along an inclined magnetic field. The two curves show cuts of the ellipsoid solution at a given time along the major axis (crosses) and along the minor axis (asterisks). The solid and dotted lines represent the analytical solutions corresponding cuts of the fitted 2D Gaussian profile. CR diffusion is treated using an explicit scheme with $D_\parallel = 100$ (in units of pc$^2$Myr$^{-1}$ and $D_\perp = 0$. Image reproduced from Hanasz and Lesch (2003), copyright by ESO</figcaption></figure>

<figure><img src="image_60.png" /><figcaption>**Fig. 35** CR propagation along an inclined magnetic field in the case CR backreaction is active. Equipartition between gas, magnetic field and background CR pressures is assumed. An over pressure of a factor 100 in CR is injected at the center at the start of the simulation. The diffusion coefficients are: $D_\parallel = 100$ and $D_\perp = 4$ (in units of pc$^2$ Myr$^{-1}$). Left: CR pressure map and magnetic field lines. Right: gas density map and velocity vectors. Image reproduced with permission from Hanasz and Lesch (2003), copyright by ESO</figcaption></figure>

Perpendicular diffusion imposes a broadening of the CR profile perpendicular to the magnetic field lines. Similar tests are proposed in Snodin et al. (2006) (see their Fig. 5).

A second type of tests used in the context of CR acceleration at SNR blast waves involves Sod shock-tube simulations (restricted to HD). Pfrommer et al. (2017) (see also Salem and Bryan 2014) derive an analytical solution of the shock-tube problem including a CR gas. Figure 36 shows a solution of a 1D Riemann shock tube problem. Three cases are shown: on the left the solutions for a shock propagating in a composite gas of thermal plasma and CRs but without any CR acceleration, in the middle the same case but now including CR acceleration, on the right a shock propagating in the thermal gas only but with CR acceleration. Butsky and Quinn (2018) propose a 1D Brio–Wu shock-tube test (hence in MHD).
<figure><img src="image_61.png" /></figure>

**Fig. 36** 1D shock-tube test problems including CRs. From top to bottom are shown: the gas mass density, gas pressure, gas velocity and shock Mach number. Left: shock propagation in a composite gas including CRs but without the effect of CR acceleration. Middle: same as the left case but now including CR acceleration. Right: same as the middle case but without background CRs. Solid lines: analytical solutions, dotted lines: numerical solutions. The different pressures are: gas (blue), CRs (orange), total (red), the pressure in the CR injected at the shock front is in green. In the left and middle cases the initial CR pressure is two times the gas pressure on the left part and equal to it on the right part. Image reproduced with permission from Pfrommer et al. (2017), copyright by the authors

### 5.3.4 Cosmic ray physics using a bi-fluid approach

We illustrate here the ability and the limits of the bi-fluid approach in treating three important aspects of CR physics: CR acceleration at shocks, CR-driven winds and CR-induced magnetic production in the Galaxy.

**Cosmic ray acceleration at shocks** After the derivation of the test-particle solution produced by DSA (see Sect. 2.2), it appeared that the pressure carried by the high energy CRs is large enough to modify the shock structure and hence the Fermi acceleration process itself. Drury and Voelk (1981), Axford et al. (1977) derived self-consistent fluid solutions including the CR pressure in the Rankine–Hugoniot conditions using a bi-fluid model. They find that the non-linear model has up to three stationary shock solutions when the adiabatic indices of the gas and CR are $\gamma_{ad} = 5/3$ and $\gamma_c = 4/3$ respectively. The selection among these solution depends on the injection efficiency (Malkov 1997b). Saito et al. (2013)$^{46}$ discuss the time-dependent stability of these solutions and confirm that the high- and low-efficiency solutions are stable against large-amplitude perturbations, while the intermediate one is not. At high-injection

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

<footnote>$^{46}$ Section 2 of Saito et al. (2013) provides a short updated review of the main assumptions and issues of the bi-fluid model.</footnote>
efficiency and high-shock Mach numbers the shock may be completely smoothed in the two-fluid approach. Malkov and Drury (2001), Malkov (1997a) show that a sub-shock must always exist in a kinetic model and hence that these smooth solutions are an artefact of the model because the maximum CR energy is always finite.

Acknowledging these issues, the main interest of the bi-fluid approach to investigate CR back-reaction resides in its simplicity and also in the fact that this formalism can be easily combined with a MHD model (Webb et al. 1986; Jun et al. 1994; Frank et al. 1995) and used to follow CR-MHD fluid non-linear dynamics with respect to the shock obliquity. The bi-fluid model also allows to combine effects of CRs and radiation (Wagner et al. 2007). It is particularly useful to diagnose the profile of $H_\alpha$ lines in Balmer-dominated shocks occurring while a SNR propagate in a partially ionized medium (Wagner et al. 2009). However, only a kinetic calculation permits to follow the energy dependence of the particle distribution and the instabilities particles may generate. Non-linear kinetic models have since superseded the bi-fluid approach in this context. Numerical non-linear kinetic models are presented in Sects. 5.4 and 5.5 and semi-analytic models in Sect. 5.7.

**CR driven winds** CRs by the interplay of the generation of waves for instance by the streaming instability can convert a part of the CR bulk momentum into fluid momentum and hence drive winds (Ipavich 1975). Breitschwerdt et al. (1991), Breitschwerdt et al. (1993) investigate both analytically and numerically the launch of CR-driven winds in the framework of the bi-fluid model. Actually as their model includes also self-generated waves it proposes a three-fluid description. The authors use the flux tube approximation where stationary MHD equations are solved along the outwardly-directed magnetic field lines. The flux tube has a variable area cross section $A(z)$, $z$ marking the height above the galactic disk. Solutions, especially the sonic point of the MHD system, are searched for as in the case of the solar wind. The numerical solutions start from the calculation of the fluid speed gradient at the critical point and then are propagated out- and downward to match the inner and outer boundaries. Stationary CR-driven outflows are obtained if the inter galactic gas pressure is low enough. Zirakashvili et al. (1996) extend this work by including galactic rotation. They show that the wind forces the gas in the halo to corotate with the galactic disk up to a distance of a few kpc. The waves generated by the CRs also contribute to heat the halo. Everett et al. (2008) adopted a similar modelling as in previous works although with different boundary condition assumptions. They find that thermal and CR pressures are equally important to drive the wind. Their solutions are found to be consistent with an injection of CRs by SNe at a standard rate although with a bit high efficiency. Uhlig et al. (2012) carry SPH simulations including a treatment of adiabatic gain/loss of the CR gas and the physics of streaming. Pakmor et al. (2016), Ruszkowski et al. (2017) add the effect of anisotropic diffusion coefficients in the disk in a bi-fluid formalism. All these works show that CRs have a systematic negative feedback impact over the star formation rate in our Galaxy while they can launch powerful winds from the galactic disk.

**Galactic magnetic dynamo** Dynamo is likely at the origin of large scale magnetic fields generated in the Galaxy (see Kulsrud 1999 and Sect. 6 in Brandenburg 2018).
Bi-fluid models have been widely used to investigate the Parker or buoyancy instability thought to participate fast dynamo processes in spiral galaxies (Parker 1992; Hanasz and Lesch 1993). CRs produce an inflation of magnetic loops anchored in the galactic disk that induces magnetic reconnection and disconnects disk and halo magnetic fields. The closed loops in the halo and the anchored loops in the disk can be subject to cyclonic rotation at the origin of an $\alpha\Omega$ dynamo process, the $\alpha$ effect resulting from the Parker instability. Notice that CRs can provoke $\alpha$ dynamo because of their current, strong enough to trigger the non-resonant streaming instability (Beresnyak and Li 2014).

A numerical investigation of the Parker instability using a bi-fluid model has been proposed by Hanasz and Lesch (2000) using a flux-tube approximation. In that approximation CR are injected from a SNR at its lifetime end in the flux tube composed of the magnetic field lines threading the SNR. CRs are injected with an over-pressure of ~30 with respect to their background pressure and then inflate the flux tube in the vertical direction (with respect to the disk), eventually leading to a flux-tube explosion in a runaway process and ultimately to the production of galactic winds.<fnref n="47" /> The authors were able to find $\alpha$ coefficients large enough to ensure efficient Parker instability including magnetic field back-reaction. Ryu et al. (2003) investigate the growth rate of the Parker instability as function of the CR diffusion coefficient and find that the growth rate decreases with the parallel coefficient whereas the perpendicular coefficient has no strong impact. Kuwabara et al. (2004) conduct a 2D analysis of the Parker and explosion instabilities. Later Hanasz et al. (2004) and Hanasz et al. (2009) propose 3D simulations of the Parker instability using the PIERNIK code in a 3D box and in global Galactic geometry respectively. The global simulations show the growth of the magnetic field strength saturates at about $t = 4$ Gyr, reaching values of 3–5 $\mu$G in the disk. The magnetic field, initially randomly oriented, shows at the end of the simulation a toroidal component forming a spiral structure with reversals in the plane of the disk.

## 5.4 The fluid-kinetic framework

Solving the Fokker-Planck equation (FPE) is usually a very demanding task for computer simulations (see Sect. 3.6). Then, coupling a FPE with a MHD code has only been scarcely attempted yet. In this section we discuss approaches developed to couple a fluid model (HD or MHD) and a kinetic model (Vlasov or Fokker-Planck or approximations of these), appropriate to describe respectively the shock and the particles in the context of CR acceleration at shocks. In Sect. 5.4.1 we discuss methods that couple (M)HD equations and the diffusion-convection equation. In Sect. 5.4.2 we present a different approach that relies on a more general Vlasov-Fokker-Planck equation for description of the evolution of the distribution function $f(p)$.

<footnote n="47">⁴⁷ Figure 1 in Lesch and Hanasz (2003) lists a series of physical processes which result from the Parker instability.</footnote>
### 5.4.1 Coupling hydro and diffusion-convection equations

Works reviewed in this section model separately the evolution of the astrophysical flow (e.g., the blast wave of a supernova remnant, or a cosmological structure formation shock), using HD or MHD equations (see Sect. 5.1), and the acceleration of energetic particles, using the diffusion-convection equation, that is the angle-averaged Fokker-Planck equation (see Sect. 3.6) of the form

$$ \frac{\partial f}{\partial t} + \frac{\partial}{\partial x} (uf) = \frac{\partial}{\partial x} \left( D \frac{\partial f}{\partial x} \right) + \frac{1}{3p^2} \frac{\partial p^3}{\partial p} \frac{\partial u}{\partial x}, \quad (93) $$

which includes advection terms in space and in momentum, and a diffusion term in space (we are here restricting ourselves to first order acceleration, and so have not included diffusion in momentum). When written in a conservative form the spatial advection term may be included in the hydrodynamical solver, using the operator-splitting technique. The right-hand side is commonly solved with Finite Differences methods (FDM), using the techniques presented in Sect. 3.6. The diffusion term is the most difficult to treat. As stressed before, explicit schemes suffer from the constraint due to a CFL condition that is quadratic in the resolution $\Delta t < \Delta x^2/2D$, and which for physically-motivated values of the diffusion coefficient imposes time-steps that are much shorter than the hydrodynamic time-step. As in Sect. 5.3, this leads to the use of implicit or semi-implicit schemes, or possibly accelerated explicit schemes. The diffusion coefficient $D$, although it should in principle be computed from $f$ itself, is commonly prescribed as a function of $p$, usually as a power-law $D(p) \propto p^{\alpha}$, with optionally a dependence on position $x$ (the Bohm coefficient reduces to such a form with $\alpha = 2$ in the non-relativistic regime and $\alpha = 1$ in the relativistic regime). Compared to the bi-fluid approach of previous Sect. 5.3, the methods presented in this section aim at describing the spectrum of particles $f(p)$, using various approximations and techniques to handle the large range of scales to be resolved.

**The piecewise power-law method** This method was introduced by Jones et al. (1999), Tregillis et al. (2001) to investigate the time evolution of relativistic synchrotron-emitting electrons in radio galaxies.<fnref n="48" /> The basic idea is to approximate the electron distribution $f(p)$ as a piecewise power-law with $N$ bins in momentum, $f(p) = \sum_i f_{0,i}(p/p_i)^{\alpha_i}$ with $i \in [1, N]$ (on a fixed, logarithmic momentum grid), and to calculate $f_{0,i}$, $\alpha_i$ as function of the position in the simulation space (here a jet). The number of relativistic electrons in the bin $i$ is $n_i = \int_{p_i}^{p_{i+1}} f(p)4\pi p^2 dp$, which can be normalized to the grid mass density of the gas. We note $b_i = n_i/\rho$. At shocks, the standard DSA theory in the test-particle limit is applied (see Sect. 2.2.1), the index is fixed to $\alpha = 3r/(r-1)$ where $r$ is the shock compression ratio evaluated

<footnote n="48">⁴⁸ An alternative approach was proposed by Micono et al. (1999) where a set of Lagrangian test particles, a set of relativistic electrons is followed along a jet including acceleration and losses. Note that simulations combining mixed Lagrangian (for energetic particles) and Eulerian (for the magnetized fluid) coordinates have also been developed in the framework of relativistic flows by Vaidya et al. (2018). The authors solve the diffusion-convection equation for relativistic electrons transported in a relativistic flow calculated using a relativistic version of the PLUTO code.</footnote>
from the shock Mach number. At shocks $b = \sum_i b_i$ is calculated solving a simple equation $db/dt = Q_{inj}/\rho$ where $Q_{inj} = \epsilon u_{sh,L}/\mu_e m_p$. Here $u_{sh,L}$ is the Lagrangian shock speed, $\mu_e$ the electron mean mass, $m_p$ the proton mass. The parameter $\epsilon$ is the fraction of thermal electrons injected as non-thermal particles. Away from shocks, the parameters $b_i$ and $\alpha_i = 3 - (b_{i+1}/b_i)/\Delta \ln(p)_i$ evolve accounting for adiabatic and radiative losses (see Miniati 2001).

Adiabatic losses, which produce a shift in energy to lower energy, are of particular importance for CR-MHD systems discussed in Sect. 5.3. These systems can be extended to account for more than one CR population. This extension is reflected in the choice of the mean diffusion coefficient corresponding to a particular CR energy range of a given population. This effect can be treated with the piecewise power-law method (Miniati 2001; Girichidis et al. 2014). The adiabatic loss rate is $d \ln(p)/dt = (\gamma_{CR} - 1) \nabla \cdot \mathbf{u}$. Here the adiabatic CR index $\gamma_{CR}$ changes depending on the CR energy, from 5/3 in the non-relativistic case to 4/3 in the relativistic case. The CR density can be computed from the relation $n_i = E_i/\langle E_i \rangle$, where $E_i$ is the CR energy in the ith bin, and the average energy per CR particle in bin i $\langle E_i \rangle = \int_{p_i}^{p_{i+1}} E(p) p^2 dp / \int_{p_i}^{p_{i+1}} p^2 dp = e_{CR,i}/n_i$ does not depend on f. The CR energy density in the bin i can then be updated accounting for the loss term

$$e_{\text{CR},i}^{n+1} = e_{\text{CR},i}^{n} + \Delta t (\Phi_{i-1/2} - \Phi_{i+1/2}) \quad (94)$$

where the energy flux is $\Phi_{i+1/2} = 4\pi/\Delta t \int_t^{t+\Delta t} b(E) E p^2 f(p)_{i+1/2}$. The final energy $E(t + \Delta t)$ is expressed in terms of the adiabatic loss term and $\Delta t$ and has to be in the interval $[E_{i-1/2}, E_{i+1/2}]$, otherwise a sub-cycling is required (see Girichidis et al. (2014) for further details).

**Including backreaction effects** In efficient DSA the evolution of the shock and of the particles are non-linearly coupled (NLDSA). In the hydro-kinetic approach the connection between the fluid and the particles is made via extra terms in the equations, that model the injection and the back-reaction of particles. First, a source term is required in the diffusion-convection equation at the shock front, that depends on the shock jump conditions, and a corresponding loss term must be added in the set of hydrodynamic equations to ensure energy conservation (and in principle conservation of mass, but the inertia of energetic particles is normally negligible). Secondly, similar to the bi-fluid approach, one has to take into account the pressure of particles $P_{CR}$ given by Eq. (9). The force $-\nabla P_{CR}$ exerted by the particles on the flow is added as a source term to the equation of conservation of momentum, and the corresponding work $-\mathbf{u} \cdot \nabla P_{CR}$ is added as a source term to the equation of conservation of energy. At this level of modeling it is the gradient of particle pressure upstream of the shock that will cause the appearance of the shock precursor,<fnref n="49" /> which in turn will produce concave particle spectra.

The first investigations of time-dependent NLDSA using this approach were performed by Falle and Giddings (1987) and Bell (1987). Knowing that the canonical result for $f(p)$ is a power-law of index $s = 4$, Falle and Giddings (1987) work with

<footnote n="49">⁴⁹ At the PIC/hybrid level of modelling of Sect. 3.1, this is done via the Lorentz force.</footnote>
the quantity $g(p) = p^4 f(p)$, and since the momentum gain is proportional to the current momentum they replace $p$ with the quantity $y = \ln p$. Equation (93) is then rewritten as

$$ \frac{\partial g}{\partial t} + \frac{\partial (ug)}{\partial x} = \frac{\partial}{\partial x} \left( D \frac{\partial g}{\partial x} \right) + \frac{1}{3} \left( \frac{\partial g}{\partial y} - g \right) \frac{\partial u}{\partial x}. \quad (95) $$

These early works assumed a small dependence of the diffusion coefficient on $p$, using resp. $\alpha = 1/4$ and $\alpha = 1/2$, because of numerical limitations. As already noted in Sect. 2.2, it is the diffusion coefficient that sets the spatial and temporal scales of the simulation. The spatial scales range from the microscopic scale where the particles decouple from the fluid (of the order of a few thermal gyration lengths) to macroscopic scales (like the radius of the supernova remnant). The resolution of the numerical grid is dictated by the diffusion of the lowest energy particles, whereas the size of the grid is dictated by the diffusion of the highest energy particles. The ratio $D(p_{max})/D(p_{min})$, and thus the number of cells, may exceed ten orders of magnitude if $D(p) \propto p$, which is too demanding in terms of memory requirements and computing time. Fortunately we need very high resolution only around the shock, since low energy particles cannot travel far away from the shock. More generally, particles of a given momentum $p$ require a certain spatial resolution over a certain extent around the shock. This leads to the implementation of adaptive mesh refinement (AMR), reviewed in Sect. 5.2.5. In the context of DSA simulations, this technique was pioneered by Duffy (1992), to be able to use the true Bohm scaling $D(p) \propto pv$, and later used by Kang et al. (2001) and by Ferrand et al. (2008). The latter authors also parallelized the scheme in the $p$-direction, and studied repeated acceleration by multiple shocks. Ferrand et al. (2014a) use a different diffusion coefficient, proposed for perpendicular shocks, while DSA simulations are usually focused on the case of parallel shocks.

The hydro-kinetic approach was mostly developed by the team of Kang and Jones and collaborators with their code CRASH, starting from Kang and Jones (1991), with a number of publications produced until this date. An example of results can be seen in Fig. 43, compared with two other methods discussed elsewhere in this review. Some important numerical developments include the following. Gieseler et al. (2000) introduced an injection scheme based on the thermal leakage model. Kang et al. (2001) implemented a grid-based AMR scheme, as well as sub-zone shock tracking, in order to address realistic diffusion coefficients at a manageable computational cost. Jones and Kang (2005) applied a Coarse-Grained Finite Momentum-Volume Scheme (CGFMV), an extension of the piecewise power-law method introduced by Jones et al. (1999) already presented above. The basic idea is to lower the numerical resolution in momentum (down to as few as two to three bins per decade), but prescribe the shape of the spectrum in each bin so as to maintain reasonable accuracy; the numerical spectrum is then no longer a piecewise constant function but a piecewise linear function. All the works cited so far were restricted to slab geometry, Kang and Jones (2006) simulated shocks in spherical geometry—still effectively one-dimensional, under the assumption of spherical symmetry—in order to study CR feedback at SNR shocks. For this the authors employed a frame comoving with the outer shock, which was found to lower the convergence requirements. Following the formalism of cosmology, they use coordinate $\tilde{r} = r/a$ where $a(t)$ is the expansion factor, with expansion rate obtained from the
measured shock speed in the Cartesian grid. Then density and pressure (or energy density) are re-scaled as $\tilde{\rho} = \rho a^3$ and $\tilde{P} = Pa^3$ (the time variable is left unchanged). For quasi-parallel plane shocks, Kang and Jones (2007); Kang et al. (2009) found a self-similar evolution and proposed analytic forms for the solutions during that stage. All these improvements allowed the team to investigate the time evolution of CR-modified shocks and particle spectra, and eventually compute the non-thermal radiation from the accelerated particles in SNRs (Edmon et al. 2011; Kang et al. 2012). The latest developments are the inclusion of prescriptions for magnetic field amplification (MFA, see Sect. 2.2.4) and of the effects of the Alfvénic drift (see Eq. (88)) by Kang (2012), Kang (2013), Kang et al. (2013), with a formalism similar to the approaches discussed later in Sect. 5.7.2.

**Other approaches in spherical geometry** To account for the spatially-dependent diffusion of energetic particles, an alternative approach to AMR is to perform a change of variables, as done by Ptuskin et al. (2010); Zirakashvili and Ptuskin (2012) and by Telezhinsky et al. (2012b, a). These simulations are performed in spherically symmetric geometry, and notably include the reverse shock and the contact discontinuity as well as the forward shock (located at radius $R_b$, $R_c$, and $R_f$ respectively).

In Ptuskin et al. (2010), Zirakashvili and Ptuskin (2012) the radius $r$ is replaced by $r/R_b$ for $r < R_b$ and $r/R_f$ for $r > R_f$ (unshocked medium), and by $(r - R_c)/(R_c - R_b)$ for $R_b < r < R_c$ and $(r - R_c)/(R_f - R_c)$ for $R_c < r < R_f$ (shocked region). Also time is replaced by a dimensionless parameter that scales as $\ln(R_f)$. The hydrodynamic equations are solved using FDM, separately in the upstream regions (using an implicit scheme since the flow is supersonic) and in the downstream regions (using an explicit scheme), and the three discontinuities between the different regions are moved manually in time. As in other methods, the diffusion-convection equation for particles is recast as a set of tri-diagonal equations. In these works the CR pressure is included in the HD equations, and the velocity used for particles includes the Alfvénic drift. An example of the results is shown in Fig. 37, for a fiducial SNR at age 1000 year.

In Telezhinsky et al. (2012b, a) the radius is first normalized to the outer shock radius: $x = r/R_f$ (a comoving coordinate, as in Kang and Jones 2006) and then transformed according to $(x-1) = (x_* - 1)^3$, so that $dx/dx_* = 3(x_* - 1)^2$, with a uniform binning in new coordinate $x_*$. The transport equation for the particles is solved using an implicit FDM, while the HD equations are solved with the VH-1 code. These simulations were made in the test-particle regime, they do not include the back-reaction of energetic particles on the hydrodynamics. Telezhinsky et al. (2012a) discuss the possible contribution to the SNR emission of CRs accelerated at the reverse shock w.r.t. those accelerated at the forward shock.

### 5.4.2 Coupling MHD and Vlasov-Fokker-Planck equations

Bell et al. (2011), Reville and Bell (2013) develop a method based on performing a spherical harmonics expansion of the distribution $f$ in momentum $p$, using the KALOS code previously developed for simulations of laser-plasma interactions (Bell et al. 2006). The authors consider a 1D problem with a planar shock (with the normal to the shock front oriented in the $x$ direction with $x < 0$ the upstream medium) and keep the
<figure><img src="image_63.png" /><figcaption>Fig. 37 Representative results from a hydro-kinetic simulation of a SNR in spherical geometry, at t = 1000 year. The supernova has typical mass 1.4 M$_{\odot}$ and energy 10$^{51}$ erg, released in a uniform medium of density 0.1 cm$^{-3}$, temperature 10$^{4}$ K, and magnetic field 5 µG. Left: radial dependencies of gas density (thick solid line), gas velocity (dotted line), CR pressure (thick dashed line), and gas pressure (dashed line). The reverse and forward shocks are visible as discontinuities in the hydro profiles. An extended CR precursor is visible ahead of each shock. Right: spectra of accelerated particles: protons at the forward shock (thick solid line), ions at the reverse shock (thick dashed line), electrons at the forward shock (×100, thin solid line), and positrons at the reverse shock (×100, thin dashed line) present from radioactive decay in the ejecta. Note the concavity of the spectra, plotted as $p^4 f(p)$. Image reproduced with permission from Zirakashvili and Ptuskin (2012), copyright by Elsevier</figcaption></figure>

information about the momentum vector **p**. They solve the following Vlasov–Fokker–Planck (VFP) equation (Bell et al. 2011)

$$ \partial_t f(x, \mathbf{p}, t) + (v_x + u)\partial_x f - \partial_x u \, p_x \, \partial_{p_x} f - \frac{u}{c} \partial_x u \, p_x \, \partial_{p_x} f + q\mathbf{v} \times \mathbf{B} \cdot \partial_{\mathbf{p}} f = \mathcal{C}(f), \quad (96) $$

where $v$ and $u$ stand for the particle and the fluid speeds, and $q$ is the particle charge. The particle distribution function is measured in the fluid rest-frame. **B** is the local magnetic field. The term $\mathcal{C}(f)$ is a collision term; in the collisionless plasmas investigated here this term is due to wave-particle interactions. It is written

$$ \mathcal{C}(f) = \partial_{\mu} v_s (1 - \mu^2) \partial_{\mu} f(x, \mathbf{p}, t), \quad (97) $$

where $\mu = \cos(\theta)$ and $v_s$ is a parameter to be scaled with respect to the particle gyro-frequency $\Omega = qB/\gamma mc$. In this model particle angular deflection is assumed to be small. This justifies the use of the term Fokker–Planck.

The particle distribution is then expanded into spherical harmonics as

$$ f(x, \mathbf{p}, t) = \sum_{l,m} f_l^m(x, p, t) P_l^{|m|}(\cos \theta) \exp(-im\phi), \quad (98) $$

where $l$ is a positive integer, $m$ is an integer between $-l$ and $+l$, $\theta$ is the particle momentum pitch-angle with respect to the shock normal, $\phi$ is the particle momentum azimuthal angle, and $P_l^m(\cos \theta)$ are the associated Legendre polynomials of order $l$. We have $f_l^{-m} = (f_l^m)^*$, the complex conjugate of $f_l^m$. By applying this expansion
to the VFP equation, one obtains a hierarchy of equations for each component $f_l^m$, not reproduced here [see Eq. (2) of Bell et al. 2011]. Usually, as in shocks the particle distribution is close to isotropy due to efficient wave-particle scattering it is only necessary to retain the first few terms $f_0^0$, $f_1^0$, and $f_1^1$. $f_0^0$ represents the isotropic part of the distribution function while $f_1^0$ represents the particle flux along the shock normal. $f_1^1$ complements the construction of the CR current. In the limit of a small ratio $u/c$ the system reduces to a simple system using the Chapman–Enskog expansion $f_l^m \sim (u/c)^l f_0^0$ (with $v \simeq c$). In the upstream medium

$$ \begin{aligned} & \frac{c}{3} \partial_x f_0^1 + u \partial_x f_0^0 = 0, \\ & c \partial_x f_0^0 + 2 \Omega_z Re(f_1^1) = -v_s f_1^0, \\ & -i \Omega_x f_1^1 - \frac{\Omega_z}{2} f_1^0 = -v_s f_1^1, \end{aligned} \quad (99) $$

while in the downstream medium the particle distribution is isotropic, so limited to the term $f_0^0$. The solution at the shock front can be obtained using the continuity of the particle distribution and its flux at the shock front, accounting for the frame transformation across the shock. Looking for a power-law solution $f \propto p^{-\gamma}$, at the leading order in $u/c$ one finds

$$ \gamma = 3 + 3 \frac{u_2}{u_{\text{sh}} - u_2} \frac{f_0^0(+\infty)}{f_0^0(x=0)}, \quad (100) $$

where $u_2$ is the fluid speed and $f_0^0(+\infty)$ the particle isotropic component far downstream. One recovers $\gamma = 4$ for a strong unmodified shock.

**Oblique shock solutions** Bell et al. (2011) apply the above formalism to an oblique flow, where the magnetic field direction is given by the angle $\theta$ with respect to the shock normal. The procedure requires to start with a guess of the particle distribution at the shock front as a power-law. The index $\gamma$ is calculated iteratively with the full set of equations for $f_l^m$. These are solved using FDM techniques (see Sect. 3.6.1) with $f_l^m$ calculated at the cell centers for even $l$ and at the cell boundaries for odd $l$. The spatial resolution is finer close to the shock front and coarser at the edge of the simulation box. Non-linear CR backreaction is not considered in this work: the background magnetic field and fluid velocity are not modified by the CR pressure. The shock velocity profile is smoothed using a hyperbolic tangent profile of width $x_s$ of 1% of the particle Larmor radius upstream. The shock index solutions show soft distributions at quasi-perpendicular fast shocks. Such configurations can occur if strong magnetic field amplification occurs after the supernova shock breakout (Marcowith et al. 2018) or/and if the shock propagates in circumstellar medium with a spiral magnetic field (Bell 2008).

**CR driven instabilities** Reville and Bell (2013) and Bell et al. (2013) couple the VFP method described above to MHD solutions in order to calculate the CR pressure backreaction over the shock solutions. The authors adopt mixed coordinate frames (Skilling

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