# Multi-scale simulations of particle acceleration in astrophysical systems

Alexandre Marcowith, Gilles Ferrand, Mickael Grech, Zakaria Meliani, Illya

Plotnikov, Rolf Walder

► To cite this version:

Alexandre Marcowith, Gilles Ferrand, Mickael Grech, Zakaria Meliani, Illya Plotnikov, et al. Multi-scale simulations of particle acceleration in astrophysical systems. Living Reviews in Computational Astrophysics, 2020, 6, pp.1. 10.1007/s41115-020-0007-6. hal-03004333

HAL Id: hal-03004333

https://hal.science/hal-03004333

Submitted on 13 Nov 2020

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Multi-scale simulations of particle acceleration in astrophysical systems

Alexandre Marcowith¹ · Gilles Ferrand² · Mickael Grech³ · Zakaria Meliani⁴ · Illya Plotnikov⁵ · Rolf Walder⁶

Received: 30 July 2019 / Accepted: 25 January 2020

© The Author(s) 2020

## Abstract

This review aims at providing an up-to-date status and a general introduction to the subject of the numerical study of energetic particle acceleration and transport in turbulent astrophysical flows. The subject is also complemented by a short overview of recent progresses obtained in the domain of laser plasma experiments. We review the main physical processes at the heart of the production of a non-thermal distribution in both Newtonian and relativistic astrophysical flows, namely the first and second order Fermi acceleration processes. We also discuss shock drift and surfing acceleration, two processes important in the context of particle injection in shock acceleration. We analyze with some details the particle-in-cell (PIC) approach used to describe particle kinetics. We review the main results obtained with PIC simulations in the recent years concerning particle acceleration at shocks and in reconnection events. The review discusses the solution of Fokker-Planck problems with application to the study of particle acceleration at shocks but also in hot coronal plasmas surrounding compact objects. We continue by considering large scale physics. We describe recent developments in magnetohydrodynamic (MHD) simulations. We give a special emphasis on the way energetic particle dynamics can be coupled to MHD solutions either using a multi-fluid calculation or directly coupling kinetic and fluid calculations. This aspect is mandatory to investigate the acceleration of particles in the deep relativistic regimes to explain the highest cosmic ray energies.

**Keywords** Kinetic physics · Magnetohydrodynamics · Turbulence · Particle acceleration

## Contents

<table><tr><td>1 Introduction</td><td></td></tr><tr><td>1.1 Layout</td><td></td></tr><tr><td>1.2 List of acronyms and notations</td><td></td></tr><tr><td>2 Astrophysical and physical contexts</td><td></td></tr></table>

<footnote>Extended author information available on the last page of the article</footnote>
2.1 Stochastic acceleration

2.2 Diffusive shock acceleration

2.2.1 Fermi processes and building power-laws

2.2.2 The transport equation and the diffusion coefficient

2.2.3 The injection problem

2.2.4 Back-reaction and non-linear effects

2.3 Shock drift and shock surfing acceleration

2.3.1 Shock drift acceleration: SDA

2.3.2 Shock surfing acceleration: SSA

2.4 Fermi acceleration process at relativistic shocks

2.4.1 General statements

2.4.2 Progress with fully kinetic simulations

2.4.3 Long term evolution

2.5 Reconnection in astrophysical flows

2.5.1 Collisional reconnection models

2.5.2 Collisionless reconnection

2.5.3 Other effects

2.6 Laser plasma experiments

2.6.1 Overview of laser facilities and characteristic parameters

2.6.2 The collisionless regime

2.6.3 Collisionless shock experiments

2.6.4 Magnetic reconnection experiments

3 Solving kinetic problems

3.1 The Vlasov–Maxwell description of a collisionless plasma

3.1.1 Governing equations

3.1.2 Initial and boundary conditions

3.2 Solving the Vlasov–Maxwell system numerically: general considerations

3.2.1 Initialization and time-loop

3.2.2 Brief discussion of Maxwell solvers

3.3 Particle-In-Cell codes

3.3.1 Method

3.3.2 Additional physics modules

3.3.3 Stability issues: relativistic flows

3.3.4 Examples

3.4 Vlasov–Maxwell codes

3.4.1 PIC codes versus Vlasov codes

3.4.2 The problem of filamentation in phase-space

3.4.3 Example of different methods for electromagnetic Vlasov codes

3.5 Hybrid methods

3.6 Solving Fokker–Planck problems

3.6.1 The Fokker–Planck equation

3.6.2 Stochastic differential equations

3.6.3 Simulating shock acceleration using a Monte-Carlo method

4 Small and meso-scale numerical particle acceleration studies

4.1 Shock acceleration numerical experiments: PIC simulations

4.1.1 Ultra-relativistic shocks

4.1.2 Mildly relativistic shocks

4.1.3 Non-relativistic shocks

4.1.4 Discussion

4.2 Kinetics of magnetic reconnection

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

4.2.2 Some key results based on kinetic simulations

4.2.3 Selected physical processes which accelerate particles

4.2.4 A critical discussion and outlook

5 Macro scale numerical particle acceleration studies

5.1 The equations of magnetohydrodynamics

5.1.1 Classical magnetohydrodynamics
<table><tr><td>5.1.2</td><td>Equation of state (EOS)</td></tr><tr><td>5.2</td><td>Numerical solutions</td></tr><tr><td>5.2.1</td><td>The Courant-Friedrichs-Lewy condition</td></tr><tr><td>5.2.2</td><td>Riemann solvers</td></tr><tr><td>5.2.3</td><td>Semi-implicit and implicit schemes</td></tr><tr><td>5.2.4</td><td>Magnetic divergence-free algorithms</td></tr><tr><td>5.2.5</td><td>Adaptive mesh refinement techniques</td></tr><tr><td>5.2.6</td><td>Errors estimator</td></tr><tr><td>5.2.7</td><td>Load balance</td></tr><tr><td>5.2.8</td><td>Standard numerical tests</td></tr><tr><td>5.2.9</td><td>Relativistic magnetohydrodynamics</td></tr><tr><td>5.3</td><td>Multi-fluid methods</td></tr><tr><td>5.3.1</td><td>Model equations</td></tr><tr><td>5.3.2</td><td>Specific numerical schemes of CR-fluid systems</td></tr><tr><td>5.3.3</td><td>Numerical tests</td></tr><tr><td>5.3.4</td><td>Cosmic ray physics using a bi-fluid approach</td></tr><tr><td>5.4</td><td>The fluid-kinetic framework</td></tr><tr><td>5.4.1</td><td>Coupling hydro and diffusion-convection equations</td></tr><tr><td>5.4.2</td><td>Coupling MHD and Vlasov-Fokker-Planck equations</td></tr><tr><td>5.5</td><td>Particle-in-cell-magnetohydrodynamics</td></tr><tr><td>5.5.1</td><td>Numerical schemes</td></tr><tr><td>5.5.2</td><td>Coupling MHD and PIC time steps</td></tr><tr><td>5.5.3</td><td>Adaptive mesh refinement for the PIC module</td></tr><tr><td>5.5.4</td><td>Numerical tests</td></tr><tr><td>5.5.5</td><td>Results on cosmic ray acceleration at non-relativistic shocks</td></tr><tr><td>5.5.6</td><td>Other astrophysical applications</td></tr><tr><td>5.6</td><td>A list of HD and MHD codes with CR physics</td></tr><tr><td>5.7</td><td>Semi-analytical approaches for cosmic ray acceleration</td></tr><tr><td>5.7.1</td><td>General considerations</td></tr><tr><td>5.7.2</td><td>Construction of a NLDSA model</td></tr><tr><td>6</td><td>Summary and conclusions</td></tr><tr><td></td><td>References</td></tr></table>

# 1 Introduction

Particle acceleration is a widespread process in astrophysical, space and laser plasmas. Acceleration results from the effect of electric fields, but supra-thermal particles gain energy because their residence time in the acceleration zone increases due to magnetic confinement. Hence, particle acceleration is an electromagnetic process. Particle acceleration can be classified into three main sub-types (Blandford 1994; Kirk 1994; Melrose 1996): acceleration at flow discontinuities among which shock waves, stochastic acceleration, acceleration by direct electric fields. This last mechanism occurs in the environment of fast rotating magnetized objects like pulsars or planetary magnetosphere. It will not be discussed in this review, interested readers can refer to Cerutti and Beloborodov (2017) and references therein for what concerns pulsar magnetospheric physics. Some recent discussions concerning particle acceleration in Jupiter, the fastest rotator among solar system planets and other giant planet magnetospheres can be found in Mauk et al. (2017), Delamere et al. (2015).

As stated above acceleration processes can be classified into different categories. Let us give here a short overview of the main mechanisms.
1. **Stochastic Fermi acceleration (SFA).** This is historically the first discovered acceleration process. SFA is at the heart of Fermi's work on the origin of cosmic rays (CRs) (Fermi 1949, 1954). In its original description particles with speed $v$ gain energy through a stochastic interaction with the convective electric field carried by magnetic clouds moving randomly at a speed $U \ll v$ in the interstellar space (see Sect. 2.2.1). This process has some well-known caveats: it is inefficient as the mean relative energy gain $\langle \Delta E/E \rangle \propto (U/v)^2$ (Parker 1958a), it produces non-universal power-law solutions (see Eq. 3), which can not explain the power-law distribution of CRs observed at Earth. A modern description of SFA includes randomly moving electromagnetic waves (Hall and Sturrock 1967), usually in the magnetohydrodynamic (MHD) limit if we want to consider the issue of CR acceleration (Parker 1955; Kulsrud and Pearce 1969). For SFA by MHD waves, the relevant speed for the scattering centers usually are proportional to the local Alfvén speed $U_A$.$^1$ SFA is more efficient if the speed of the scattering center is close to the speed of light (Marcowith et al. 1997; Gialis and Pelletier 2004) or for particles with speeds close to $U_A$, or other plasma characteristic wave phase speeds as it is likely the case for low-energy CRs which propagate in the interstellar medium (ISM) (e.g., Ptuskin et al. 2006), in the solar corona (e.g., Miller et al. 1996; Pryadko and Petrosian 1998) or in the heliosphere (Zhang and Lee 2013). SFA is a multi-scale process if the scattering waves are distributed over large wave-number bands as it is the case in turbulent flows. SFA is an important process at the origin of particle acceleration and gas heating in hot corona which develop around compact objects (see Sect. 3.6).

To finish, we also mention the betatron-magnetic pumping process (Parker 1958a). In this process, particles propagate along a magnetic field slowly varying with time. As the field strength increases, the particle momentum increases due to the conservation of the adiabatic invariant $p_⊥^2/B$.² If the particle suffers some scattering of its pitch-angle or some elastic collisions, a decrease of the magnetic field strength while keeping particle isotropization by pitch-angle scattering produces a net gain in energy during a magnetic field variation cycle.

2. **Shock or shear-flow acceleration.** One way to cure the inefficiency of the SFA is to allow scattering centers to have a mean direction of motion (Parker 1958a; Wentzel 1963). Then the mean relative energy gain scales as $\langle \Delta E/E \rangle \propto (U/v)$. This case occurs for a shock because of the advection of the scattering centers towards the shock front but also in the configuration of a shearing, for instance in jets. We do not explicitly for now discuss the case of shear-flow acceleration, the interested reader can refer to Rieger and Duffy (2006), and we consider hereafter the case of shock acceleration. Shear-flow acceleration will be reviewed in a forthcoming version of the text. The acceleration process is more efficient if particles can reside for a sufficient amount of time around the shock front (Kirk 1994). In fact, particle acceleration at shock waves covers three basic different processes (see Kirk 1994; Treumann and Jaroschek 2008b; Marcowith et al. 2016): diffusive shock

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

<footnote>$^1$ $U_A = B/\sqrt{4\pi\rho_i}$, where $B$, $\rho_i$ are respectively the local magnetic field strength, the ion mass density.</footnote>

<footnote>$^2$ $p_⊥ = p \sin \alpha$, where $p$ and $\alpha$ are the particle momentum and pitch-angle (the angle between the particle velocity and magnetic field) respectively.</footnote>
acceleration (DSA), shock drift acceleration (SDA), shock surfing acceleration (SSA), all described below.

Before, let us introduce some elements of vocabulary associated with different descriptions of the shock front. First, particle acceleration requires the shock to be collisionless, i.e., to be mediated by electromagnetic processes rather than collisions otherwise collisions being the fastest process force the shocked particle distribution to be Maxwellian. This condition is usually fulfilled in astrophysical and space plasmas. At a macroscopic level a shock wave is characterized by a discontinuity in the thermodynamical variables of the flow. A shock occurs when the flow is supersonic, with a sonic Mach number $M_s = U_{sh}/c_s > 1$, where $U_{sh}$ is the shock speed in the upstream medium restframe. Rankine–Hugoniot conditions give the jump in the flow density, velocity, pressure, temperature and entropy at the shock front (Treumann and Jaroschek 2008a). At a microscopic level a shock front is a complex, dynamical structure where multi-scale instabilities can develop (Marcowith et al. 2016). The so-called supercritical magnetized shock front is composed of three sub-structures: the foot, a bump in gas and magnetic field pressures due to the accumulation of ions reflected at the ramp, the ramp which marks the fast rise of the electromagnetic field potential and gas density and finally the overshoot-undershoot produced by the gyromotion of reflected ions moving in the post-shock gas. Shocks in astrophysics also include precursors in the upstream gas which can have very different origins: radiation, mixture of ions and neutrals or CRs. The size of these precursors make them usually impossible to explore numerically with the shock front structure as a single complex dynamical system.

As stated above particle acceleration in shocks proceeds through three different mechanisms. In DSA particles repeatedly gain energy by crossing the shock front back and forth (Drury 1983). SDA results from the effect of the convective electric field $\mathbf{E} = -\mathbf{U}/c \times \mathbf{B}$ upstream the shock front due to the motion of the flow at a speed $\mathbf{U}$ (Kirk 1994; Decker and Vlahos 1985). The particle guiding-center drifts due to the effect of the electric field and to the gradient of the magnetic field in the ramp. SSA results from the trapping of the particle at the shock front because of the combined effect of shock potential raise at the ramp and the convective upstream electric field (Sagdeev 1966). We will come back with more details on these mechanisms in Sects. 2.2 and 2.3.

3. Magnetic reconnection (REC). Magnetic reconnection is the process which transfers magnetic field energy into kinetic energy in an explosive event by re-arranging the magnetic field topology. The most simple 2D picture is sketched in the left panel of Fig. 1. Separatrices (green, dashed lines) divide the 2D plane into 4 different regions: in the left region, the magnetic field connects points A and B. In the right region, the field connects points A' and B' and is oriented in opposite direction to the field in the left region. No field is present in the upper and lower region between the separatrices. In 2D, due to Ampère's law, a current pointing normal to the plane is necessarily present between the oppositely oriented fields. The REC process then leads to a re-arrangement of the field lines, lowering the magnetic energy. The field now connects the points A and A' (B and B' respectively). These field lines are highly bent and will relax, accelerating the plasma upwards and
<figure><img src="image_4.png" /><figcaption>**Fig. 1** The re-arrangement of the field topology in magnetic reconnection in a 2D model. Before the event (Image 1), points A and B (A' and B' respectively) are located on the same field line. After the event (3), field lines now are connecting points A and A' (B and B' respectively). Strongly accelerated outflow is driven in the directions where the highly bent magnetic field lines are relaxing. The dashed green lines are called separatrices, lines which separate the regions of field which are topologically not connected. Image adapted from Melzani (2014)</figcaption></figure>

downwards. The term magnetic reconnection was first coined in Dungey (1958) and was later adopted by the community. More details about the REC process are provided in Sect. 2.5.

REC induces a transfer of magnetic energy into: heat, plasma and particle acceleration and hence radiation (Gonzalez and Parker 2016; Priest 1994). Particle acceleration in reconnection sites can either occur by a direct acceleration in electric fields in the current sheet, or because of Fermi first order acceleration in the plasma converging towards the reconnection zone or if particles are trapped in a contracting magnetic islands (de Gouveia Dal Pino and Kowal 2015). The physics of particle acceleration in kinetic reconnection is discussed in Sect. 4.2.

Acceleration mechanisms, as we have seen from the above rapid descriptions, intrinsically involve multi-scale processes which bring particles from the thermal to supra-thermal speeds. In astrophysics these processes have to explain the CR spectrum observed at Earth which extends at least over 15 orders of magnitude in energy (from MeV to ZeV) and more than 30 orders of magnitude in flux (see Fig. 2). Note that in space plasmas, maximum energies reached by the energetic particles are more modest but still supra-thermal, and the particle distributions cover about 5 orders of magnitude (from keV to GeV) (see, e.g., Zharkova et al. 2011 in the context of solar flares). The investigation of particle acceleration then requires different numerical approaches to probe the different inter-connected scales involved in the process of acceleration. Multiple techniques are also required as actually it is not possible to account for such large dynamical spatial, time and energy scales even with modern computers. It is the main object of this review to address these different techniques.

## 1.1 Layout

This review is organized as follows. Section 2 addresses the scientific context. It describes the main acceleration processes at work in astrophysical plasma systems.

<figure><img src="image_5.png" /></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 2** The cosmic ray spectrum observed at the Earth multiplied by $E^{2.6}$. Image reproduced with permission from Patrignani and Particle Data Group (2016), copyright by Regents of the University of California</figcaption></figure>

It also a short review on the on-going experimental efforts to reproduce collisionless shocks, magnetic reconnection and particle acceleration in laser-plasma-based experiences. The next sections treat the different numerical approaches to investigate particle acceleration from microscopic scales to macroscopic scales. Section 3 describes the different numerical methods adapted to the description of plasma kinetics. Section 4 discusses particle acceleration and transport at micro- and meso-plasma scales. Section 5 describe numerical techniques developed to follow macroscale dynamics and detail recent results on particle acceleration and transport in astrophysical flows. We conclude in Sect. 6.

**How to read this review** The scientific questions and the numerical experiments developed to investigate them are entangled. We have decided to describe this complex modelling in two steps. The first step presents a general description of the main acceleration processes in astrophysical plasmas. This presentation is the main purpose of Sect. 2. Notice that we complement it by a dedicated section addressing some recent studies of particle acceleration at collisionless shocks and magnetic reconnection in the context of laser plasmas in Sect. 2.6. The second step describes technical numerical aspects. They are presented in Sects. 3.4 to 3.5 and in all sub-sections of Sect. 5. Sections 3, 4, 5 then include discussions which connect the numerical work and scientific questions exposed in Sect. 2.

## 1.2 List of acronyms and notations

All quantities are in cgs Gaussian units.
<table><thead><tr><td>Acronym name</td><td>Definition</td></tr></thead><tbody><tr><td>AGN</td><td>Active galactic nucleus</td></tr><tr><td>AMR</td><td>Adaptive mesh refinement</td></tr><tr><td>CFL</td><td>Courant–Friedrichs–Lewy</td></tr><tr><td>CR</td><td>Cosmic Ray</td></tr><tr><td>DCE</td><td>Diffusion-convection equation</td></tr><tr><td>DSA</td><td>Diffusive shock acceleration</td></tr><tr><td>EP</td><td>Energetic particle</td></tr><tr><td>FDM</td><td>Finite difference method</td></tr><tr><td>FVM</td><td>Finite volume method</td></tr><tr><td>FP</td><td>Fokker–Planck</td></tr><tr><td>GRB</td><td>Gamma-ray burst</td></tr><tr><td>HD</td><td>Hydrodynamics</td></tr><tr><td>ISM</td><td>Interstellar medium</td></tr><tr><td>MFA</td><td>Magnetic field amplification</td></tr><tr><td>MHD</td><td>Magneto-hydrodynamics</td></tr><tr><td>NLDSA</td><td>Non-linear diffusive shock acceleration</td></tr><tr><td>PDE</td><td>Partial differential equation</td></tr><tr><td>PIC</td><td>Particle-in-cell</td></tr><tr><td>PWN</td><td>Pulsar wind nebula</td></tr><tr><td>REC</td><td>Magnetic reconnection</td></tr><tr><td>SDA</td><td>Shock drift acceleration</td></tr><tr><td>SFA</td><td>Stochastic Fermi acceleration</td></tr><tr><td>SNR</td><td>Supernova remnant</td></tr><tr><td>SSA</td><td>Shock surfing acceleration</td></tr></tbody></table>

We recall here the definition of the different quantities used to construct the main parameters involved in shock acceleration processes: *B* is the magnetic field strength, *ρ* is the gas mass density, *ρ*<sub>i</sub> is the ion mass density, *γ*<sub>ad</sub> is the gas adiabatic index, *v* and *p* are the charged particle speed and momentum and *Ze* is the charge.

<table><thead><tr><td>Notation</td><td>Definition</td></tr></thead><tbody><tr><td>θ<sub>B</sub></td><td>Shock magnetic field obliquity<br>(Angle between field lines and shock normal)</td></tr><tr><td>U<sub>sh</sub></td><td>Shock velocity in the upstream (observer) frame</td></tr><tr><td>U<sub>A</sub> = B/√4πρ<sub>i</sub></td><td>Local Alfvén speed</td></tr><tr><td>M<sub>A</sub> = U<sub>sh</sub>/U<sub>A</sub></td><td>Alfvénic Mach number</td></tr><tr><td>c<sub>s</sub> = √γ<sub>ad</sub>P/ρ</td><td>Local sound speed</td></tr><tr><td>M<sub>s</sub> = U<sub>sh</sub>/c<sub>s</sub></td><td>Sonic Mach number</td></tr><tr><td>σ = B<sup>2</sup>/4πρc<sup>2</sup></td><td>Local magnetization<br>(Ratio of the upstream magnetic pressure to the upstream gas kinetic energy) (shocks)</td></tr><tr><td>σ = ω<sub>ci</sub><sup>2</sup>/ω<sub>pi</sub><sup>2</sup></td><td>Local magnetization<br>(Ratio of the square of cyclotron to plasma frequencies) (reconnection)</td></tr><tr><td>β<sub>p</sub></td><td>Plasma parameter</td></tr><tr><td>r<sub>L</sub> = pv/ZeB</td><td>Gyro-radius (Larmor radius)</td></tr></tbody></table>

Notice of caution: There is not a unique way to define the magnetization parameter *σ*. In shock studies, *σ* is the ratio of the magnetic pressure to the kinetic energy of the ambient gas, all quantities being measured in the upstream rest-frame. In relativistic

<figure><img src="image_7.png" /></figure>
shock studies, the magnetization parameter is sometimes defined as $\sigma = B^2/4\pi\Gamma U c^2$, where $\Gamma = \sqrt{1+U^2}$ is the Lorentz factor of the shock and $U$ is the four-velocity of the flow (Marcowith et al. 2016). In reconnection studies, $\sigma$ is the ratio of the square of cyclotron to plasma frequencies. It is related to the ratio $r_A = U_A/c$ by $\sigma = r_a^2/(1-r_A^2)$ (e.g., Sironi and Beloborodov 2019).

## 2 Astrophysical and physical contexts

Following the introductory remarks in Sect. 1, we present below an overview of the astrophysical and physical contexts where the numerical tools discussed in this review are actively developed. We aim here at a short description of the basic concepts necessary to describe particle acceleration. In particular, we show that particle acceleration involves a large range in scale/time/energy which justifies the use of very different numerical techniques detailed in the next sections.

First, in Sect. 2.1 we briefly overview the mechanism of stochastic acceleration (SA). Then the three next sections cover different aspects of the physics of particle acceleration at collisionless shocks. In Sect. 2.2 we provide a general and rather detailed presentation of the physics of diffusive shock acceleration (DSA) which is one of the main frameworks to study particle acceleration in astrophysical systems. Beyond a standard description of the process itself we discuss specific issues connected with the acceleration of cosmic rays (CRs) at fast astrophysical shock waves: the injection problem and non-linear back-reaction of CRs over the flow solution. These two difficulties require the development of specific scale-dependent numerical techniques described in the next sections. Section 2.3 is a short presentation of the other two shock acceleration processes, namely the shock drift acceleration (SDA) and the shock surfing acceleration (SSA) which are especially relevant for particle injection in the DSA process. Section 2.4 discusses the specific case of Fermi acceleration at relativistic shocks and the development of micro turbulence at these shock fronts. Magnetic reconnection (REC) is discussed in some detail in Sect. 2.5, where we present the most relevant vocabulary necessary to understand particle acceleration in reconnection structures. Section 2.6 reviews the most important undergoing or planned laser experiments to study particle acceleration. This rapidly growing field of research starts now to investigate astrophysically relevant conditions for particle acceleration at collisionless shocks and magnetic reconnection. Notice that we decided to not include any review of acceleration processes in space plasmas, this will deserve a special section in a forthcoming version.

It should be stressed that *by no means this section is intended to be exhaustive.* It has to be understood as a short introduction to the scientific cases where the different simulation techniques discussed hereafter are developed. For each type of acceleration/transport mechanism we refer the interested reader to more complete dedicated reviews.
## 2.1 Stochastic acceleration

As discussed in Sect. 2.2.1 stochastic acceleration occurs because, on average, energetic particles at a speed $v$ interact with scattering centers moving at a speed $U$ more often through head-on collisions than through rear-on collisions if $v \gg U$. This results in a broadening of the particle distribution and an increase of the mean particle energy (Melrose 1980). In astrophysical plasmas the scattering centers often<fnref n="3" /> can be described as plasma waves, and when we deal with high-energy CRs these waves can be described using the MHD approximation (Parker 1955; Sturrock 1966; Kulsrud and Ferrari 1971). But it is necessary to go beyond MHD if we want to consider the acceleration of non-relativistic or mildly relativistic particles (Marcowith et al. 1997; Pryadko and Petrosian 1998).

As explained in Sect. 1, well-known caveats prevent the interpretation of the CR spectrum observed at the Earth as resulting from stochastic acceleration by MHD waves: (1) the non-universality of the distribution of accelerated particles, (2) a weak relative energy gain at each wave-particle interaction scaling as $(U/v)^2$. The second issue can be partly overcome if we consider the case of low energy (sub-GeV) CR propagation in the ISM, as in that case the ratio $v/U$ drops. Still, an important problem results in the prohibitive amount of ISM turbulence necessary to re-accelerate the low energy end of CR spectrum (Ptuskin et al. 2006; Thornbury and Drury 2014; Drury and Strong 2017). Nevertheless, SA has been invoked to be an important source of turbulence damping and particle acceleration in solar flares (Petrosian 2012), in active corona above the accretion discs around compact objects (Dermer et al. 1996; Liu et al. 2004; Belmont et al. 2008; Vurm and Poutanen 2009), in SNRs or their associated superbubbles (Bykov and Fleishman 1992; Kirk et al. 1996; Marcowith and Casse 2010; Ferrand and Marcowith 2010), in galaxy clusters (Brunetti and Lazarian 2007), or in the case the wave phase (Alfvén) speed gets close to the speed of light as can be the case in AGNs (Henri et al. 1999), in GRBs (Schlickeiser and Dermer 2000), or in pulsar winds (Bykov et al. 2012).

## 2.2 Diffusive shock acceleration

DSA is probably the favored production mechanism of CRs. It is thought to be a natural outcome of collisionless shocks, and so is believed to be at work in astrophysical shocks at all scales, active in the bow shocks in the solar system, in the blast waves of supernova remnants (SNRs) or in the internal shocks in the jets of gamma-ray bursts (GRBs) or active galactic nuclei (AGNs). DSA is rooted in the early ideas of Fermi (1949, 1954); it was developed independently in the late 1970s by Krymskii (1977), Axford et al. (1977), Bell (1978a, b), Blandford and Ostriker (1978), see Drury (1983), Jones and Ellison (1991), Malkov and Drury (2001) for comprehensive reviews. A key feature of DSA is that it produces power-law distributions as a function of energy (although this

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

<footnote n="3">³ Note however that, similarly to Fermi’s original ideas, some models invoke finite amplitude waves or shocks as the origin of stochastic energy exchanges with energetic particles, see in different astrophysical contexts e.g. Bykov and Toptygin (1987), Achterberg (1990), Pelletier and Marcowith (1998), Gialis and Pelletier (2004).</footnote>
spectrum can be altered by non-linear effects), which is similar to the CR spectrum as observed from the Earth, modulated by CR transport and escape from the Milky Way. DSA requires two ingredients to accelerate particles: a converging flow (the shock wave), and scattering centers (perturbations of the magnetic field). In this mechanism individual microscopic particles can be accelerated up to very high energies because they are interacting a large number of times with the macroscopic shock discontinuity before escaping the system. Again, one major difficulty in simulating this process becomes apparent: DSA is intrinsically a multi-scale problem.

### 2.2.1 Fermi processes and building power-laws

In his original model, Fermi considered the interaction of charged particles with moving magnetized clouds. In the cloud frame, the particle (of velocity $v$) is elastically deflected around the $B$ field. In the Galactic frame (with respect to which the cloud is moving at velocity $U$), the energy $E$ of the particle changes according to

$$ \frac{\Delta E}{E} = -2 \frac{v \cdot U}{c^2}. \qquad (1) $$

The effect depends on the geometry of the encounter: for head-on collisions ($v \cdot U < 0$) the particle gains energy ($\Delta E > 0$), whereas for overtaking collisions ($v \cdot U > 0$) it loses energy ($\Delta E < 0$). The exchange is mediated by the magnetic field, even though $B$ does not appear in the formula ($\Delta E$ is nothing but the work of the Lorentz force exerted on the particle by the electric field $E$ induced by the moving $B$). For a random distribution of moving clouds, after many interactions the particle experiences a net energy gain, because head-on collisions are more likely. This is only an average gain, hence the name stochastic acceleration, and it scales as $\beta^2$ where $\beta = v/c$, hence the name second-order Fermi acceleration (or simply Fermi II). Fermi himself realized that this process was probably not efficient enough to produce the bulk of Galactic CRs. Now if somehow only face-on collisions occur, then the energy gain is systematic, hence the name regular acceleration, and it scales as $\beta$, hence the name first-order Fermi acceleration (or simply Fermi I).

A shock wave (the S in DSA) provides such a configuration: both the upstream and the downstream medium see the opposite side arriving at the same speed $\Delta U = \frac{r-1}{r} V_{sh} (= \frac{3}{4} U_{sh}$ if $r = 4$) where $r$ is the compression ratio and $U_{sh}$ is the speed of the shock (with respect to the unperturbed upstream medium). Let us further assume that magnetic turbulence in the vicinity of the shock efficiently scatters off the particles (leading to the D in DSA), so that they are effectively isotropized on each side of the shock, meaning that their mean velocity follows the local flow velocity.<fnref n="4" /> Then, the particles experience a regular Fermi acceleration, the clouds being replaced by a reflecting wall moving at velocity $\Delta U$. Averaging Eq. (1) over all angles, one gets a mean energy gain

<footnote n="4">⁴ More precisely, the velocity that matters for particles is that of the magnetized waves present in the flow, the difference with respect to the fluid flow is called the Alfvén drift, in the case the main scattering waves are Alfvén waves, see Sects. 5.3 and 5.7.2.</footnote>
$$ \left(\frac{\Delta E}{E}\right) = \frac{4}{3} \frac{\Delta U}{c} = \frac{4}{3} \frac{r-1}{r} \frac{U_{\text{sh}}}{c} (= \beta_s \text{ if } r = 4). \qquad (2) $$

By considering the duration of a complete reflection from the opposite medium, and the probability of escape from the acceleration region, one can derive the final distribution in energy of particles. To do so, two approaches are possible: a microscopic approach, where one considers the fate of individual particles, and a kinetic approach, where one reasons with their distribution function as a function of energy or rather momentum (see details of the calculations for both approaches in Drury (1983) and references therein). This basic choice will also apply to the numerical methods presented in this review. From a general point of view an acceleration mechanism can be characterized by its acceleration time $\tau_{\text{acc}}$ (defined so that particles are accelerated at a rate $\partial E/\partial t = E/\tau_{\text{acc}}$) and its escape time $\tau_{\text{esc}}$ (defined so that particles escape the accelerator at a rate $\partial N/\partial t = N/\tau_{\text{esc}}$). If particles are injected at an energy $E_0$ with a rate $Q(E_0)$, after a time t a number density $N(E)dE = Q(E_0)\exp(-t(E)/\tau_{\text{esc}})dt$ of particles will have escaped. Now energy and time are linked by $t(E) = \tau_{\text{acc}} \ln(E/E_0)$, inserting this time in the previous relation leads to the steady-state solution

$$ N(E) \propto E^{-s} \quad \text{with} \quad s = 1 + \frac{\tau_{\text{acc}}}{\tau_{\text{esc}}}. \qquad (3) $$

In the limit where escape never occurs ($\tau_{\text{esc}} = \infty$), the hardest spectrum one can obtain is $N(E) \propto E^{-1}$. In DSA the ratio $\tau_{\text{acc}}/\tau_{\text{esc}}$ turns out to be independent of E, and so one gets a power-law distribution, of index

$$ s = \frac{r+2}{r-1} (= 2 \text{ if } r = 4). \qquad (4) $$

The spectral index is controlled by the compression ratio of the shock, and so is a universal value for strong shocks.

### 2.2.2 The transport equation and the diffusion coefficient

Assuming the particle distribution is isotropic (to first order) in momentum p, and considering here for simplicity a plane-parallel shock along direction x, in the kinetic description we may work with the quantity $f = f(x, p, t)$, which is defined so that the number density is $n(x, t) = \int f(x, p, t) 4\pi p^2 dp$, and which obeys the convection-diffusion equation (CDE)<fnref n="5" />:

$$ \frac{\partial f}{\partial t} + \frac{\partial(Uf)}{\partial x} = \frac{\partial}{\partial x} \left( D \frac{\partial f}{\partial x} \right) + \frac{1}{3} \frac{\partial p^3 f}{\partial p} \frac{\partial U}{\partial x}. \qquad (5) $$

On the right-hand side of the equation, the second term represents advection in momentum, powered by the fluid velocity divergence $\partial U/\partial x$, while the first term models the

<footnote n="5">⁵ This equation is obtained by averaging the full Fokker-Plank equation (FPE) over particle pitch-angle and azimuthal gyromotion angle (see Sect. 3.6).</footnote>
spatial diffusion of the particles, resulting from their scattering off magnetic waves, and described by a diffusion coefficient $D(x, p)$. This coefficient, together with the shock speed $U_{sh}$, sets the space- and time-scales of DSA. Upstream of the shock front, particles can counter-stream the flow up to a distance

$$ \ell_{\text{prec}}(p) = \frac{D(p)}{U_{\text{sh}}} \qquad (6) $$

which sets the scale of what is called the CR precursor region. The acceleration timescale, defined so that $dp/dt = p/t_{acc}$, goes as

$$ t_{\text{acc}}(p) \propto \frac{D(p)}{U_{\text{sh}}^2} \qquad (7) $$

where the proportionality factor is of order 8–20 depending on the shock obliquity and the Rankine–Hugoniot conditions linking up- and downstream magnetic field strengths, see Reynolds (1998). For DSA to work requires $\ell_{prec}$ to be less than the accelerator’s size, and $t_{acc}$ to be less than the accelerator’s age, which puts limits on the maximum momentum $p_{max}$ that the particles can reach (for particles that radiate efficiently like electrons, $p_{max}$ may also be limited by losses). The diffusion law is thus a critical ingredient in DSA. This aspect will however be treated in very different ways in different kinds of numerical simulations. In kinetic approaches, where one solves one equation of the kind of Eq. (5), the diffusion coefficient $D(x, p)$ or some equivalent quantity must be specified. In microscopic approaches, where one directly integrates the equation of motion of individual particles, the diffusion coefficient may actually be measured from the observed paths of an ensemble of particles.

Computing the value of the diffusion coefficient from theory is a difficult problem (see Shalchi 2009b for a review). Very generally, the diffusion coefficient can be expressed as $D = \frac{1}{3}\ell v$ where again $v$ is the particle velocity and $\ell$ its mean free path. When charged particles are deflected by Alfvén waves, $\ell$ is inversely proportional to the energy density $\delta B^2$ in waves present with the resonant wavelength $\lambda \simeq r_g$, where $r_g = \frac{pc}{qB}$ is the particle gyroradius. A special case of interest is the so-called “Bohm limit” (see e.g. Kang and Jones 1991; Berezhko and Völk 1997; Bell 2013) reached when $\ell \simeq r_g$, that is when the particles are scattered within one gyroperiod, meaning that the turbulence is random ($\delta B \sim B$) on the scale $r_g$. This constitutes a lower limit on the value of the (parallel) diffusion coefficient, and so on the acceleration time-scales. In that case, $D \propto pv$ so that

$$ D_B(p) = D_0 \frac{p^2}{\sqrt{(1 + p^2)}} \qquad (8) $$

where one can evaluate $D_0 \simeq 3 \times 10^{22}/B \text{ cm}^2 \text{s}^{-1}$ with $B$ in $\mu G$, and $p$ is expressed in $m_p c$ units.

Historically the Bohm limit has been widely favored in the literature, in its true form in Eq. (8) or with a free normalization and only keeping the “Bohm scaling” in $p$, and using the exact dependence on $p$, or only the relativistic scaling $D(p) \propto p$,
or a parametrized scaling of the form $D(p) = D_0 p^\alpha$ with free index $\alpha$ (commonly informally called “Bohm-like” coefficients). This choice of the Bohm limit may be in part due to the fact that it is the most favorable case, and also may just stem from habit given the lack of clear theoretical alternatives. Analytically it has only been derived under the assumption of strong turbulence Shalchi (2009a), and numerical studies have found that it does not generally hold (Casse et al. 2002; Candia and Roulet 2004; Parizot 2004). The validity of this assumption in the context of supernova remnant studies has been regularly questioned (Kirk and Dendy 2001; Parizot et al. 2006), and in the context of interplanetary shocks other models have been used according to the turbulence properties and shock obliquity (Dosch and Shalchi 2010; Li et al. 2012). In any case, when it comes to simulations, a key aspect is how strongly the diffusion coefficient depends on the particle’s energy, given relations (6) and (7).

### 2.2.3 The injection problem

DSA is a bottom-up acceleration mechanism, whereby (a fraction of) the particles from a plasma get boosted to very large energies. The particles are accelerated from a non-thermal distribution that extends beyond the thermal distribution of the plasma (often assumed to be a Maxwellian). The way these two populations are connected is a delicate problem. The discussion of DSA above assumes that particles are sufficiently energetic that they can leap over the shock wave and perceive it as a discontinuity, meaning that their mean free path in the magnetic turbulence is already larger than the physical width of the shock wave, which is typically of the order of a few Larmor radii of thermal ions. The way by which particles from the background plasma enter the acceleration process is known as the “injection” mechanism. The general idea, referred to as “thermal leakage” is that particles heated at the shock may be able to re-cross the shock to start the DSA cycle (Malkov and Voelk 1995; Malkov 1998). The efficiency of the injection determines the fraction of the available shock energy that is channeled into energetic particles. It is widely expected to vary as a function of parameters such as the shock obliquity, although there is no firm agreement yet on what are the most favorable configurations.

Injection is treated in much different ways according to the level of the numerical modeling. In kinetic approaches that decouple the non-thermal population (obeying Eq. 5) and the thermal population (obeying classical conservation laws), the injection process is parametrized. The simplest way to do this is to postulate that some fraction $\eta$ of the particles crossing the shock enter the acceleration process, at some momentum $p_{\text{inj}}$ above the typical thermal momentum. The requirement that the particles power-law matches the plasma Maxwellian at $p_{\text{inj}}$ actually implies that these quantities are related, as shown in Blasi et al. (2005). A more advanced approach is to use a “transparency function” to inject particles at the shock, as done by Gieseler et al. (2000). In contrast, in the Monte-Carlo simulations of the kind of Ellison and Eichler (1984) no formal distinction is made between the thermal and non-thermal populations, allowing for a more consistent treatment of injection, for a given scattering law. Only PIC simulations are able to directly address the formation of the collisionless shock concomitantly with the energization of particles, and it has been only very recently that the computational power has become sufficient to see the DSA
<figure><img src="image_9.png" /><figcaption>Fig. 3 Sketch of the physical components and their couplings in the diffusive shock acceleration mechanism</figcaption></figure>

power-law naturally emerge (e.g., Caprioli and Spitkovsky 2014a)—although still on
very small space- and time-scales compared to any astrophysical object of interest.
The results associated to these simulations are described with more details in Sect. 4.
This again illustrates the need for a model at different scales and their entanglement,
a given numerical approach often relying on results obtained from other approaches
for the aspects it cannot describe.

2.2.4 Back-reaction and non-linear effects

DSA at astrophysical shocks involves three kinds of actors: energetic particles, a
plasma flow, and magnetic waves (see Fig. 3). Charged particles are being injected
from the plasma and accelerated at the shock, thanks to their confinement by magnetic
turbulence. In our discussion so far we have assumed a prescribed background plasma
flow and magnetic turbulence, that is, we were implicitly discussing the “test particle”
regime. But if the acceleration process is efficient (meaning that a substantial fraction
of the available energy ends up into particles), then the particles will play a role in the
dynamics of the plasma and in the evolution of the magnetic field. This in turn will
affect the way they are being accelerated, so that the DSA process becomes non-linear
(NLDSA). The time-dependent problem is intractable analytically in the general case,
which is the reason why studies of (efficient) particle acceleration rely on numerical
techniques, as described in this review. To end this section on DSA, we summarize
the main aspects of the two back-reaction loops: of the particles on the plasma flow,
and of the particles on the magnetic turbulence.

Even before DSA theory was established, Parker (1958b) noted that CRs modify the
medium in which they propagate: being relativistic they lower the overall adiabatic
index of the flow. We can make this more precise by considering the CR pressure,
defined as
$$P_{cr} = \int_{p} \frac{p v}{3} f(p) 4\pi p^2 dp = \frac{4\pi}{3} m_{p} c^2 \int_{p} \frac{p^4}{\sqrt{1 + p^2}} f(p) dp \quad (9)$$

(where in the right expression momenta are expressed in $m_{p}c$ units). This quantity can grow without limit in the linear regime. But as CRs diffuse upstream of the shock in an energy dependent way, a gradient of CR pressure forms in the precursor, which produces a force acting on the plasma. This CR-induced force pre-accelerates the plasma ahead of the shock front, leading to the formation of a smooth, spatially extended velocity ramp upstream of the shock front. The shock itself is thus progressively reduced to a so-called “subshock”, whose compression ratio is $r_{sub} < 4$, while the overall compression ratio $r_{tot}$ (measured from far upstream to far downstream) becomes $> 4$ from mass conservation—the plasma is more compressible when particles become relativistic, and even more so when the particles escape the system (Berezhko and Ellison 1999). As particle of different energies can explore regions of different extent ahead of the shock, they will feel different velocity jumps, so that the spectral slope defined by Eq. 4 becomes energy-dependent. The spectrum is thus no longer the canonical power-law, but gets concave. Particles of low energy ($p \ll m_{p}c$) only sample the sub-shock, feeling a compression $r_{sub}$ that produces a slope larger than the canonical value of 2 (that is, a steeper spectrum), whereas particles at the highest energies sample the whole shock structure, feeling a compression $r_{tot}$ that produces a slope that is smaller (that is, a flatter spectrum). So one of the historically most attractive features of DSA—its ability to naturally produce power-law spectra—cannot strictly hold when it is efficient.

Turning to magnetic turbulence, it was also observed early on (e.g., Skilling 1975b) that the CRs can generate themselves the waves that will scatter them off. They can indeed trigger various instabilities by streaming upstream of the shock, which generates magnetic turbulence, which is then advected to the shock front and the downstream region. Denoting by $W(\mathbf{k}, x, t)$ the power spectrum of the magnetic waves where $\mathbf{k}$ is the wave vector, its evolution obeys a transport equation of the general form (here written for simplicity along one dimension). Assuming for simplicity that $U$ is constant we have:

$$\frac{\partial W}{\partial t} + U \frac{\partial W}{\partial x} = \Gamma_g - \Gamma_d, \quad (10)$$

where $\Gamma_g$ is the growth rate of the waves, which is dictated by the particles, and $\Gamma_d$ is their damping rate in the plasma. For a CR-induced streaming instability the growth rate $\Gamma_g$ scales as the gradient in the CR pressure. Using hybrid MHD+particle simulations, Lucek and Bell (2000), Bell (2004) showed that the seed magnetic field can be amplified by up to two orders of magnitude, which is important because in turn the magnetic field controls the diffusion and confinement of the particles and thus the maximum energy that they can reach (Bell and Lucek 2001). This discovery prompted a slew of work on this complex topic (see Schure et al. 2012 for a review). Two regimes of streaming instabilities can operate, a resonant instability at work when the Larmor radius of a particle matches the wavelength of a perturbation, and a non-resonant instability driven by the current of particles (see also Gary 1991). The amplified field saturates at an energy density $\delta\mathbf{B}^2$ that scales as $U_{\text{sh}}^2$, or possibly even $U_{\text{sh}}^3$ (Vink 2012). Other longer wavelength instabilities can be triggered too. Combined with
<figure><img src="image_10.png" /><figcaption>**Fig. 4** Particle trajectories along the shock front due to SDA or SSA mechanisms. The flow is directed along the *x*-direction, and particles are injected at the left of the box. The magnetic field is perpendicular to the plane of the plot, as revealed by the gyro-motion of the particles. In the SDA process, the particle is forced to cross the shock several times. In the SSA process, the particle moves along the shock front. Image reproduced with permission from Shapiro and Üçer (2003). Copyright by Elsevier</figcaption></figure>

observational evidence for high magnetic fields at the shocks of young supernova remnants (Reynolds et al. 2012), this has lead to the current view that magnetic field amplification (MFA) is a critical ingredient of DSA. This effect has been integrated in numerical simulations in different ways depending on the level of description of the particles and waves. For instance, in the kinetic description of DSA, one can in principle compute the diffusion coefficient $D(x, p, t)$ self-consistently from the cosmic-rays distribution $f(x, p, t)$. Obtaining a fully consistent description of the dynamical evolution of the particles, the plasma flow, and the magnetic turbulence, is still a work in progress.

## 2.3 Shock drift and shock surfing acceleration

These two acceleration processes rely on the effect of the convective electric field $E = -U/c \times B$ induced by magnetized fluid motions towards the shock. The difference between the two processes results from the way the particles are either confined at the shock front in the case of shock surfing or move up- and downstream in the case of shock drift (Hudson 1965). Figure 4 shows different trajectories adopted by particles due to these two different processes.

### 2.3.1 Shock drift acceleration: SDA

Acceleration associated to the drift of the particle's guiding center depends strongly on whether the shock is super- or subluminal. The super- or subluminal character of a shock depends on the speed of the intersection point of the upstream magnetic field with the shock front: in superluminal shocks this speed is larger than the speed of light, in subluminal shocks it is smaller.
In subluminal shocks it is always possible to find a frame where the convective electric field vanishes (so where the fluid velocity lies along the magnetic field line direction), this is the so-called de Hoffmann-Teller (HT) frame (Kirk 1994). In the HT frame particle energy is conserved. As an upstream field line intersects the shock the particle guiding center drifts along the shock and can either be transmitted or reflected at the shock front because the magnetic field is compressed there. The energy gain is the highest for particles reflected at the shock front (Decker 1988). This effect is similar to a reflection at the edge of a magnetic bottle. A calculation assuming that adiabatic theory applies uses a Lorentz transformation between a frame at which the shock is stationary to the HT frame to derive the energy gained by a particle reflected at the shock. Averaging over initial particle pitch-angles gives a ratio of the particle energy after the reflection to the initial energy of

$$ \left\langle \frac{E_f}{E_i} \right\rangle = \frac{1 + \sqrt{1 - b}}{1 - \sqrt{1 - b}}, \qquad (11) $$

where $b = B_u/B_d$ is the ratio of the upstream to downstream magnetic field strengths. For $b = 0.25$ (for a compression of the magnetic field by a factor 4) we find a maximum ratio $\langle E_f/E_i \rangle_{\max} \simeq 15.8$ and only half of this for a transmitted particle.

In superluminal shocks, a configuration obtained in perpendicular shocks, the adiabatic invariant $p_{\perp}^2/B$ is conserved even while the particle crosses the shock, as long as its speed is much larger than the shock speed (Webb et al. 1983; Whipple et al. 1986). The maximum value of the ratio $E_f/E_i$ is $1/\sqrt{b} = 2$. Superluminal shocks are the rule in relativistic flows.<fnref n="6" /> Begelman and Kirk (1990) investigate the condition under which SDA process operates at relativistic perpendicular shocks associated with the synchrotron emission of radio galaxy hot spots. As the flow speed gets closer to the speed of light the condition for the adiabatic theory breaks down. Begelman and Kirk (1990) propose an alternative method by following individual particle orbits. Due to shock dynamics, particles can cross the shock front at maximum three times before being advected downstream.

### 2.3.2 Shock surfing acceleration: SSA

Shock surfing is produced when a particle is trapped between the shock electrostatic potential $e\phi$ (+e is the particle charge in case of a proton) which appears at the shock ramp and the upstream Lorentz force along the shock normal which carries the particle back to the front. Original ideas about this mechanism can be found in Sagdeev (1966), Sagdeev and Shapiro (1973). The particle is accelerated under the action of the convective electric field until its kinetic energy along the shock normal exceeds $e\phi$. The trapped particle accelerates essentially along the shock front like a surfer's motion along a wave; hence Katsouleas and Dawson (1983) name this process *shock surfing* acceleration. The SSA mechanism is often invoked at quasi-perpendicular shocks as a pre-acceleration process. It allows to inject particles beyond the energy threshold

<footnote n="6">⁶ The shock is subluminal only if the obliquity angle between the upstream magnetic field and the shock normal is $\theta_u \le 1/\Gamma_{sh}$, where $\Gamma_{sh}$ is the shock Lorentz factor.</footnote>
for DSA to operate (Zank et al. 1996; Lee et al. 1996). Particles gain energy as long as they stay trapped at the shock front. The acceleration ceases if they escape either upstream if the magnetic field has some obliquity or downstream if the Lorentz force exceeds the electrostatic force at the shock layer. Particle acceleration can operate also in the relativistic regime where particles with small initial speeds are trapped the longest. Shapiro and Üçer (2003) for instance find that particles can have 10 bounces and reach Lorentz factors ~ 2 through this mechanism.

## 2.4 Fermi acceleration process at relativistic shocks

We now discuss how particle acceleration is affected when the flow itself is relativistic.

### 2.4.1 General statements

If we consider a relativistic shock front moving with a Lorentz factor $\Gamma_{sh} = 1/\sqrt{(1 - (U_{sh}/c)^2)}$, the relative energy gain as the particle is doing a shock crossing cycle (e.g., up-down-upstream) can be obtained from relativistic kinematics by imposing a double Lorentz transformation between the upstream and downstream rest frames. The relative energy variation between the final $E_f$ and the initial $E_i$ particle energies is (Gallant and Achterberg 1999)

$$ \frac{\Delta E}{E} = \frac{(E_f - E_i)}{E_i} = \Gamma_r^2 (1 - \beta_r \mu_{u \to d}) (1 + \beta_r \mu'_{d \to u}) - 1, \quad (12) $$

where unprimed and primed quantities mark upstream and downstream rest frame quantities respectively. The relative Lorentz factor between upstream and downstream is $\Gamma_r = \Gamma_{sh}/\sqrt{2}$. The mean energy gain is obtained by averaging over the cosine $\mu_{u \to d}$ and $\mu'_{d \to u}$ of the penetration angles of the particle from upstream to downstream and downstream to upstream with respect to the direction of the boost. In the most optimistic case a high relative energy gain $\Delta E/E \sim \Gamma_{sh}^2$ can be achieved (Vietri 1995). However, this gain is restricted to the first cycle if the initial particle distribution is isotropic. The particle distribution after one cycle becomes highly anisotropic, beamed in a cone of size $1/\Gamma_{sh}$ and due to particle kinematics, the average relative gain drops to ~ 2 for the next crossings (Gallant and Achterberg 1999). Particle deflection in the cone can either proceed through its motion in a uniform magnetic field in the absence of scattering waves or by scattering with resonant waves with $kr_g \sim 1$. Resonant scattering occurs only if the amplitude of the magnetic perturbations is small enough (Achterberg et al. 2001). The shock particle distribution in the test-particle limit shows a universal energy spectrum $N(E) \propto E^{-2.2}$ whatever the deflection upstream if scattering is effective downstream. This result is consistent with the index of the relativistic electron distribution producing synchrotron radiation in GRB afterglow (Waxman 1997). Note that this result has been obtained for an isotropic turbulence downstream. A more general formulation in terms of shock speed gives an energy index (Keshet and Waxman 2005) of
<figure><img src="image_11.png" /><figcaption>Fig. 5 Spectral index of the shock particle distribution as function of the shock Lorentz factor. Solid lines are the solutions obtained in Keshet and Waxman (2005) compared with the solutions of Kirk et al. (2000) plotted with symbols. A strong shock solutions with the Jüttner/Synge equation of state is shown in solid line and crosses, a strong shock with fixed adiabatic index $\gamma_{ad} = 4/3$ is shown in dashed line and x-marks, and a shock with a relativistic gas where $\beta_u\beta_d = 1/3$ is shown in dash-dotted line and circles. Image reproduced from Sironi et al. (2015). Copyright by Springer</figcaption></figure>

$$s = \frac{\beta_u - 2\beta_u\beta_d^2 + \beta_d^3 + 2\beta_d}{(\beta_u - \beta_d)}, \quad (13)$$

where $\beta_{u/d}$ is the upstream/downstream fluid velocity normalized to c; $s = 2+2/9 \simeq 2.2$ is recovered in the ultra-relativistic limit ($\beta_u \to 1$ and $\beta_d \to 1/3$). The value of the ultra-relativistic index has also been assessed by numerical simulations using a Monte-Carlo method (Bednarz and Ostrowski 1998; Achterberg et al. 2001) or a semi-analytical method based on the derivation of eigenfunctions in the particle pitch-angle cosine of the solution of the diffusion-convection equation (Kirk et al. 2000). However, this spectrum is not properly universal in the sense that the index depends on the geometry of the turbulence (Lemoine and Revenu 2006).

Figure 5 shows the index of the shock downstream particle distribution as a function of the shock Lorentz factor from mildly relativistic to ultra-relativistic regimes for different equation of state of the relativistic gas (the quantity plotted is $s+2$ with our notation).

It appears that unless some particular turbulence develops around the shock front, Fermi acceleration associated with repeated shock crossings does not operate at relativistic shocks because of particle kinematic condition to cross the shock front (Begelman and Kirk 1990; Lemoine et al. 2006; Pelletier et al. 2009). This fact results from that relativistic shocks are generically perpendicular unless the magnetic field upstream is oriented within an angle $1/\Gamma_{sh}$ along the shock normal. If the background turbulence around the shock is absent or if its coherence length is larger than the particle's gyroradius it can be shown that while returning to the shock from downstream to upstream the particle is unable to do more than one cycle and a half before
being advected downstream (Lemoine et al. 2006). There are two necessary condi-
tions for an efficient scattering and particle acceleration: (1) the turbulence which
develops around the shock has to be strong, with a perturbed magnetic field such that
$\delta B/B_0 > 1$, where $B_0$ is the background upstream magnetic field strength, (2) the
perturbations have to be at scales smaller than the particle gyroradius (with respect
to the total magnetic field) (Pelletier et al. 2009). One drawback of these conditions
is that, as the micro-turbulence develops with a coherence scale smaller than $r_g$, the
spatial diffusion coefficient scales as $D \propto r_g^2$, and as the acceleration timescale scales
as $t_{\text{acc}} \propto D/U_{\text{sh}}^2 \propto E^2$ even when $U_{\text{sh}} \to c$, the time required to reach extremely high
energies can become very large. This is the main reason which explains why ultra-
relativistic GRB shocks $\Gamma_{\text{sh}} \gg 1$ cannot be at the origin of CRs with energies $\sim 100$
EeV (Lemoine and Pelletier 2010; Plotnikov et al. 2013; Reville and Bell 2014).

**One note on the origin of the micro-turbulence** We have seen that the onset of micro-turbulence is necessary for the Fermi process to operate at relativistic shocks. What type of micro-turbulence develops depending on shock velocity regimes and upstream medium properties?

The characteristic scale on which the micro-turbulence can develop is given by the CR precursor scale $\ell_{\text{prec}}$. It is either set by the regular gyration of particles reflected by the shock front, in which case $\ell_{\text{prec}} = r_g / \Gamma_{\text{sh}}^2$, or by diffusion, in which case $\ell_{\text{prec}} = r_g^2 / \ell_c$, where $\ell_c$ is the turbulence coherence scale. As the energy spectrum is softer than $E^{-2}$, low energy particles carry the bulk of free energy and generate most of magnetic perturbations. Hence, the scale over which micro-instabilities can develop in relativistic shocks is small. This restricts the number of instabilities to a few (Lemoine and Pelletier 2010). The nature of the dominant instability depends on two main parameters: the shock Lorentz factor $\Gamma_{\text{sh}}$ (or momentum $\Gamma_{\text{sh}}\beta_{\text{sh}}$ for mildly relativistic shocks), and the upstream magnetization $\sigma_u = B_0^2 / (4\pi(\Gamma_{\text{sh}}(\Gamma_{\text{sh}}-1)n^2mc^2))$, where $n^*$ is the ambient (upstream) proper gas density composed of particles of mass $m$ (Sironi et al. 2015; Marcowith et al. 2016). In weakly magnetized shocks ($\sigma_u \le 10^{-3}$) the dominant instability is the electromagnetic filamentation or Weibel instability (Weibel 1959). Filamentation/Weibel instabilities grow due to the presence of two counter-streaming population of particles and produce modes in the direction perpendicular to the streaming direction (Fried 1959; Bret 2009). Plotnikov et al. (2013) discuss also the case of the oblique two-stream instability which can have a competitive growth rate with respect to the filamentation instability. Finally, the Buneman instability has been discussed as a source of electron heating in relativistic shock precursors (Lemoine and Pelletier 2011). At higher magnetization ($0.1 > \sigma_u > 10^{-3}$) a current driven instability either in subluminal (Bell 2004, 2005; Reville et al. 2006; Milosavljević and Nakar 2006) or superluminal shocks (Pelletier et al. 2009; Casse et al. 2013; Lemoine et al. 2014) can develop in the precursor. If $\sigma_u > 0.1$ then the gyration of particles in the background magnetic field gains in coherence and the shock is mediated by the synchrotron maser instability. This instability produces a train of semi-coherent large amplitude electromagnetic waves that escapes into the upstream medium (Gallant et al. 1992; Hoshino et al. 1992; Plotnikov and Sironi 2019). The interaction of this wave with the background plasma is a source of an efficient electron
pre-heating up to equipartition with protons (Lyubarsky 2006; Sironi and Spitkovsky 2011). The question of whether this wave can lead to a significant particle acceleration is debated (e.g., Hoshino 2008; Iwamoto et al. 2017; Lyubarsky 2018).

## 2.4.2 Progress with fully kinetic simulations

Until late 2000s, most progress on the understanding particle acceleration at relativistic shocks were supported by Monte-Carlo simulations (e.g., Bednarz and Ostrowski 1998; Lemoine and Pelletier 2003; Ellison and Double 2004; Lemoine and Revenu 2006) or semi-analytic approaches (Kirk et al. 2000; Achterberg et al. 2001; Keshet and Waxman 2005). As pointed out by Bykov and Treumann (2011), recent advances were possible by employing fully kinetic PIC simulations. Self-consistent build-up of Fermi process was observed and a survey in which parameter space it operates was done (Spitkovsky 2008b; Martins et al. 2009; Sironi and Spitkovsky 2009, 2011; Sironi et al. 2013; Plotnikov et al. 2018). For instance, these simulations demonstrate that the Weibel-filamentation instability dominates in controlling the shock structure in weakly magnetized shocks, as predicted by Medvedev and Loeb (1999) and Gruzinov and Waxman (1999). Particle acceleration is correlated to the efficiency of triggering this instability. Typically, the non-thermal particles contain about 1% of total particle number and about 10% of total energy. The tail develops into a power-law with spectral slope $s \simeq 2.4$ that is close to the semi-analytic prediction of 2.2. The maximum energy of particles evolves in time as $E_{max} \propto \sqrt{t}$ (Sironi et al. 2013) due to the small-scale nature of the magnetic turbulence (see above). In the small-angle scattering regime, the spatial diffusion coefficient of particles is $D \propto E^2$, unless the external magnetic field imposes a saturation that sets the maximal particle energy to be $E_{max} > e\delta B^2 / B_0 l_c$ (Plotnikov et al. 2011; Marcowith et al. 2016). For typical parameters in ultra-relativistic shocks with $\Gamma_{sh} \sim 100$ propagating in the ISM with $B_0 \sim 3 \mu G$, this energy does not exceed $10^{16}$ eV. Section 4 presents more detailed discussions of these studies.

## 2.4.3 Long term evolution

The fate of the micro turbulence and, more generally, the long-term evolution of the weakly magnetized shocks remains the major unanswered question in relativistic (but also in non-relativistic) shock physics. As this micro turbulence is composed of initially short wavelength perturbations, these are expected to be rapidly damped by Landau damping downstream (Gruzinov 2001; Chang et al. 2008; Lemoine 2015). One possibility to overcome this effect would be to have some amounts of inverse cascade to generate large wavelengths or to rely one the perturbations generated upstream by high-energy particles and transmitted downstream. These possibilities, and others are detailed in Sect. 3.2 of Sironi et al. (2015).

# 2.5 Reconnection in astrophysical flows

Magnetic reconnection can occur in a collisional or in a collisionless plasma. The bulk part of the particles are accelerated to order of the Alvén speed and heated in the

<figure><img src="image_12.png" /></figure>
same process. A fraction of the particles can be accelerated to much higher velocities and form a power-law up to very high Lorentz factors. The power-law slopes can be harder than power-laws produced by a Fermi process, e.g. in collisionless shocks. In astrophysics, reconnection (REC) is a highly important process to accelerate particles.

In this introductory section, we briefly summarize some important concepts of magnetic reconnection. For any further details, we refer to recent reviews on the subject by Zweibel and Yamada (2009), Zharkova et al. (2011) (solar flares), Melzani (2014), Gonzalez and Parker (2016), Jafari and Vishniac (2018). A more deep insight into subjects directly related to this review, the acceleration of particles, will be given in Sect. 4.2.

**Astrophysical objects where reconnection takes place** Well-known REC sites in the solar systems are the upper chromosphere and the corona of the Sun as well as the magnetotail and the magnetopause of planets. There, predominately electrons are accelerated to very high speeds. REC may be partly responsible for heating the solar corona and thus for the existence of the solar wind. But highly accelerated particles are also a severe threat for spacecrafts and astronauts and even aircraft passengers. They are at the source of geo-magnetic storms which severely endanger communication and power grids on Earth. The demand for a better understanding of space-weather is one reason why in recent years the effort to understand REC has intensified and brought decisive new insights.

In outer space, REC was found to play a crucial role in the understanding of high-energy objects such as pulsars and their winds and nebulae, as well as magnetars, (micro-)quasars, and GRBs. In most of these objects, REC is partly a driver of their dynamics. Besides shocks and wave-turbulence, REC can accelerate particles to highest energies under such conditions. These particles and their interaction with the environment also inevitably contribute to the emission spectrum of such objects. In addition, they may be at the source of the production of high-energy neutrinos observed on Earth. REC in such objects is mostly relativistic in that the energy stored in the associated magnetic fields exceeds the rest mass energies of electrons and protons.

Figure 6 shows the parameter space of magnetic REC present in astrophysics and relates it to concrete systems. As can be taken from the figure, collisional and non-collisional REC equally contribute to the overall picture of magnetic REC in space. Indicated are as well different other regimes which will be discussed below. Fusion devices like ITER and TFTR and experimental setups like MRX, NGRX, MST, VTF are also shown.

## 2.5.1 Collisional reconnection models

**Sweet-Parker model** The first theory of magnetic reconnection was presented by (Sweet 1958) for a *collisional* plasma with resistivity $\eta$; $J = \eta E$. Parker (1957) worked out the scaling relations presented below. It was soon clear that it is not always applicable, as this model predicts too slow events as compared to observations. The question how to make REC fast is still today in the center of the discussion and not generally solved.
<figure><img src="image_13.png" /><figcaption>Fig. 6 Different REC regimes, as derived by Daughton and Roytershteyn (2012) (left panel) and Ji and Daughton (2011) (right panel). $S_L$ is the Lundquist number defined by Eq. (17), $L_{sp}$ the macroscopic system size, $\rho_i$ the ion Larmor radius in the asymptotic magnetic field (including the guide field), not to be confused with the ion mass density. The red curve is computed for $\beta_p = 0.2$ and a REC rate of $R = 0.05$. Images reproduced with permission from Daughton and Roytershteyn (2012), copyright by Springer; and from Ji and Daughton (2011). Copyright by AIP</figcaption></figure>

Assume a collisional plasma with a certain resistivity. Then the induction equation and Ohm's law become, with $\mathbf{U}$ the flow velocity,

$$ \frac{\partial \mathbf{B}}{\partial t} = \nabla \wedge (\mathbf{U} \wedge \mathbf{B}) + \eta \nabla^2 \mathbf{B}; \quad \mathbf{E} + \frac{\mathbf{U}}{c} \wedge \mathbf{B} = \frac{\mathbf{J}}{\sigma}. \qquad (14) $$

The non-ideal terms are the resistive diffusivity $\eta = c^2/4\pi\sigma$ in the induction equation and the resistive current in Ohm's law expressed in terms of the electrical conductivity $\sigma$. The non-ideal induction equation makes clear that there is a competition between the diffusion of the magnetic field (governed by the resistive time-scale $\tau_R$) and the ideal evolution (governed by the Alfvénic time scale $\tau_A$). This balance is expressed by the magnetic Reynold's number, $R_m \equiv \ell \cdot \tilde{U}/\eta$ with $\ell$ a characteristic length scale and $\tilde{U}$ a typical velocity of the system. In ideal MHD $R_m \gg 1$ and REC is suppressed; whenever $R_m \ll 1$ field diffusivity wins and REC becomes possible though not mandatory.

Figure 7 describes a 2D steady situation. Plasma from an outer ideal region flows in parallel to the x-direction towards a dissipation region, which has a length-scale, $L$, and a thickness, $\delta$. The inflow velocity is just given by the $E \times B$-drift in the plasma. In the outer ideal region ($R_m \gg 1$), the plasma is frozen to the magnetic flux. This is no longer true in the diffusion region where the resistivity is dominant: the plasma decouples from the magnetic field. This opens the possibility that the field reconnects and plasma is expelled in the z-direction. These outflows are called exhausts.

Applying mass and energy conservation, non-compressibility, and that the field energy is dominant at the inflow and the kinetic energy of the particles at the outflow, two important relations follow:
<figure><img src="image_14.png" /><figcaption>Fig. 7 Sweet-Parker model of magnetic reconnection. See details in the text. Image adapted from Melzani (2014)</figcaption></figure>

$$U_{\text{out}} = \sqrt{2} \frac{B_0}{\sqrt{4\pi m n_{\text{in}}}} = \sqrt{2} U_{\text{A,in}}, \qquad (15)$$

$$\frac{\delta}{L} = \frac{U_{\text{in}}}{U_{\text{A,in}}} = M_{\text{a,in}}, \qquad (16)$$

where $U_{A,in}$ is the Alfén speed and $M_{a,in}$ the Alfénic Mach-number of the inflow. Thus, the outflow speed in the exhausts is of order of the Alfénic speed of the inflow. Assuming non-forced REC, the inflow is just $E \times B$-drift, $U_{in} = cE_y/B_0 = \eta/\delta$ and thus generically very sub-Alfénic. The Lundquist-number, $S_L$, and the REC rate $R$ are defined as

$$S_L \equiv \frac{L U_{\text{A,in}}}{\eta} \sim \left(\frac{L}{\delta}\right)^2 \sim \left(\frac{U_{\text{A,in}}}{U_{\text{in}}}\right)^2 \sim (M_{\text{a,in}})^2, \qquad (17)$$

$$R \equiv \frac{U_{\text{in}}}{U_{\text{out}}} \sim \frac{\delta}{L} \sim 1/S_L^{1/2}. \qquad (18)$$

The Lundquist-number is equal to the magnetic Reynolds-number, $R_m$, for the case where the typical velocity is equal the Alfvénic speed of the inflow. Highly conducting plasmas as found in astrophysics have high Lundquist numbers: laboratory plasma experiments typically have Lundquist numbers between $10^2$ and $10^8$. In astrophysics, they are higher, up to $10^{20}$ (Fig. 6, right panel) and thus Sweet-Parker reconnection rates very low—in contrast with what is observed.

The ratio between the incoming and outgoing energy flux in Sweet-Parker reconnection is $\propto 1/S_L$. So, energy is indeed transferred from the magnetic field (incoming flux) to particles (outgoing flux). In the diffusive region with de-magnetized particles, the electric field can accelerate them. But there are other, most probably even more important acceleration mechanisms as will be discussed in Sect. 4.2.

REC rates of such high Lundquist numbers are much too low as compared to what is observed. This result can be translated to time-scales. The magnetic Reynolds-number can also be expressed as $R_m = \tau_R/\tau_A$, where $\tau_R$ is the resistive diffusion time-scale and
<figure><img src="image_15.png" /><figcaption>**Fig. 8** Illustration of the fast REC Petschek model. Image adapted from Melzani (2014)</figcaption></figure>

$\tau_A$ Alfvén time-scale. In astrophysical environments with high Lundquist numbers it is thus found that $\tau_R \ll \tau_A$. The typical Sweet–Parker REC time-scale is $\sim \tau_A^{1/2} \tau_R^{1/2}$. This indicates that Sweet–Parker REC is indeed faster than resistive diffusion of the magnetic field (scaling with $\tau_R$). However, it is much too slow when compared to REC observed in astrophysical plasmas which scales as 10–100 $\tau_A$. For instance, in flares of the solar corona, $S_L \sim 10^8$ (see Fig. 6), $U_A \sim 100$ km s$^{-1}$, and $L \sim 10^4$ km. Thus, the Sweet–Parker-timescale is a few tens of days. Observed is a magnetic energy release within a few minutes to an hour. This major discrepancy is known as the fastness problem of Sweet–Parker REC.

On the other hand, numerical simulations based on resistive MHD as well as experiments such as the MRX, the Magnetic Reconnection Experiment (Ji et al. 1998) are in good agreement with the Sweet–Parker model. Clearly, the Sweet–Parker model has its deficits, in that it neglects dimensionality and any time-dependence, as well as viscosity, compressibility, downstream pressure, and, in particular, turbulence and is strictly valid only for a collisional plasma. There have been numerous papers addressing the fastness problem. Only in recent years, significant progress in this question has been made, see below.

**Petschek model** Petschek (1964) proposed a REC model in which the reconnection rate is nearly independent of the Lundquist number, $v_{in}/v_{A,in} \approx \pi/8 \ln S_L$ and thus REC is fast. The trick is to add, in close neighborhood to the separatrices slow shocks (left panel of Fig. 8) into the configuration. In this way, particles can be accelerated without having to pass through the inner dissipation region with resistive dissipation. Instead, magnetic energy can be conversed to kinetic particle energy in the shocks.

However, all resistive MHD simulation are in agreement with the Sweet–Parker-model unless a localized anomalously large resistivity is used, mimicking that the mean free particle path becomes larger than the reconnection layer. Otherwise, shocks

<figure><img src="image_16.png" /></figure>
<figure><img src="image_17.png" /><figcaption>**Fig. 9** Illustration of the fast REC stochastic model by Lazarian and Vishniac (1999). Image adapted from Melzani (2014)</figcaption></figure>

are not observed in MHD simulations. Therefore, the Petschek-model is likely not a model of resistive MHD—though this is still a controversial question. However, recent PIC simulations in a collisionless plasma show an X-point and separatrix-structure in reconnection, which resembles somehow the Petschek-model (Higashimori and Hoshino 2012; Liu et al. 2012; Lapenta et al. 2015)—at least at scales larger than kinetic scales. There is also observational evidence that in the reconnection region of Earth’s magnetotail, slow shocks are present (Eriksson et al. 2004). This discussion will be resumed in Sect. 4.2.

**Turbulence** External or self-generated in the REC process—seems to be the key process which allows resistive MHD REC to be fast. As indicated in Fig. 9 turbulent fluctuations allow to form many, much smaller scaled, reconnection spots along the global length, L, of the sheet. As shown in Lazarian and Vishniac (1999), the REC becomes thus much faster and is independent of the exact REC mechanism at each of these spots (Sweet–Parker, collisionless, ...). The exact result depends, however, on the nature of the turbulence and its fluctuation. Numerical simulations show good agreement with the analytic result (Kowal et al. 2009, 2012).

Simulations show that a Sweet–Parker like current sheet generates islands above a critical Lundquist number $S_c \sim 10^4$ (Daughton and Roytershteyn 2012). In relativistic flows, this critical number may be higher, $S_c \sim 10^8$ (Zanotti and Dumbser 2011). This is confirmed by newer investigations and linked to an extremely fast growing tearing instability of the current sheet (Del Zanna et al. 2016; Papini et al. 2018). This limit is indicated as the green line in the left panel of Fig. 6.
Lapenta (2008) showed that a Sweet-Parker sheet setup in a Harris or force-free equilibrium sheet develops slow REC. On a much longer time-scale, tearing modes start to fragment the sheet and several X-points form. The exhausts of these X-points generate turbulence leading to multiple short lived REC regions, popping up randomly, frequently and at multiple locations simultaneously. Consequently, fast REC sets in. Similar findings for 3D resistive reconnection are presented by Oishi et al. (2015). By linking such self-generated turbulence with external turbulence, Lapenta and Lazarian (2012) formulate a united approach. So one may, with still some care, conclude that also collisional, resistive REC is fast, at least under certain conditions.

## 2.5.2 Collisionless reconnection

On length scales shorter than the ion inertial length $c/\omega_{p,i}$ where $\omega_{pi} \equiv \sqrt{4\pi n_i Z^2 e^2/m_i}$ is the ion plasma frequency, ions decouple from electrons and the magnetic field becomes frozen into the electron fluid rather than the bulk plasma. Consequently, other terms than just resistivity start to contribute to the Ohm's-law. For instance, based on a two-fluid non-relativistic plasma model, Melzani (2014) derives a more complex Ohm's-law for electrons:

$$ \underbrace{\mathbf{E} + \frac{\mathbf{v}_i}{c} \wedge \mathbf{B}}_{\text{E-field in the ion plasma frame}} = \underbrace{\frac{1}{n_e e} \mathbf{J} \wedge \mathbf{B}}_{\text{Hall term}} - \underbrace{\frac{m_e}{e} \left( \frac{\partial \mathbf{v}_e}{\partial t} + \mathbf{v}_e \cdot \nabla \mathbf{v}_e \right)}_{\text{electron bulk inertia}} - \underbrace{\frac{1}{n_e e} \nabla \cdot \mathbf{P}_e}_{\text{e}^{-} \text{thermal inertia}} + \underbrace{\frac{\chi}{(n_e e)^2} \mathbf{J}}_{\text{e-i collisions}} + \underbrace{\frac{\chi_e}{n_e e} \nabla^2 \mathbf{v}_e}_{\text{e-e collisions}} \quad (19) $$

Here, $n_e$ is the electron number density, $v_i$, $v_e$ the ion and the electron velocity respectively, c the speed of light, e the elementary charge, $\chi$ accounts for the effect of collisions between electrons and ions which, in general can be anisotropic and depend on the magnetic field orientation. $\chi_e \nabla^2 \mathbf{v}_e$ describes the electron viscosity due to electron-electron collisions. $P_e$ is the pressure tensor

$$ P_e = \int d^3\mathbf{v} m_e (v_a - \bar{v}_a)(v_b - \bar{v}_b) \quad (20) $$

with $a, b = x, y, z$, and $\bar{v}$ is the mean velocity. Electron inertia, both thermal and bulk, now contribute to the non-ideal terms. In particular, if the plasma is completely collisionless, ($\chi = \chi_e = 0$), these are the only contribution of non-idealness of the plasma.

The sketch in the left panel of Fig. 10 shows that the dissipation region now is subdivided into a larger ion dissipation region with a size of $\delta_i$ and a smaller electron dissipation region, sized to $\delta_e$. Here, $\delta_{i,e}$ denotes the ion and electron inertial length.

On these scales, the Hall effect becomes important, because now the magnetic field lines are advected with the electrons while the ions no longer follow this motion. The Hall term is not responsible for REC as it appears when the magnetic flux is
<figure><img src="image_18.png" /><figcaption>**Fig. 10** Collisionless REC in a electron-ion plasma. Left panel: sketch of REC on a scale smaller than the ion inertial length. The diffusion regions for ions are much larger than that for electrons. The current sheet is more like an X-point than a double y-point. Ion trajectories normally do not pass through the electron non-ideal region. Image adapted from Melzani (2014). Right panel: X-point, exhausts and islands from a collisionless electron-ion PIC-simulation of REC using a mass-ratio $m_i/m_e = 25$. Visibly, the ion diffusion regions is about a factor of 5 ($\delta_k \sim \sqrt{m_k}$, $k = e, i$) larger than the electron diffusion region. Image adapted from Melzani et al. (2014a)</figcaption></figure>

still frozen to the motion of electrons. However, there is a debate whether it may contribute to the fastness of REC as it allows to accelerate electrons to higher speeds, increasing the bulk inertia. As can be taken from the right panel of Fig. 10 this two-layer picture derived from a two-fluid model is quite accurately reproduced by full kinetic simulations though the two-fluid model will not provide the full picture as effects like wave turbulence, Landau-damping and particle acceleration to speeds much higher than the Alfvén speed. Fully collisionless REC is found to be always fast. It will be further addressed in Sect. 4.2.

### 2.5.3 Other effects

**Dimensionality** REC in 3D shows a variety of new features. In some cases still a separatrix-like reconnection as in 2D can be observed, but there are also many other cases. It is not the place to discuss this here in detail. A good summary can be found in Melzani (2014).

**Fat tails and high-energy power-laws** Magnetic reconnection is efficient to accelerate particles, both in the collisional and collisionless regime. The typical speed of accelerated particles is the local Alfvén speed. If the flow is highly magnetized this speed can be close to the speed of light. But kinetic simulations have revealed that the distribution function of accelerated particles have fat tails and power-laws up to very large relativistic Lorentz factors (Cerutti et al. 2013; Melzani et al. 2014b; Sironi and Spitkovsky 2014; Werner et al. 2018; Ball et al. 2018). Different acceleration mechanisms are at work which will be discussed in Sect. 4.2.

**Driven reconnection** Reconnection sites are normally embedded in a large scale environment which is dynamic as well: jets, accretion disks, stellar winds, stellar atmospheres and coronae. Some of these environments are turbulent flows. As seen above,
this can decisively accelerate the REC process. But also directed large scale flows—as compared to the diffusion region or X-point where REC actually happens—can significantly accelerate REC in that they provide significantly higher inflow velocities. Therefore, much more magnetic flux can be carried from larger scales to the reconnection site. The timescale of the forcing also proves to be important (Pei et al. 2001; Pritchett 2005; Ohtani and Horiuchi 2009; Klimas et al. 2010; Usami et al. 2014, 2018).

**Multi-scale and multi-physics problem** As was seen so far, REC is a multi-scale problem. Large scale MHD flows can have a significant impact on the rate and the energetics of REC. Another scale is the transition to a diffusive regime which 'prepares' for REC, e.g., a Sweet-Parker reconnection sheet. Such sheets may break apart, introducing even smaller scales. This cascade in scales likely ends on kinetic scales. There also the physics may change, from a collisional to a collisionless regime. Another important point is the scale-difference in mass between electrons and ions, which also translates into differently scaled diffusive regions, the ion-diffusion region being about 42.85 (≡√*m*<sub>p</sub>/*m*<sub>e</sub>) times bigger than the electron diffusion region. And the different spatial lengths translate into equally different temporal scales. Magnetization and with it the ratio between an inertial length and the gyroradius yet complicates the situation.

But one has to address also other physical processes which influences the REC process. Outstanding here are radiative processes like synchrotron emission which directly changes the gyroradius. In an environment which is rich of photons, Compton scattering and Bremsstrahlung become important. More and more such processes are being addressed (Kirk and Skjæraasen 2003; Jaroschek and Hoshino 2009; Cerutti et al. 2013; Beloborodov 2017; Uzdensky 2016; Werner et al. 2019).

Multi-scale, multi-physics simulations demand for special techniques which are now in the course of being developed (Tóth et al. 2005; Daldorff et al. 2014; Tóth et al. 2012; Innocenti et al. 2013; Markidis et al. 2014; Ashour-Abdalla et al. 2015; Rieke et al. 2015; Lapenta et al. 2016; Tóth et al. 2016; Makwana et al. 2017; Lapenta et al. 2017; Lautenbach and Grauer 2018; Gonzalez-Herrero et al. 2018; Usami et al. 2018). We will come back to the issue in Sect. 4.2.

## 2.6 Laser plasma experiments

Over the past four decades, tremendous progress in the development of high-energy and high-power laser systems has brought the scientific community with the possibility to reproduce, in the laboratory, various scenarios relevant to astrophysics, space physics and planetology. This opened a new avenue for the development of so-called *Laboratory Astrophysics*, a field of growing activity that federates several communities [among which but not restricted to astrophysicists and (laser-)plasma physicists] and relies on the joint development of novel experimental and numerical capabilities.

In this section, we briefly review some key experiments focusing on the study of collisionless shocks and magnetic reconnection in laser-created plasmas<fnref n="7" /> The reader,

<footnote n="7">⁷ Other experimental facilities such as the Z-pinch machines (Remington et al. 2006) or the Large Plasma Device (LAPD) at UCLA (CA, USA) (Gekelman et al. 1991, 2016) also offer interesting opportunities for laboratory astrophysics, but are not discussed in this review.</footnote>
interested in other branch of laboratory astrophysics using laser-plasma experiments, will find interesting material covered in the review articles by Ripin et al. (1990), Rose (1994), Takabe et al. (1999), Remington et al. (1999), Remington et al. (2006), Takabe et al. (2008). These reviews cover topics ranging from warm dense matter, to equation of states and their application to planetology, opacities relevant to stellar interiors, or experiments investigating the hydrodynamics and magnetohydrodynamics of supernovae and (collisional) shocks.

In addition to presenting some of the main experimental results on collisionless shocks and magnetic reconnection, this section also aims at providing the reader with the characteristic parameters and conditions that can be created in the laboratory. To do so, we first introduce, in Sect. 2.6.1, the two main classes of lasers used for laboratory astrophysics. Then, in Sect. 2.6.2 we discuss the conditions that have to be met to ensure that collisional effects can be neglected. Finally, Sects. 2.6.3 and 2.6.4 summarize some of the key experimental results obtained on collisionless shock formation and magnetic reconnection.

## 2.6.1 Overview of laser facilities and characteristic parameters

Two main classes of laser systems are today used to support laboratory astrophysics research. First, high-energy density laser facilities delivering long (nanosecond) energetic (from few kJ up to 10s of kJ) light pulses have already allowed to reproduce various astrophysics-relevant scenarios, from warm dense matter studies, to the physics of hydrodynamic (radiative or not) shocks (Remington et al. 2006). Second, ultra-high intensity laser facilities deliver short (from few tens of femtoseconds to few picoseconds) light pulses that once focused onto a target allow to reach very high intensities. Even though laboratory astrophysics studies on this second class of laser systems is still in its infancy, recent developments of petawatt (and multi-petawatt) laser systems worldwide open new possibilities.

In this section, we report on some of the prominent laser facilities that are currently operating or will soon operate. Figure 11 lists these facilities as a function of the delivered energy and peak-power (the corresponding pulse durations are also indicated).

*High-energy density lasers* The development of high-energy density (HED) laser systems delivering energies of few tens of kilo-Joule (kJ) up to the Mega-Joule (MJ) over few to tens of nanoseconds (ns) has been strongly pushed forward by inertial confinement fusion programs (Atzeni and Meyer-Ter-Vehn 2004). Most of the experimental work that will be discussed in what follows has been performed on such laser systems. Various HED laser systems are today available, most of which are multi-beam facilities. Each beam can deliver ns pulses with few to 10 kJ (i.e., hundreds of beams are used on MJ-class laser systems) that, once focused onto target, allow to reach moderately high intensities of $10^{13}$ to a few $10^{15}$ W/cm$^2$. HED laser technology is based mainly on Nd:Glass amplifiers, which provide light beams at a (central) wavelength of ~ 1.05 µm, but often use frequency doubling of tripling techniques, so that the operating wavelength can be decreased to ~ 0.53 µm (doubling) or ~ 0.35 µm (tripling).
<figure><img src="image_19.png" /><figcaption>**Fig. 11** Prominent laser facilities presented as a function of the delivered laser energy and peak power. High-energy density (HED) laser facilities are shown in blue, ultra-high-intensity (UHI) laser facilities in green. The laser characteristics (energy, power and typical pulse duration) are indicative</figcaption></figure>

Among the prominent facilities are—at the multi-kJ level—the LULI 2000 laser<fnref n="8" /> in France, the Orion<fnref n="9" /> and VULCAN<fnref n="10" /> facilities in the UK, the GEKKO XII facility<fnref n="11" /> in Japan, and the Omega laser<fnref n="12" /> in Rochester, US. Two mega-joule-class lasers are also operating or under construction: the National Ignition Facility (NIF)<fnref n="13" /> in Livermore, California, started operating in the early 2010s. The Laser-MegaJoule (LMJ)<fnref n="14" /> is still under construction in the South-West of France. Note that the MJ-energy level is achieved by combining hundreds of 10 kJ nanosecond laser beams.

*Ultra-high intensity lasers* High-power ultra-high intensity (UHI) laser facilities provide light pulses of moderate energy (from few tens of Joule to few kilo-Joule) but of a very short duration (from tens of femtoseconds to a few picoseconds) that, when focused onto a target, allow to reach tremendous intensities (beyond 10¹<fnref n="8" /> W/cm²). At such intensities, electrons rapidly—in less than an optical cycle—become relativistic, and such UHI laser can help drive extremely fast, potentially relativistic, flows of plasmas.

Among the UHI lasers that are today considered for laboratory astrophysics studies, many are coupled to HED facilities. This is the case for instance of the PETAL<fnref n="15" /> and

<footnote n="8">⁸ https://portail.polytechnique.edu/luli/fr/installations/luli2000 (in French).</footnote>

<footnote n="9">⁹ https://www.awe.co.uk/what-we-do/science-engineering-technology/orion-laser-facility/.</footnote>

<footnote n="10">¹⁰ https://www.clf.stfc.ac.uk/Pages/Vulcan-laser.aspx.</footnote>

<footnote n="11">¹¹ http://www.ile.osaka-u.ac.jp/eng/facilities/gxii/index.html.</footnote>

<footnote n="12">¹² http://www.lle.rochester.edu/omega_facility/omega/.</footnote>

<footnote n="13">¹³ https://lasers.llnl.gov/.</footnote>

<footnote n="14">¹⁴ http://www-lmj.cea.fr/.</footnote>

<footnote n="15">¹⁵ http://www.enseignementsup-recherche.gouv.fr/cid99515/petawatt-aquitaine-laser-petal.html (in French).</footnote>
NIF-ARC<fnref n="16" /> petawatt-class lasers coupled to the LMJ and NIF facilities, respectively, which deliver petawatt-level light pulses with energy of a few kJ and duration in the picosecond range. LULI 2000, VULCAN and ORION also have short (picosecond) pulse beamlines that deliver energy up to a few 100s J.

Other UHI facilities are also available that are not coupled to HED laser systems. This is the case e.g. of the femtosecond laser GEMINI<fnref n="17" /> in UK. Delivering 15 J in 30 fs, GEMINI has been used e.g. to produce dense electron-positron clouds which were used to drive current instabilities in a Helium plasma (Warwick et al. 2017). Let us further note that the most powerful laser is today the CoReLS<fnref n="18" /> in Gwanju, South-Korea, that has recently delivered 4 PW pulse (Nam et al. 2018). In addition the Apollon<fnref n="19" /> laser (in construction on the Plateau de Saclay, 20 km south of Paris, France) and the ELI project<fnref n="20" /> aim at reaching the unprecedented power of 10 PW within the next few years (in light pulses of a few tens to few 100s fs). Laboratory astrophysics studies are envisioned on these facilities.

## 2.6.2 The collisionless regime

As previously stated, this section focuses on collisionless laser-plasma experiments and on the physics of collisionless shocks and magnetic reconnection in particular. The first observations of a collisionless coupling in laser-created plasmas date back to the early 1970s (Cheung et al. 1973), and very early the question collisionality effects arose, see, e.g., the work by Dean et al. (1971) and following exchange (Wright 1972; Dean et al. 1972).

The first (theoretical) investigations of collisionless shock experiments actually addressed this issue (Drake and Gregori 2012; Park et al. 2012; Ryutov et al. 2012). These works proposed the first designs and scaling laws to reproduce electrostatic or Weibel-mediated shocks (see Sect. 4 for complementary definitions) in laser-plasma experiments, and addressed the potential effects of particle collisions (and how to mitigate them) in counter-streaming plasma flows.

Of particular importance are collisions in between counter-streaming ions of the two flows (inter-flow collisions) that can have a dramatic effect on the shock formation. Indeed, and as stressed by Drake and Gregori (2012), the mean-free-path for ion-ion collisions measures the length over which an ion (subject to multiple scatterings/collisions) sees its velocity deflected by 90°. Hence, collisional effects will effectively isotropize the flow over a characteristic given by this mean-free-path and collisional shocks are known to develop on this spatial scale. Conducting a collisionless shock experiment thus requires that this (inter-flow) collision *greatly* exceeds the characteristic length of shock formation. As will be detailed in the following Sect. 2.6.3, it turns out that this condition can be “easily” met in electrostatic and to some extent magnetized shock experiments. In the case of Weibel-mediated shocks, entering the

<footnote n="16">¹⁶ https://lasers.llnl.gov/science/photon-science/arc.</footnote>

<footnote n="17">¹⁷ https://www.clf.stfc.ac.uk/Pages/Laser-system-Gemini.aspx.</footnote>

<footnote n="18">¹⁸ https://www.ibs.re.kr/eng/sub02_03_05.do.</footnote>

<footnote n="19">¹⁹ http://cilexsaclay.fr/Phocea/Vie_des_labos/Ast/ast_technique.php?id ast=9 (in French).</footnote>

<footnote n="20">²⁰ https://eli-laser.eu/.</footnote>
collisionless regime requires extremely fast (several 1000s km/s) flows overlapping for
a sufficient time accessible only on the most energetic (MJ-class) laser systems (Park
et al. 2012).

In addition to inter-particle collisions, internal collisions between ions of the same
flow can also be of importance, in particular as the ion temperature in the flow is
quite low (at least before shock formation). This issue is briefly addressed by Drake
and Gregori (2012) and Ryutov et al. (2012). Yet their impact on the development of
instabilities such as the ion Weibel instability, or shock formation remains unclear.
While it is difficult to claim that these (internal) collisions may not strongly modify
the physics of shock formation, this may be checked by careful numerical modeling,
e.g. relying on kinetic (Particle-In-Cell) simulations including collisional effects (see
Sect. 3.3.2).

Last, Drake and Gregori (2012) discussed the possible impact of electron-ion col-
lision on the dissipation of the magnetic structures that play a central role in the
formation of Weibel-mediated shocks; and on longer time on particle acceleration.
The authors showed that such collisions may indeed impact the small scale magnetic
structures, but will most likely no impact the larger scale structures that develop on
the scale of the ion skin-depth and thus the formation of Weibel-mediated shocks.

2.6.3 Collisionless shock experiments

As previously mentioned, evidences of collisionless processes in the presence of
counter-streaming laser-produced plasmas were reported is the early 1970s. The first
reported observation of a collisionless shock in a laser-created plasma<sup>21</sup> dates back
to Bell et al. (1988). This experiment was conducted on the VULCAN laser (Ross
et al. 1981) at the Rutherford Laboratory (UK) where two laser pulses, each deliv-
ering 120J over 18 ns (FWHM), were focused in a 50 µm-diameter spot onto a flat
carbon target. The resulting laser-produced ablation plasma had a density ~ 10<sup>18</sup> cm<sup>-3</sup>
and velocity of a few 100s of km/s. It collided with an obstacle (located 250 µm away
from the ablated target). The experiment led to the formation of density structures that
were interpreted as collisionless bow shocks. In this experiment, all mean-free-paths
were larger than the mm, while the width of the observed shock front ranged from 0.01
to 0.05 mm. The nature of the shock—either electrostatic or weakly magnetized—was
however not fully defined.

*Electrostatic shocks* Following this pioneering work, and since the late 2000s in particular, collisionless electrostatic shocks have been abundantly produced in laser-plasma experiments. These later developments were accompanied by both strong developments in diagnostics, and the use of kinetic (Particle-In-Cell) simulations to support the experimental effort.

For instance, Romagnani et al. (2008) demonstrated the creation of an electrostatic
shock following the sudden expansion of a plasma into a rarefied gas. In this experiment
carried out on the LULI 100 TW laser facility, one laser pulse with duration 470 ps
and energy of a few tens of J was focused onto a Tungsten or Aluminium foil. Quickly

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

<footnote n="21">21 Collisionless shock waves were obtained in plasma experiments, albeit not using lasers, since the mid-1960s [see Strokin (1985) and references therein]. Already these studies where motivated by space-plasmas and astrophysics.</footnote>
<figure><img src="image_21.png" /><figcaption>**Fig. 12** Left: typical proton radiographies of laser-driven electrostatic shocks. Region I shows strong modulations associated with the ablating plasma, regions II and III show different structures that are interpreted as shock waves propagating ahead of the ablating plasma, while the modulated pattern in Region IV is located ahead of the shock front and possibly associated with a reflected ion bunch. The arrow indicates the laser beam direction. Right: shock structure 150 (b), 160 (c) and 170 ps (d) after the beginning of the interaction. The different times are accessible, for a single shot, by selecting protons with different energies [11.5 MeV for panel (b), 10 MeV for panel (c) and 9 MeV for panel (d)] as their time of flight from their source to the shock structure is different. Images reproduced with permission from (left) Romagnani et al. (2008) and (right) from Ahmed et al. (2013), copyright by APS</figcaption></figure>

heated, the ablated foil expanded in the surrounding media and drove the formation of a collisionless electrostatic shock wave, about 1 mm away from the target, that was propagating at a velocity close to the ion acoustic velocity ~200–400 km/s. The shock was diagnosed using proton radiography (Borghesi et al. 2001), a technique that is now central to the study of collisionless shock in laser-plasma experiments. It relies on the deflection (in the electromagnetic fields developed at the shock front) of protons created by a second ultra-short (~300 fs) ultra-intense (≥ 10¹⁸ W/cm²) laser pulse. The proton radiography is recorded onto dosimetrically calibrated radiochromic films (RCFs), as shown in Fig. 12.

Other experiments (Kuramitsu et al. 2011; Ahmed et al. 2013; Morita et al. 2013) have similarly reported the formation of electrostatic collisionless shock waves using ablating plasmas, either in direct interaction with a standing (background) plasma, or in counter-streaming plasma flows.

**Magnetized shocks** The first observation of a magnetized collisionless shock motivated by astrophysics studies was claimed by Niemann et al. (2014) combining the use of a 25 ns - 200 J laser and the LAPD. In this experiment, the LAPD was used to produce a large scale (17 m × 0.6 m) low density (10¹²–10¹³ cm⁻³) and temperature (Tᵢ = 1 eV, Tₑ = 6 eV) hydrogen plasma embedded in an external magnetic field B₀ = 300 G. The 10¹³ W/cm² ns laser pulse was fired at a solid polyethylene target embedded inside the magnetized plasma, which launched a denser (8 × 10¹⁶ cm⁻³) slightly warmer (Tₑ ~ 7.5 eV) carbon ion plasma at a velocity of ~500 km/s directed perpendicular to the magnetic field. The interaction of this super-Alfvénic plasma with the ambient (LAPD) plasma led—through a collisionless coupling—to the formation of a magnetic piston and then to the formation of self-sustained magnetosonic shock, supported by the ambient ions and propagating away from the piston at a velocity of ~370 km/s (corresponding to an Alfvénic Mach number Mₐ ~ 2). The reported
<figure><img src="image_22.png" /><figcaption>**Fig. 13** Magnetized shock experiment at the Large Plasma Device (LAPD, University of California Los Angeles). a Magnetic (B) field structure as a function of position (x is the direction of the carbon plasma flow, transverse to the direction of the background magnetic field) and time. b Temporal evolution of the magnetic field at x = 35 cm with (black) and without (red) the ambient plasma. c Spatial profile of the magnetic at two different times, 0.3 µs (before shock formation, dashed light blue line) and 0.7 µs (after shock formation, solid black line). Image reproduced with permission from Niemann et al. (2014), copyright by AGU</figcaption></figure>

measurements (shock velocity and magnetic field compression $B/B_0 \sim 2$) were found to be consistent with Rankine–Hugoniot conditions as well as with two-dimensional, collisionless, simulations performed using an electromagnetic Darwin code (Winske and Gary 2007). Note that, as illustrated in Fig. 13, for this particular experiment, the use of the LAPD allowed to follow the magnetic piston and shock formation over large spatial (few tens of cm) and temporal (few microseconds) scales, well beyond what is usually accessible using HED or UHI laser systems.

The first laboratory observation of a laser-driven high-Mach-number magnetized collisionless shock was reported by Schaeffer et al. (2017a). This experiment was conducted on the Omega EP laser facility at Rochester (US) and built up on the concept of magnetic piston used in, e.g, the previous experiment by Niemann et al. (2014). However, it relied solely on the use of the HED laser Omega EP and allowed to create super-critical magnetized shocks with (magnetosonic) Mach number $M_{ms} \equiv u_{sh}/c_{ms} \sim 12$, with $u_{sh} \sim 700$ km/s the measured shock velocity ($c_{ms}^2 = u_A^2 + c_s^2$, $u_A$ and $c_s$ being the (upstream) Alfvén and ion acoustic velocities, respectively). To

<figure><img src="image_23.png" /></figure>
<figure><img src="image_24.png" /><figcaption>**Fig. 14** Magnetized shock experiment at the Omega EP laser facility. (Left) simulation setup. (Right) panels (a–e) present the angular filter refractometry measurements for different configurations: a no background plasma, no external magnetic field (no shock is observed); b presence of a background plasma but no external magnetic field; c with both background plasma and external magnetic field but only one piston plume; d, e in the presence of the two piston plumes and background plasma, with external magnetic fields in parallel and anti-parallel configuration, respectively. f Proton radiography signal revealing the strong magnetic field compression at the location of the shock structure. Images adapted from Schaeffer et al. (2017a, b)</figcaption></figure>

do so, various beams of the Omega EP laser were used. A first beam, with energy 100 J and duration 1 ns, was focused (at intensity ~ $10^{12}$ W/cm$^2$) onto a CH target thus producing a background plasma. A second and third beam, with energy 1.5 kJ and duration 2 ns, were then focused onto two opposing CH targets, leading to the production of two counter-propagating ablation plasmas. Even though the details of the resulting three plasma flows are not fully documented in either (Schaeffer et al. 2017a) or the companion paper (Schaeffer et al. 2017b), the density and temperature of the overlapping plasma were estimated to ~ $6 \times 10^{18}$ cm$^{-3}$ and $T_e$ ~ 15 eV, respectively. The whole set-up was embedded into a 80 kG perpendicular magnetic field produced by a pulsed current passing through Copper wires located behind the 2 opposing CH targets (for the readers convenience, the experimental set-up is reproduced in the left panel of Fig. 14). The resulting magnetic piston and shock structures are evidenced in the right panels of Fig. 14. The importance of the background plasma (created by the first 100 J-beam) is made clear by comparing panel (a) to panels (b–e). Without background plasma, panel (a), no shock is observed. In the presence of a background plasma, panels (b–e), a shock-like structure is observed in all cases, strongest in the presence of the external magnetic field [panels (c–e)], but still present when no external magnetic field is applied [panel (b)]. The authors advance the possibility, supported by PIC simulations, that in this latter case, the Biermann-Battery process was responsible for seeding a large scale magnetic field even though no external one is applied. Note also, that in panel (c), the authors report the creation of a shock when only one of the 1.5 kJ-beam was used, demonstrating that only one piston plume interacting with the background plasma was needed to produce a magnetized shock. Last, panel (f) reports the measurement obtained using proton radiography and reveals a strong magnetic field compression at the location of the shock structure.
Note that the overall experimental campaign strongly relied on advanced diagnostics [shadowgraphy, angular filter refractometry (Haberberger et al. 2014), and proton radiography] as well as the combined use of hydrodynamic (for the plasma characterization) and PIC (for the shock formation and evolution) simulations.

Since these experiments, the effort in producing and studying magnetized shocks has continued, e.g., exploring the possibility to produce parallel shocks (Weidl et al. 2017), or to use the Biermann-Battery process to magnetize the plasmas (Umeda et al. 2019).

Weibel-mediated shocks Weibel-mediated collisionless shocks are certainly the most sought after collisionless shocks in laser-based laboratory astrophysics experiments. Even though important progress have been made in the last decade, recreating such shocks in the laboratory has not yet been achieved. The main difficulty in recreating such shocks stems from the need to achieve flow densities that are, on the one hand, sufficiently small to ensure that one operates in the collisionless regime and, on the other hand, large enough for the ion Weibel instability to develop and the resulting magnetic turbulence to build up. The combined experimental and theoretical effort in this endeavor has been started since the early 2010s and HED NIF-class laser systems, allowing to produce plasma flows with densities of a few $10^{19}$ cm$^{-3}$ and velocities of several 1000s km/s, have been identified as the most promising path toward collisionless Weibel-mediated shock formation.

The first step toward creating plasma flows relevant for such studies was taken by Park et al. (2012), demonstrating the possibility to drive plasma flows with velocities of ~ 1000 km/s and densities of ~ $10^{18}$ cm$^{-3}$ from plasma ablation at the Omega laser facility. The production of large-scale electromagnetic structures (Kugland et al. 2012) and later clear demonstration of Weibel-type ion filamentation instabilities (Fox et al. 2013; Huntington et al. 2015) in the presence of counter-streaming plasmas were obtained by irradiating a pair of opposing plastic (CH) foils with few kJ, few ns laser pulses on the Omega EP laser system. Figure 15 reproduces the schematic experimental set-up used by Huntington et al. (2015) together with a typical proton radiography measurement of the magnetic field filamentary structures following from the development of the ion Weibel instability. The region imaged by the proton radiograph is about 3 mm wide, and the filamentary structures have a typical width of ~ 150–300 µm, consistent with the ion skin-depth for the reported plasma density of a few $10^{18}$ cm$^{-3}$.

The plasmas flows achievable at Omega are unfortunately too low density, and short life, to allow creating a Weibel-mediated collisionless shock. A recent theoretical model and 2D PIC simulations by Ruyer et al. (2016) indeed predict that the isotropization necessary for shock formation may be achievable if the two counter-streaming flows overlap over a length of the order of at least:

$$L_{\text{iso}} \approx 35 [m_i/(Zm_e)]^{0.4} (c/\omega_{pi}) \rightarrow 5 \text{ cm} \times \left(\frac{A}{Z}\right)^{0.9} \sqrt{\frac{10^{19} \text{ cm}^{-3}}{n_0}}, \quad (21)$$

where $c/\omega_{pi}$ is the ion skin-depth associated to the plasma flow density $n_0$, $m_i$ ($m_e$) the ion (electron) mass, and $Z$ ($A$) the ion charge (mass) number. Conversely, the two flows shall overlap for a time of the order of $\tau_{\text{iso}} = L_{\text{iso}}/v_0 \sim$, with $v_0$ the relative flow

<figure><img src="image_25.png" /></figure>
<figure><img src="image_26.png" /><figcaption>Fig. 15 Demonstration of the ion Weibel instability in counter-streaming plasmas on the Omega EP facility. Simulation setup: two ablation plasmas are formed in counter-streaming configuration by irradiating two opposing plastic targets. A $^3$He-D target is imploded to drive fusion reactions that allow for the production of 3 MeV and 14.7 MeV protons that are used to radiograph the electromagnetic fields developed in the region where the two counter-streaming plasma overlap. The resulting proton radiograph (here taken about 5 ns after the beginning of interaction using the 14.7 MeV protons) shows evidence of the growth of filamentary structures associated to the Weibel magnetic fields. Image reproduced by permission from Huntington et al. (2015), copyright by Macmillan</figcaption></figure>

velocity. A necessary, yet not sufficient, condition for maintaining the collisionless regime imposes that the isotropization length given by Eq. (21) is much smaller than the characteristic ion-ion collision mean-free-path (Park et al. 2012):

$$ \lambda_{mfp} \approx 5 \text{ cm} \times \frac{A^2}{Z^4} \left( \frac{v_0}{1000 \text{ km/s}} \right)^4 \left( \frac{10^{19} \text{ cm}^{-3}}{n_0} \right), \quad (22) $$

for which are considered only collisions between ions of the two counter-streaming flows. These estimates allow to infer that Weibel-mediated collisionless shocks may be achieved in the presence of (hydrogen) plasma flows with density of the order of $10^{19} \text{ cm}^{-3}$, colliding at a relative velocity of at least a few 1000s km/s, provided these flows overlap over distance of a few cm during a few to tens of ns. Such conditions can be met only at the most energetic, MJ-class laser systems such as NIF.

The first experiments to recreate Weibel-mediated shocks at NIF have been started in the framework of the *Discovery Science* program. The first experimental results were reported by Ross et al. (2017) focusing on how to tune the experimental conditions to access the collisionless regime. These experiments, for which no external magnetic field was used, considered two solid targets, made of a mixture of Carbon and either hydrogen or deuterium (CH or CD), each irradiated by several beams allowing to deliver energies of ~ 250 kJ (much larger than accessible e.g. at Omega) per foil during about 5 ns. The resulting ablation plasma flows at velocities of ~ 1000 km/s and (ion) densities of a few $10^{19} \text{ cm}^{-3}$ in the interaction (overlapping) region. This work demonstrated that if the foils (in opposing configuration) were sufficiently distant
from one another, collisional effects could be strongly mitigated due to the reduction of
the plasma flow density in the overlapping region. Most importantly, this work reported
evidence of a collisionless (collective) heating in the flow interaction region, which
the authors associate with the nonlinear stage of the Weibel instability and thus to the
early stage of shock formation. This result suggests that the scientific community is on
the verge of producing Weibel-mediated collisionless shocks in the laboratory. A new
experiment was actually conducted at NIF in the last months increasing the driving
laser beam energy to ~ 500 kJ delivered to each foil. This experiment is expected to
lead to shock formation, and may also gives the first signs of particle acceleration in
the shock. To this date, the results of this last experimental campaign have not been
announced.

*Prospective numerical studies*
The possibility to drive collisionless shocks in the
laboratory has prompted the laser-plasma community to investigate various laboratory
configurations to drive collisionless collective processes and shocks *in silico*. Indeed,
various numerical experiments have been performed using Particle-In-Cell (PIC) sim-
ulation. Even if some of these numerical experiments consider laser parameters not yet
within our reach (ultra-short ps-level, energy at the 100s J level) other have addressed
conditions that are or will soon be achievable on the forthcoming extreme light facil-
ities such as Apollon or ELI.

Fiuza et al. (2012) put forward the possibility to drive—through Weibel-like instabilities—a collisionless shock in a dense target using an ultra-intense light pulse. This scenario was revisited by Ruyer et al. (2015) that demonstrated the dominant role of laser-driven hot electrons in the shock formation. More recently, Grassi et al. (2017) showed that tuning the laser-plasma interaction configuration can help mitigate the hot electron production so that shock formation can be driven by the ion Weibel instability, as expected in astrophysical scenarios.

In addition, dense electron-positron flows have been produced in laser-plasma
experiments (Chen et al. 2015; Sarri et al. 2015), offering the opportunity to study
pair-plasma processes in the laboratory, and motivated various numerical experiment.
Using QED-PIC simulation, Lobet et al. (2015) demonstrated the possibility to drive
ultra-relativistic, counter-propagating electron-positron pair plasmas using extreme
light pulses (with intensity beyond 10<sup>23</sup> W/cm<sup>2</sup>, 100s kJ and duration of few tens of
fs). Ultra-fast isotropization and thermalization (a first step in shock formation) were
observed in the simulation, and associated to both the Weibel instability and a remark-
able contribution from synchrotron emission by the ultra-relativistic leptons in the
strong (Weibel) magnetic fields. More recently, the (collisionless) interaction of midly
relativistic pair jets with background (electron-ion) plasma was also investigated in
kinetic simulations (Dieckmann et al. 2018a, b). Remarkably, these studies are not only
motivated by astrophysics (Dieckmann et al. 2019) but also by recent experiments that
demonstrated the growth of a current-driven instability developing during the inter-
action a quasi-neutral pair beam with a background (electron-ion) plasma (Warwick
et al. 2017).
<figure><img src="image_27.png" /><figcaption>**Fig. 16** Adapted from the work by Nilson et al. (2006) presenting the first magnetic reconnection experiment with laser-created plasmas. **a** Experimental set-up using the two-spot configuration. **b** Proton-radiography measurements showing in dark the regions where protons are recorded. **c** Shadowgraphy measurements indicating the formation of two jets at the reconnection layer</figcaption></figure>

### 2.6.4 Magnetic reconnection experiments

Laser-plasma experiments also provide a test bed for magnetic reconnection studies.<fnref n="22" /> The first evidences for magnetic reconnection in a laser-plasma experiment were reported by Nilson et al. (2006). This experiment was performed at the VULCAN laser facility in the UK, and relied on a now popular set-up consisting in firing two laser pulses at a solid (here Aluminium or Gold) target (see Fig. 16a). The interaction of each pulse leads to plasma ablation and expansion associated with the generation of an azimuthal magnetic field through the Biermann-Battery process. In between the two pulses, the magnetic fields driven by the two pulses are in an anti-parallel configuration, and a reconnection layer can form.

This experiment was carried out in the HED regime of interaction, each laser pulse of the VULCAN facility delivering 200J over 1 ns in a 30–50 µm focal spot (the corresponding laser intensity is moderate ~ 10¹⁵ W/cm²). Various complementary diagnostics were used, highlighting features consistent with magnetic reconnection. (i) Proton radiography (Borghesi et al. 2001) allowed to probe the generated magnetic fields. A typical measurement is reproduced in Fig. 16b; and an additional analysis of these measurements is given in Willingale et al. (2010). Light regions correspond to regions free of protons, i.e., to regions where the strong Biermann-Battery magnetic field (estimated to be of the MG-level in this experiment few 100s of ps after the ablated plasma started expanding) is present.

The presence of a strong proton signal (dark region) in between the two lighter blobs was identified as the reconnection layer, where opposite magnetic field lines can reconnect and lead to a null-magnetic field region. (ii) In addition, the interaction region was also probed by a short (10ps) light pulse allowing to produce a shadowgraphy [as well as an interferogram (not shown here)] of the interaction region. Such a shadowgram is reproduced in Fig. 16c and shows the formation - on the ns-timescale-

<footnote n="22">²² In this review, we focus once more on laser-plasma experiments. Yet, other experiments have been conducted on various devices: the Magnetic Reconnection Experiment (MRX) at Princeton (Yamada et al. 1997), the LAPD at UCLA (Gekelman et al. 2010), Z-pinch machines (Hare et al. 2017). Magnetic reconnection is also known to affect Tokamak experiments (Goetz et al. 1991). See also Howes (2018) for a review of various laboratory experiments for space plasma physics.</footnote>
of two distinct jets [see original paper and Nilson et al. (2008) for more details] with velocities of ~500 km/s, which would not be expected should only hydrodynamic processes govern the plasma evolution. (iii) Finally, Thomson Scattering measurements (not shown here) showed that, while the electron temperature was of the order of 800 eV (at 1.5 ns) and 700 eV (at 2.5 ns) in the ablated plasmas, a much higher electron temperature ~1.7 keV (at 1.2 ns) was measured in what was identified as the reconnection layer (where there is no laser). Such high electron temperatures were also put forward as a result of magnetic reconnection; and were consistent with supporting hybrid simulations.

Since this first experiment, similar results were obtained on other HED laser facilities. Li et al. (2007) report on an experiment performed at the OMEGA laser facility using a slightly more energetic 1ns laser pulse (500J), with a spot diameter of ~800 µm (corresponding to an intensity of ~ $10^{14}$ W/cm$^2$). This experiment benefited from a high-quality proton radiography (monoenergetic 14.7 MeV protons were produced by fusion reactions from an imploded $D^3$-He target), which allows the authors to probe the changes in the magnetic field topology as magnetic reconnection proceeds. See also Rosenberg et al. (2015b,a). Similarly, Zhong et al. (2010) reported on a similar experiment carried out on the Shenguang II (SG II) laser facility in Shanghai, China. In this experiment, four laser beams (1ns, few 100s J, 50–100 µm-wide spots) are used to drive the plasma expansion in the two-spot configuration previously discussed, but shining the lasers on the front and back side of the target simultaneously. This experiment also put the accent on scaling their results with respect to reconnection outflows in solar flares. As an example, relying on X-ray imaging of the interaction region, the authors could demonstrate a change in the directionality of the jets due to an asymmetry in the driving laser intensities,<fnref n="23" /> as shown in Fig. 17.

Another experimental set-up was also proposed by Fiksel et al. (2014) and conducted on the OMEGA EP laser. In contrast with the previous experiments, this new set-up relies on (i) an head-on configuration with two targets (irradiated by kJ, ns laser beams), (ii) current-carrying conductors placed behind the two targets to create an external (up to 80 kg) magnetic field imposed perpendicular to the expanding plasma flows and designed such that the field a x-type null point in between the two targets; and (iii) the presence of a background plasma created by a third (100J, 1ns) laser beam. Proton radiography measurements indicate the formation and collision of magnetic ribbons, pileup of magnetic flux and reconnection which are found to be in remarkably good agreement with 2D PIC simulations (that include particle collisions).

While the previous experiments were conducted in a collisional regime, more recent experiments have focused on the collisionless regime. Dong et al. (2012) conducted an experiment on the SG II laser, using a similar set-up then previously presented by Zhong et al. (2010), but firing the laser beams (450J on each target) at two Al-targets separated by 150–240 µm, which, even though the collisionless nature of the reconnection region is not fully addressed, lead the authors to claim the study of a

<footnote n="23">²³ See also the work by Rosenberg et al. (2015a) for a study of asymmetric reconnection OMEGA laser facility.</footnote>
**Fig. 17** X-ray imaging showing two bright spots in the Al target were the expanding plasmas are heated by the laser beams. Here the asymmetry in the laser intensities (in regions B1 and B2) leads to an inclination of the upward flow. Below the Al-target, a Cu-target is placed. The downward flow impinges on this target and results in a hot X-ray source. Image reproduced with permission from Zhong et al. (2010), copyright by Macmillan

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

structure of collisionless reconnection. These authors also report on the ejection of a plasmoid that, when it rapidly propagates away, deforms the reconnected magnetic field and generate a secondary current sheet. This process seems to be well reproduced by PIC simulations, and the primary reconnection event is found to be associated with well-collimated plasma outflows containing high-energy (MeV) electrons.

In addition, Raymond et al. (2018) conducted experiments, on both the OMEGA EP facility and the HERCULES laser at University of Michigan, in a regime where magnetic reconnection was not only collisionless but also driven by relativistic electrons. This was made possible by using short pulse laser beams (20 ps for the OMEGA EP laser, 40 fs for the HERCULES laser) in a configuration otherwise similar to that (the two-spot experiment) initially proposed by Nilson et al. (2006). Using short pulses indeed allowed to reach ultra-high intensities (~ $10^{18}$ W/cm$^2$ on OMEGA EP and ~ $2 \times 10^{19}$ W/cm$^2$ on HERCULES), thus allowing to enter the relativistic regime of laser-plasma interaction. Figure 18 shows the typical experimental set-up as well as a typical X-ray imaging where the two heated and expanding plasmas can be seen, together with a reconnection layer in between. In addition, the authors report the formation of a non-thermal (few MeV) electron population whenever reconnection is expected, consistent with supporting 3D PIC simulations.
<figure><img src="image_29.png" /><figcaption>**Fig. 18** Experimental set-up relying on ultra-high-intensity short pulse laser beams allowing to probe relativistic reconnection. A typical X-ray imaging shows the location of the two hot expanding plasmas and, in between, the reconnection layer. Image reproduced with permission from Raymond et al. (2018), copyright by APS</figcaption></figure>

# 3 Solving kinetic problems

This section reviews the model equations used to describe particle kinetic physics, i.e., the dynamics of charged particles in the configuration and momentum space under the effect of electro-magnetic forces. The section is divided as follows: in Sect. 3.1 we describe the Vlasov-Maxwell system of equations, then in Sect. 3.2 we discuss the numerical methods developed to follow the dynamics of such a system. Section 3.3 discusses the particle-in-cell (PIC) technique used to study solutions of the Vlasov-Maxwell system. In Sect. 3.4 we provide a discussion on the comparison between PIC and Vlasov approaches. Section 3.5 briefly describes hybrid methods where a fluid approximation is introduced for the electronic component whereas kinetic (PIC) techniques are used to describe ions. In Sect. 3.6 we specifically discuss the Fokker-Planck description of kinetic problems. The Fokker-Planck approach is particularly well-adapted to investigate cosmic ray propagation. Finally, we give a particular focus on Fokker-Planck simulations developed in the context of the study of the radiative transfer in hot plasmas around compact objects.

## 3.1 The Vlasov-Maxwell description of a collisionless plasma

### 3.1.1 Governing equations

Let us consider a plasma composed of various species, labeled *s*, corresponding to particles with mass $m_s$ and charge $q_s$. The kinetic description of this plasma relies on the representation of each species *s* by its (one-particle) distribution function $f_s(t, \mathbf{x}, \mathbf{p})$, $f_s(t, \mathbf{x}, \mathbf{p})d^3\mathbf{x}d^3\mathbf{p}$ measuring, at any time $t$, the (probable) number of particles of species *s* in a volume element $d^3\mathbf{x}d^3\mathbf{p}$ at a position $(\mathbf{x}, \mathbf{p})$ in phase-space $(\mathbf{x} \text{ and } \mathbf{p})$

<figure><img src="image_30.png" /></figure>
standing for the spatial and momentum coordinates, respectively). In the absence of collisions, the evolution of the distribution functions $f_s$ satisfies the Vlasov equation:

$$ \frac{\partial}{\partial t} f_s + (\mathbf{v} \cdot \nabla) f_s + (\mathbf{F}_s \cdot \nabla_p) f_s = 0, \qquad (23) $$

where $\mathbf{v} = \mathbf{p}/(m_s\gamma)$ is the velocity corresponding to a particle of momentum $\mathbf{p}$ and Lorentz factor $\gamma = \sqrt{1 + \mathbf{p}^2/(m_s c)^2}$ ($c$ is the speed of light in vacuum), and $\mathbf{F}_s$ is the force acting on the species particles. As this work focuses on electromagnetic plasmas, this force is the Lorentz force exerted by the collective electric $\mathbf{E}$ and magnetic $\mathbf{B}$ fields:

$$ \mathbf{F}_s = q_s \left( \mathbf{E} + \frac{1}{c} \mathbf{v} \times \mathbf{B} \right). \qquad (24) $$

It is important to stress that, in the Vlasov equation, the electromagnetic fields are collective fields, also referred to as macroscopic fields in the sense that they do not account for the microscopic variations developing at the particle scale. Hence collisions are not considered in this description.

The electric $\mathbf{E}(t, \mathbf{x})$ and magnetic $\mathbf{B}(t, \mathbf{x})$ fields are in general functions of both space and time, and satisfy Maxwell's equations<fnref n="24" />:

$$ \nabla \cdot \mathbf{E} = 4\pi\rho, \qquad (25a) $$

$$ \nabla \cdot \mathbf{B} = 0, \qquad (25b) $$

$$ \nabla \times \mathbf{E} = -\frac{1}{c} \frac{\partial \mathbf{B}}{\partial t}, \qquad (25c) $$

$$ \nabla \times \mathbf{B} = \frac{4\pi}{c} \mathbf{J} + \frac{1}{c} \frac{\partial \mathbf{E}}{\partial t}. \qquad (25d) $$

The electromagnetic fields act onto the plasma through the Lorentz force (24), and are in turn modified by the plasma through the total charge and current densities, $\rho = \sum_s \rho_s$ and $\mathbf{J} = \sum_s \mathbf{J}_s$, respectively, where each species charge and current densities are defined as:

$$ \rho_s(t, \mathbf{x}) = q_s \int d^3 p \, f_s(t, \mathbf{x}, \mathbf{p}), \qquad (26a) $$

$$ \mathbf{J}_s(t, \mathbf{x}) = q_s \int d^3 p \, \mathbf{v} f_s(t, \mathbf{x}, \mathbf{p}). \qquad (26b) $$

The coupled system of Eqs. (23) and (25), together with the Lorentz force (24) and definitions of the charge and current densities (26) form the Vlasov-Maxwell model. It provides a self-consistent, kinetic description for the evolution of a collisionless plasma and the associated collective electromagnetic fields.

<footnote n="24">²⁴ In some cases where only electrostatic fields are important, only Poisson Eq. (25a) may be used. The Vlasov-Poisson description is of particular importance for the description of cold plasmas in particular, as considered e.g., for plasma propulsion, see Ref. (Boeuf and Garrigues 2018).</footnote>
### 3.1.2 Initial and boundary conditions

The Vlasov-Maxwell model relies on a system of partial differential equations and thus requires initial and boundary conditions. The initial condition of the system (at time $t=0$) consists first in defining the initial distribution functions $f_s(t=0, \mathbf{x}, \mathbf{p})$ for all species $s$ of the system. One usually considers equilibrium$^{25}$ distribution functions, and Maxwellian or Maxwell-Jüttner distribution functions (drifting or not) are often considered.$^{26}$ The initial electromagnetic fields spatial distribution $\mathbf{E}(t=0, \mathbf{x})$ and $\mathbf{B}(t=0, \mathbf{x})$ also needs to be prescribed. At $t=0$, these fields have to satisfy Eq. (25a) and Eq. (25b), respectively. Hence, $\mathbf{B}(t=0, \mathbf{x})$ has to be divergence free, while $\mathbf{E}(t=0, \mathbf{x})$ can be either divergence free (e.g., if an external electric field is considered) or has to be computed from Poisson's Eq. (25a) using the initial distribution functions $f_s(t=0, \mathbf{x}, \mathbf{p})$ to compute the initial charge density.

Various boundary conditions (BCs) can be considered and will strongly depend on the physics at hand. First, BCs on the distribution functions can be used to reflect, thermalize already existing particles or inject new particles at the border of the spatial domain. In addition, when directly solving the Vlasov equation in phase-space (see Sect. 3.4 on so-called Vlasov codes), BCs on the momentum component have to be considered. Similarly, electromagnetic fields can be reflected, absorbed or injected at the domain border by prescribing the correct BCs for the electric and magnetic fields.

## 3.2 Solving the Vlasov-Maxwell system numerically: general considerations

Computer simulation is an outstanding tool for solving the Vlasov-Maxwell system of equations together with the prescribed initial and boundary conditions, and most of today's kinetic simulations of plasmas rely on massively parallel tools to do so. In what follows, we present two of the main numerical approaches to solve this system. The first method is used in so-called Vlasov codes, while the second is used in so-called Particle-In-Cell (PIC) codes. The main difference between the two methods lies in the way they solve the Vlasov equation. Otherwise, both methods follow the same procedure which rely on discretizing the fields onto a spatial grid (henceforth referred to as the simulation grid), advancing the distribution function and then updating the associated charge and current densities onto the simulation grid. This procedure is here briefly detailed and summarized in Fig. 19.

### 3.2.1 Initialization and time-loop

First, the initialization step consists in prescribing the distribution functions for all species $s$ at time $t=0$ together with the initial electric and magnetic fields. Again, one should stress that both fields must satisfy Eqs. (25a) and (25b). Thus the initial electric field can be obtained by solving Poisson Eq. (25a) and adding any divergence-free (e.g., external) electric field. For the magnetic field, one can start either from a

<footnote>$^{25}$ At least in the sense of hydrodynamic equilibrium.</footnote>

<footnote>$^{26}$ The loading of a species with drifting Maxwell-Jüttner distribution in Particle-In-Cell codes should be handled with some care, as discussed e.g., in Refs. (Melzani et al. 2013; Zenitani 2015).</footnote>
<figure><img src="image_31.png" /><figcaption>**Fig. 19** Schematic presentation of the numerical procedure used to solve the Maxwell–Vlasov system of equations. *If the computation of the current densities onto the simulation grid is done in such a way that charge is conserved, and considering that the initial electromagnetic fields satisfy Eqs. (25a) and (25b), solving Maxwell–Ampère (Eq. (25d)) and Maxwell–Faraday (Eq. 25c) is sufficient to ensure that Eqs. (25a) and (25b) remain satisfied (within the machine precision) at all times. **Direct integration here refers to advancing the distribution function on a grid in phase-space (x, p), as further discussed in Sect. 3.4*</figcaption></figure>

zero magnetic field or any non-zero divergence-free magnetic field that will act as an external field applied to the system.

One then enters the time loop of the numerical solver. This time loop consists in advancing the various quantities (defined on the simulation grid) from a timestep $n$ (time $t_n = n \Delta t$) to the next timestep $n+1$ (time $t_{n+1} = t_n + \Delta t$). Various methods are available to do so, some of which rely on defining different quantities at either integer or half-integer timesteps to ensure a centering of the numerical time derivatives. For the sake of simplicity, we will not account for this subtlety here.

The first step in the time-loop consists in advancing the distribution functions for all species $s$. Knowing the electromagnetic field at timestep $n$, the distribution functions at timestep $n+1$ are computed either by direct integration (this is the case in Vlasov codes, see Sect. 3.4) or by advancing so-called macro-particles which are in effect discrete element of the distribution function (this is the case in PIC codes, as described in Sect. 3.3).

The updated distribution functions are then used to compute, onto the simulation grid, the updated charge and current densities (step 2).

These densities are then used, in a third step, to advance the electromagnetic fields from time step $n$ to $n+1$. If the current density deposition onto the grid (step 2) conserves the charge,<fnref n="27" /> solving Maxwell–Ampère and Maxwell–Faraday Eqs. (25d) and (25c), respectively, is sufficient to ensure that the electromagnetic fields remain

<footnote n="27">²⁷ In the sense that it satisfies, at the machine precision, the continuity equation $\partial_t \rho + \nabla \cdot \mathbf{J} = 0$.</footnote>
Maxwell-consistent. To solve these equations, different Maxwell solvers are available (and are discussed in Sect. 3.2.2).

With the fields updated, the loop gets back to step 1 and is run as long as required to reach the final timestep of the simulation.

### 3.2.2 Brief discussion of Maxwell solvers

Various numerical methods (so-called solvers) can be used to solve Maxwell's equations. Here we briefly introduce two methods that are most popular in plasma simulation.

The Finite-Difference Time-Domain (FDTD) method is a time-honoured approach to solve Maxwell's equations (Taflove 2005). It relies on a finite-difference discretization of the partial time derivatives and curl operators in Maxwell–Ampère and Maxwell–Faraday's equations. Most important is that all differential operators are centered. Centering in space requires the use of a staggered grid, the so-called Yee grid, different components of the electric and magnetic fields being defined at different positions in space (onto the grid). Centering in time requires that the electric and magnetic fields are advanced in a leap-frog way, e.g., solving first Maxwell–Ampère equation to advance the electric field then using the updated electric field to solve Maxwell–Faraday's equation and advance the magnetic field. A major advantage of the FDTD method stems from its local nature, which allows for its easy and effective (scalable) implementation in massively parallel environments. A major drawback of the method is that it is subject to numerical dispersion, the numerical electromagnetic wave propagating with velocities potentially smaller than $c$ (Birdsall and Langdon 1985; Nuter et al. 2014). This effect is in part responsible for the spurious numerical Cherenkov instability, see Sect. 3.3.3.

Pseudo-Spectral methods, on the other hand, allow to solve Maxwell's equations with an extraordinary level of precision and correctly capture the dispersion relation of electromagnetic waves. They consist of advancing the electromagnetic fields in Fourier space (for the spatial coordinates) while relying on an (explicit) finite-difference for the time derivatives (Liu 1997; Vay et al. 2013).

The increased precision allowed by pseudo-spectral methods however comes with the cost of global communications associated with the use of Fourier transforms over the entire simulation domain. These global communications have been a major impediment to the adoption of pseudo-spectral methods in massively parallel environments. Recently, Vay et al. (2013) have proposed a domain-decomposition method that allows for the efficient parallelization of pseudo-spectral solvers. This method takes advantage of the finiteness of the speed of light and relies only on local (over subdomains much smaller than the entire simulation domain) fast Fourier transforms and communications between neighbouring subdomains. Vincenti and Vay (2018) have demonstrated that this method may allow for unprecedented scalability of pseudo-spectral solvers over tens to hundreds of thousands of computing elements (cores).
## 3.3 Particle-In-Cell codes

The Particle-In-Cell (PIC) method was introduced in the mid-1950s by Harlow and collaborators to solve fluid dynamics problems (see Harlow 2004 and references therein). Following the pioneering works of Buneman (1959), Dawson (1962), Birdsall and Fuss (1969) and Langdon and Birdsall (1970) (see Dawson 1983 and Verboncoeur 2005 for a history of the development of PIC codes), PIC codes have become a central tool for plasma simulation (Birdsall and Langdon 1985). Indeed, the simplicity of the PIC method together with the possibility to implement it efficiently in a massively parallel environment have established PIC codes as the most popular tool for the kinetic simulation of plasmas.

### 3.3.1 Method

The PIC method differs from the Vlasov-code approach in the way it solves the Vlasov equation, and by extension, the way it computes the current densities on the simulation grid. In PIC codes, the distribution function $f_s$ is approximated as a sum over N macro-particles:

$$f_s(t, \mathbf{x}, \mathbf{p}) = \sum_{p=1}^{N} w_p S(\mathbf{x} - \mathbf{x}_p(t)) \delta(\mathbf{p} - \mathbf{p}_p(t)), \quad (27)$$

where $\delta(\mathbf{p})$ is the Dirac delta-distribution, $S(\mathbf{x})$ is the so-called particle shape-function, and $w_p, \mathbf{x}_p(t)$ and $\mathbf{p}_p(t)$ are the $p^{th}$ particle numerical weight, position and momentum, respectively. The macro-particles can be regarded as *walkers* in Monte-Carlo simulations, and the PIC method as a Monte-Carlo procedure for solving the Vlasov equation (Lapeyre et al. 2003). The initial state of each species of the plasma is obtained from a random sampling of the distribution function $f_s$ at time $t = 0$, and the Vlasov equation is then solved following the macro-particles/walkers motion through the influence of the collective electromagnetic fields.

Indeed, introducing the discretized distribution function (27) in Vlasov equation (23), one can show (see, e.g., Derouillat et al. 2018) that solving Vlasov equation reduces to solving, for all macro-particles $p$, their equations of motion:

$$\frac{d\mathbf{p}_p}{dt} = \frac{q_s}{m_s} \left( \mathbf{E}_p + \frac{\mathbf{v}_p}{c} \times \mathbf{B}_p \right), \quad (28a)$$

$$\frac{d\mathbf{x}_p}{dt} = \mathbf{v}_p = \frac{\mathbf{p}_p}{m_s \gamma_p}, \quad (28b)$$

with $\gamma_p = \sqrt{1 + \mathbf{p}_p^2/(m_s c)^2}$ the $p^{th}$ macro-particle Lorentz factor and where we have introduced the electric and magnetic fields seen by the macro-particle:
$$ \mathbf{E}_p = \int d^3x S(\mathbf{x} - \mathbf{x}_p) \mathbf{E}(\mathbf{x}), \qquad (29a) $$

$$ \mathbf{B}_p = \int d^3x S(\mathbf{x} - \mathbf{x}_p) \mathbf{B}(\mathbf{x}). \qquad (29b) $$

Correspondingly, the species charge and current densities on the simulation grid can be obtained by direct deposition onto the grid. Yet, such a direct deposition would not in general satisfy the charge conservation equation, and the electric fields then would required to be corrected to ensure that they verify the Poisson equation (Mardahl and Verboncoeur 1997). Some current deposition strategies that conserve the charge have, however, been proposed. A popular charge conserving deposition scheme has been proposed by Esirkepov (2001) for PIC code relying on the FDTD Maxwell solver. The macro-particle equations of motion Eq. 28 are most commonly solved using Boris pusher (Birdsall and Langdon 1985). It is a second-order leap-frog integrator where the updated particle momentum is computed knowing the electromagnetic fields at the position of the macro particle as:

$$ \mathbf{p}_p^{(n+\frac{1}{2})} = \mathbf{p}_p^{(n-\frac{1}{2})} + \frac{q_s}{m_s} \Delta t \left[ \mathbf{E}_p^{(n)} + \frac{\mathbf{v}_p^{(n+\frac{1}{2})} + \mathbf{v}_p^{(n-\frac{1}{2})}}{2c} \times \mathbf{B}_p^{(n)} \right], \qquad (30) $$

and the updated particle position is computed as:

$$ \mathbf{x}_p^{(n+1)} = \mathbf{x}_p^{(n)} + \Delta t \frac{\mathbf{p}_p^{(n+\frac{1}{2})}}{\gamma_p}, \qquad (31) $$

Several alternative solvers where recently developed (e.g., Vay 2008; Higuera and Cary 2017) presenting sometimes better accuracy than the traditional Boris algorithm.

### 3.3.2 Additional physics modules

In their most basic implementation (detailed above), PIC codes describe collisionless plasmas through the self-consistent evolution of the particle distribution functions and collective (macroscopic) electromagnetic fields. Additional physics modules can easily be implemented in PIC codes to account for additional processes. Here we provide some references for some of this processes: field ionization (Nuter et al. 2011), collisions and collisional ionization (Nanbu 1997; Nanbu and Yonemura 1998; Pérez et al. 2012), high-energy photon (synchrotron or inverse Compton) emission and its back-reaction (Duclous et al. 2011; Lobet et al. 2016; Niel et al. 2018), pair production in a strong electromagnetic [Breit-Wheeler process (Duclous et al. 2011; Lobet et al. 2016)] or Coulomb [Trident and Bethe-Heitler processes (Martinez 2018)]. These later processes are of outmost importance for extreme plasma physics as will soon be encountered on multi-petawatt laser facilities (see, e.g., Sect. 2.6.1), but also at play in the most extreme astrophysical environments around, e.g., neutron stars and black holes (Uzdensky et al. 2019).
### 3.3.3 Stability issues: relativistic flows

It is worth mentioning that there is one important numerical issue when one deals with relativistic flows in PIC simulations: the *spurious Cherenkov instability* (e.g., Godfrey 1974). This instability results from the resonance of the light-wave modes with the streaming beam when electromagnetic fields are defined on a discrete Eulerian grid. As the numerical light-wave mode is affected by the finite-difference scheme, especially at high-k (Birdsall and Langdon 1985), this resonance is nonphysical. The instability appears as well in spectral codes but has a different signature (Godfrey and Vay 2015). It is practically very difficult to avoid in long-term simulations of relativistic flows or shocks. Nevertheless various methods have been proposed that allow to mitigate, or at least delay the onset of the instability. Some of these methods rely on digital filtering of electromagnetic fields and/or current densities (Greenwood et al. 2004; Vay et al. 2011); modifying the numerical stencil of the FDTD solver (Lehe et al. 2013; Grassi 2017) or upgrading to a semi-implicit scheme (Pukhov 2019); special patching of the most unstable modes (e.g., Godfrey and Vay 2014; Li et al. 2017); artificially increasing the speed of light in the Maxwell solver (Nuter and Tikhonchuk 2016) or solving the PIC equations in Galilean coordinates (Lehe et al. 2016). Yet, there is no definitive solution to remove it completely. Even if the most unstable modes are ‘cleaned’ or stabilized, the difficulty arises in long term evolution from coupling of the secondary aliasing modes with low-wavenumber oblique modes of the streaming plasma, that is very hard to remove without touching important physical scales (Dieckmann et al. 2006). Despite this difficulty, several studies managed to push simulations beyond $10^4 \omega_{pi}^{-1}$ allowing to extract important results from simulations (see Sect. 4).

### 3.3.4 Examples

Various PIC codes are today available and used for astrophysics or space plasma applications.<fnref n="28" /> Some of these codes are freely distributed under free-software licenses, this is the case of EPOCH<fnref n="29" /> (Arber et al. 2015) PICCANTE<fnref n="30" /> (Svattoni et al. 2015), SMILEI,<fnref n="31" /> (Derouillat et al. 2018) TRISTAN-MP<fnref n="32" /> (Spitkovsky 2005), and ZELTRON<fnref n="33" /> (Cerutti et al. 2013). Among other proprietary codes used for astrophysics applications are A-PART (Melzani et al. 2013), CALDER (Lefebvre et al. 2003), OSIRIS (Fonseca et al. 2002), and PHOTON-PLASMA (Haugbølle et al. 2013). Finally, other PIC codes rely on more advanced numerical schemes; e.g., the implicit code iPIC3D (Markidis et al. 2010) or the SLURM code for modeling magnetized fluids or plasmas (Olshevsky et al. 2019).

<footnote n="28">²⁸ We restrict our presentation to electromagnetic PIC codes that have been applied to astrophysics and/or space plasma physics studies.</footnote>

<footnote n="29">²⁹ https://gitlab.com/arm-hpc/packages/wikis/packages/EPOCH.</footnote>

<footnote n="30">³⁰ https://github.com/ALaDyn/piccante.</footnote>

<footnote n="31">³¹ www.maisondelasimulation.fr/smilei.</footnote>

<footnote n="32">³² https://github.com/ntoles/tristan-mp-pitp.</footnote>

<footnote n="33">³³ http://ipag-old.osug.fr/~ceruttbe/Zeltron/index.html.</footnote>
## 3.4 Vlasov-Maxwell codes

### 3.4.1 PIC codes versus Vlasov codes

Particle-In-Cell codes reduce the problem of solving Vlasov equation to solving the (ordinary differential) equations of motion of many macro-particles. It follows that the main advantages of the PIC method are its conceptual simplicity, its robustness and easy implementation on (massively) parallel super computers. The simplicity of the PIC method also allows for PIC codes to be multi-purpose simulation tools: a single PIC code can address various problems from basics plasma physics, astrophysics studies to the modelling of laser-plasma experiments.

However, due to the introduction of a finite number of macro-particles, PIC simulation suffers from the highly exaggerated level of noise. This well known short-coming of the PIC method makes it less adapted to treating problems for which regions of phase-space where the distribution function assumes small values (e.g., in its tail) can impact the physics at play.

In contrast, Vlasov codes which directly integrate the (partial differential) Vlasov equation on a grid in phase-space are virtually noise-free, and are thus an interesting alternative to PIC codes whenever low noise simulation is required.<fnref n="34" /> The fine description allowed by Vlasov codes however comes with the cost of increased numerical complexity. As a result, Vlasov codes are in general much less multi-purpose than PIC codes, and are usually developed to tackle a definite class of problems.

### 3.4.2 The problem of filamentation in phase-space

One impediment in the development and adoption of Vlasov codes is their computational cost and memory requirement when dealing with all 6 dimensions of phase-space. This problem can however be mitigated by reducing the number of dimensions e.g., by relying on conservation laws and symmetries of the problem (see, e.g., Manfredi et al. 1995; Feix and Bertrand 2005). It also becomes less exacting with the fast development of modern high-performance computing.

A more stringent limitation to the development of Vlasov codes stems from the numerical effort necessary to *directly* solve the Vlasov equation, and to the problem of filamentation in phase-space in particular. Indeed, the time evolution of the distribution function in the Vlasov equation is associated with its breaking—filamentation —into increasingly small structures in phase-space, and thus to strong gradients of the distribution function. When discretizing the distribution function onto a grid with finite resolution, handling these gradients becomes numerically inaccurate, and can lead to spurious oscillations, numerical instabilities and inaccurate rendering of conserved quantities (e.g., non-positive distribution functions). Dealing with this issue greatly contributes to the numerical complexity behind Vlasov codes’ development.

<footnote n="34">³⁴ A theoretical discussion on the relative efficiency of the PIC and (direct) Vlasov approaches to treat a given problem is presented in (Feix and Bertrand 2005).</footnote>
### 3.4.3 Example of different methods for electromagnetic Vlasov codes

The numerical complexity of Vlasov codes has—since the seminal work by Cheng and Knorr (1976) introducing the time-splitting technique<fnref n="35" />—led to the development of various techniques to *directly* integrate the Vlasov equation onto a grid in phase-space. It is thus beyond the scope of this brief review to detail these techniques and we here restrict our presentation to some electromagnetic Vlasov codes and their applications. Review articles by (Filbet and Sonnendrücker 2003; Büchner 2007; Ghizzo et al. 2009; Palmroth et al. 2018) discuss various techniques, that range from finite-volume type methods (Fijalkow 1999; Filbet et al. 2001), to semi-Lagrangian methods (Sonnendrücker et al. 1999), and spectral methods (Klimas 1987).

Each method has its own advantages and limitations, and Vlasov codes are usually designed to tackle a specific class of physical problems. Ghizzo et al. (1990) for instance developed a relativistic electromagnetic 1D Vlasov code to study stimulated Raman scattering; while Shoucri et al. (2015) considered the problem of stimulated Brillouin scattering. These codes actually used the conservation of canonical momentum to reduce the number of dimension in velocity/momentum-space.<fnref n="36" /> A somewhat similar approach was used to study the breaking of a relativistic Langmuir wave (Grassi et al. 2014) as well as laser-driven (electrostatic) shock acceleration of ions and ion turbulence (Grassi et al. 2016).

A (non-relativistic) Eulerian Vlasov–Maxwell solver was developed by Mangeney et al. (2002). It was applied to various studies ranging from wave propagation in magnetized plasmas (Califano and Lontano 2003), to the study of the nonlinear kinetic regime of the Weibel instability (Califano et al. 2002). An off-spring of this solver is the hybrid (kinetic ions, fluid electrons) Vlasov code (Valentini et al. 2007) used in particular to tackle turbulence studies in either 2D3V (see, e.g., Cerri et al. 2017) and 3D3V (see, e.g., Cerri et al. 2018) geometries.

Another Eulerian Vlasov–Maxwell model was developed by Umeda et al. (2009) and applied to the study of various instabilities, such as the Kelvin–Helmholtz instability (Umeda et al. 2014) or the collisionless Rayleigh–Taylor instability (Umeda and Wada 2016).

The semi-Lagrangian method introduced by Sonnendrücker et al. (1999) (see also Crouseilles et al. 2010) has also led to a new kind of Vlasov codes. As an example, a relativistic semi-Lagrangian Vlasov–Maxwell solver (VLEM) was recently developed by Sarrat et al. (2017). It was used to tackle problems related to streaming instabilities in plasmas, such as the current Weibel-filamentation and two-stream instabilities; and operates in 1D3V, 2D2V and 2D3V geometries.

## 3.5 Hybrid methods

In this approach thermal electrons are taken to be a massless, neutralizing and are treated as a magnetized fluid. Ions (thermal or even non-thermal) are treated using a

<footnote n="35">³⁵ The time-splitting technique separates advection in (real) space and velocity-space.</footnote>

<footnote n="36">³⁶ Other methods have been developed that take advantage of the existence of canonical invariants to solve the Vlasov equation; see e.g., Liseikina et al. (2004), Inglebert et al. (2011).</footnote>
PIC approach. The advantage of this method is to eliminate Debye-scale physics while still catching microscopic phenomena.

In hybrid codes, ion positions are advanced using the Boris pusher as in PIC codes (see Sect. 3.3.1). Electron dynamics is the one of a massless fluid then

$$m_e n_e \frac{d\mathbf{v}_e}{dt} = 0 = -e n_e \left( \mathbf{E} + \frac{\mathbf{v}_e}{c} \times \mathbf{B} - \nabla \cdot \bar{\mathbf{P}}_e \right). \quad (32)$$

This combined with the Ampère law for a non-relativistic flow, hence neglecting the displacement current leads to an equation for the electric field

$$\mathbf{E} \approx -\frac{\mathbf{v}_i}{c} \times \mathbf{B} - \frac{1}{e n_e} \nabla \cdot \bar{\mathbf{P}}_e - \frac{1}{4\pi q_i n_i} (\nabla \times \mathbf{B}) \times \mathbf{B}, \quad (33)$$

where $\bar{\mathbf{P}}_e$ is the electron pressure rank 2 tensor. This method will not be reviewed here, the interested reader is invited to read recent references on the subject: Lipatov (2002), Kunz et al. (2014).

We note here the case of the DHYBRID code (Gargaté et al. 2007). This code is explicit fully parallelized code and it uses MPI. DHYBRID solves the dynamics of non-thermal particles based on a PIC approach. The code has been used in the context of particle acceleration and transport at collisionless shocks, some of its results are presented in Sect. 4.1.

## 3.6 Solving Fokker-Planck problems

The Fokker-Planck equation (FPE) is one of the most important equation in kinetic physics. It describes the evolution in the phase space of the particle distribution function $f(\mathbf{r}, \mathbf{p}, t)$ under the effect of a diffusion process with small increments in which initial conditions are lost (a.k.a. a Markov process). In this review we are interested in collisionless plasmas, in that case, particle diffusion results from the process of scattering off plasma waves. However, note that FPEs are also well studied in the context of collisional plasmas. We refer the interested reader to Wang et al. (2008) for the description of numerical treatments of the collision operator. As is concerning high-energy particles, the FPE describes processes which develop over scales explored by these particles, it is also adapted to the study of macroscopic processes in astrophysical plasmas detailed in Sect. 5. The interested reader can advantageously consult Risken (1989) for an overview of the properties of the FPE.

For a system of energetic particles in a magnetic field oriented along the z axis, we can write the FPE as (Schlickeiser 2002):

$$\partial_t f + v\mu\partial_z f - \epsilon\Omega_s\partial_\phi f = \frac{1}{p^2}\partial_x \left[ p^2 (D_{xx}\partial_x f + af) \right] + q(\mathbf{r}, \mathbf{p}, t) \quad (34)$$

where the diffusion process runs over the variables: x, y, z, p, $\mu$, $\phi$, hence we have 25 diffusion coefficients $D_{xx}$,37 and $\mu$ and $\phi$ are the particle pitch-angle cosine and

<footnote n="37">37 The Cartesian coordinates x, y, z mark here the position of the particle’s guiding center.</footnote>
azimuthal gyration angle respectively. Here, particles of charge q and mass m are relativistic (with speeds v ≃ c) and gyrate around a magnetic field of strength B with a gyrofrequency Ωs ≃ c/rL. We note ϵ = q/sgn(q). The term a(p, r, t) describes the momentum change of the particle either due to loss or acceleration and q(r, p, t) represents particle injection and/or escape. Although it should be kept in mind that the FPE is deduced from the more general Vlasov equation, we focus below on numerical solutions of this equation. Often, in the context of CR physics, the FPE is not directly solved but rather the convection-diffusion equation (CDE). The CDE results from the former by an averaging procedure over φ and µ in the case fast scattering processes build a gyrotropic and an isotropic distribution.

### 3.6.1 The Fokker-Planck equation

We start by studying 1D diffusion problems as is the case for stochastic acceleration. In that case the FPE can be simplified as

$$ \partial_t F(p, t) = \frac{1}{p^2} \partial_p \left[ p^2 (D_{pp} \partial_p F + a(p)F) \right] - \frac{F}{\tau_{\text{esc}}} + Q(p, t). \quad (35) $$

Here, the particle distribution $F(p, t) = \int f(r, p, t)d^3rd\mu d\phi$ is averaged over the space volume and is assumed to be isotropic (it fulfills the diffusion-convection limit) and diffusive escape is treated by the means of an escape timescale $\tau_{\text{esc}}(p)$, the loss/gain $a(p, t)$ term is also averaged. This equation can be solved using finite difference schemes (Park and Petrosian 1996).

**Boundary conditions** As stated by Park and Petrosian (1996) any boundary condition which is a linear combination of F and $F'(p) = \partial_p F$ is viable for Eq. (35) if the points $p_1$ and $p_2$ at which they are taken fulfill $0 < p_1 < p < p_2 < \infty$. So we end up with two types of boundary conditions either with no particle $F(p_1) = F(p_2) = 0$ or with no flux $\phi(p_1) = \phi(p_2) = 0$ at the boundaries, where $\phi(p) = -(p^2 D_{pp}F + a(p)F)$. The choice of one condition with respect to the other depends on the specific problem under investigation. A drawback of the no-particle condition is that it does not respect the particle number conservation.

**Numerical schemes** A simple way to solve Eq. (35) is to use an explicit finite difference method (FDM) with fluxes evaluated at grid midpoints, namely

$$ \frac{F_{j+1}^{n+1} - F_j^n}{\Delta t} = -\frac{1}{p_j^2} \left( \frac{\phi_{j+1/2}^n - \phi_{j-1/2}^n}{\Delta p} \right) - \frac{F_j^n}{\tau_{\text{esc}}(p_j)} + Q^n(p_j). \quad (36) $$

Time is discretized as $\Delta t = t_{n+1} - t_n$ and momentum is discretized following a constant logarithmic mesh where $\Delta p_j/p_j = \text{constant}$. We write $\Delta p_j = (p_{j+1} - p_{j-1})/2$. The fluxes are calculated at midpoints defined as $p_{j+1/2} = (p_{j+1} + p_j)/2$. The coefficients entering in the flux calculation are evaluated as, e.g., $a_{j+1/2} = (a(p_{j+1}) + a(p_j))/2$ instead of a direct evaluation at $p_{j+1/2}$. For an explicit scheme the CFL condition (see Sect. 5.2.1) $\Delta t/\Delta p_j^2 < p_j^2/D_{pp,j}$ usually produces prohibitively small time steps.
<figure><img src="image_32.png" /><figcaption>**Fig. 20** Time-dependent solutions of the FP problem $\partial_t F = \partial_p (p^3 \partial_p F - p^2 F) - F + \delta(p - p_0)\delta(t)$ over a logarithmically-spaced mesh in momentum. The injection momentum is $p_0$. Three different numerical boundaries are compared: $p_1 = 10^{-2}$ and $p_2 = 10^2$ (short-dashed lines); $p_1 = 10^{-2.5}$ and $p_2 = 10^{2.5}$ (medium-dashed lines); $p_1 = 10^{-3}$ and $p_2 = 10^3$ (long-dashed lines) The steady state numerical solutions were obtained at $t = 10$ in normalized units. Image reproduced with permission from Park and Petrosian (1996), copyright by AAS</figcaption></figure>

Semi-implicit and implicit methods can be used to circumvent this problem; they are obtained by changing $n$ to $n + 1/2$ and $n + 1$ in the RHS of Eq. (36) respectively. These methods lead to the derivation of a tridiagonal system of equations that can easily be solved once the boundary conditions are selected. For a given class of method, schemes then differ by the way the flux is calculated. One efficient implicit method is due to Chang and Cooper (1970) and a well-known semi-implicit calculation is the Crank-Nicholson method (see Press et al. 2002). These methods are second order in time and second order in momentum for a uniform grid and first order in momentum for a non-uniform grid. Park and Petrosian (1996) show that the no-flux condition and the implicit Chang-Cooper scheme ensure positive solutions of 1D FP problems contrary to the Crank-Nicholson method (see an example of a solution of a FP problem in Fig. 20). While accounting for losses in $\phi(p)$ it is useful to adapt the time step to the dominant loss timescale. Donnert and Brunetti (2014) use a time step $\Delta t = 1/2 \min(t_{\text{loss}}(p_j))$, where $t_{\text{loss}}(p) = a(p)/p$.

**A note on particle transport and stochastic acceleration in hot plasmas** FDM are widely used in CR physics but they are also used to solve radiative transfer problems in hot plasmas which develop in corona or in jets associated with compact objects. The radiative transfer in accretion disk corona can be treated in 1D assuming some particular geometry (usually slab-type or spherical) for the source of high-energy particles which allows to derive an escape probability and hence a simple expression for $\tau_{\text{esc}}$ in Eq. (35). In this approach the coupled system of electron-positron plasma and its associated photon field can be described by a set of FPEs. However, a major

<figure><img src="image_33.png" /></figure>
difficulty to simulate such plasma systems is that they involve non-local processes in momentum.<fnref n="38" /> This is for instance the case for the Compton scattering process (Nayakshin and Melia 1998). In that case, the fluxes are expressed in terms of integrals over lepton and photon populations (Belmont et al. 2008; Vurm and Poutanen 2009; Marcowith et al. 2013).<fnref n="39" /> Aside from Compton scattering, integrals also result from the calculation of other processes: pair production and annihilation, Coulomb losses, synchrotron losses. These codes solve the diffusion problem usually using Chang–Cooper-type methods. However, as noticed by Belmont et al. (2008), the radiative transfer problem requires a high accuracy in momentum to preserve a high level of particle number and energy conservation. The Chang–Cooper method which is only first order accurate on non-uniform grid needs to be modified using both grid center and faces and calculating the momentum derivatives as $\partial_p F = (F_{j+1/2} - F_{j-1/2})/\Delta x_j$ and $\partial_p^2 F = (F_{j+1} - F_j)/\Delta p_{j+1/2} - (F_j - F_{j-1})/\Delta p_{j-1/2}$. The integral parts have to be calculated using specific treatments. First, the different elements of the cross sections are stored and then interpolated during the course of the runs. Then, boundary conditions are different for the FP and the transfer parts. The transfer part has a wall-type boundary condition which includes a modification of the differential cross section [see Belmont et al. (2008) for details]. Finally, integrals used to calculate the Compton process can be treated differently depending on the photon energy with respect to the electron energy from a continuous process at low energy leading to a derivative term and to a full integral calculation in the Klein–Nishina limit. Compton scattering of photons involves the same kind of treatment and is applied depending on the electron energy [see (Vurm and Poutanen 2009) for details].

### Multi-variable FPE
Astrophysical or space plasmas usually involve multi-dimensional diffusion-advection processes. The study of CR propagation in the Milky Way requires to account for several complex effects: CR spallation reactions, radioactive decay, anisotropic diffusion with respect to the background magnetic field direction, etc. Specific numerical tools have been developed to handle this complexity.<fnref n="40" /> Most of CR transport codes use multi-dimensional finite difference methods, this is the case for GALPROP and DRAGON which adopt a Crank–Nicholson method. In multi-dimensional problems this method leads to a non-tridiagonal system of equations to solve. It is solved usually adopting an iterative procedure like the Gauss–Seidel relaxation method (Press et al. 2002). The integration is adapted to the specific diffusion problem by starting from a large time step and reducing it as the stationary solution is reached (Strong and Moskalenko 1998). Both GALPROP and DRAGON codes also use an operator splitting technique to handle multi-dimensional diffusion problems (Press et al. 2002). The technique of operator splitting consists in splitting the time integration in Eq. (35) or

<footnote n="38">³⁸ Some processes are also in principle non-local in space, this is the case for instance of the inhomogeneous synchro-Compton effect in jets, i.e., the Inverse Compton scattering of low energy photons generated by synchrotron radiation by a population of relativistic electrons (see Ghisellini et al. 1985).</footnote>

<footnote n="39">³⁹ A version of such radiative transfer codes adapted to GRBs can be found in Vurm et al. (2011).</footnote>

<footnote n="40">⁴⁰ We point towards the corresponding code websites: DRAGON: https://github.com/cosmicrays, GALPROP: https://galprop.stanford.edu/code.php, PICARD: http://astro-staff.uibk.ac.at/~kissmrbu/Picard.html Let us also mention the semi-analytical tool USINE, https://dmaurin.gitlab.io/USINE/, where the FPE is solved using a path integral method.</footnote>
its multi-dimensional generalization into a succession of simpler operations involving N different operators $L_i$, such that

$$ \partial_t F(p, t) = \sum_{i=1}^{N} L_i F(p, t). \qquad (37) $$

Each operator contributes to move the solution from $F^n$ to $F^{n+1}$ as $F^{n+1} = \prod_{i=1}^{N} L_i F^n$ where each finite difference operator solves a part of the numerical problem. One difficulty with this method is that operator actions do not commute, hence one usually has to proceed with trials with a guess of the correct solution to select the correct operator ordering. The DRAGON code designed in cylindrical coordinates uses a series of operators associated to each relevant transport process. For instance the operator associated with diffusion along galactic vertical height $z$ is $L_z = D_{zz}\partial_z^2 F(z, r, p, t) + \partial_z D_{zz}\partial_z F(z, r, p, t)$. Similarly other operators are derived for diffusion along other space variable $r$ (the galactic radius), momentum loss or advection (Evoli et al. 2017). Then each derivative is treated using a Crank–Nicholson scheme. The PICARD code uses a different numerical approach as it first solves a stationary problem and also as the momentum evolution is treated using an integration instead of a FDM (Kissmann 2014).

Multi-dimensional numerical solutions of FP problems is an active research field with a rich variety of solvers based on three main approaches: finite difference methods as we just discussed (FDM), finite element methods (FEM) or path integrals techniques. For most of them they still wait to be applied in the context of astrophysical or space plasma research.

### 3.6.2 Stochastic differential equations

Stochastic differential equations or SDEs are a very efficient way to solve complex multi-dimensional Fokker–Planck problems with simple numerical schemes, although SDE schemes can become themselves rather complex. We invite the interested reader to consult some monographs cited in Strauss and Effenberger (2017). These authors provide an overview of the use of SDE in the fields of DSA, CR transport in the ISM and space plasmas. The interested reader can consult this complete review to what concerns space plasmas problems. Below we bring a complementary discussion on the use of SDE in the context of shock acceleration. The intrinsic idea behind SDE is to derive a set of equations of motion which reproduce the random walk in each of the stochastic variables which describe the phase space evolution of a particle. Kruells and Achterberg (1994) demonstrate the equivalence between a FPE and a set of SDEs. The simplest SDE scheme is the Ito first order explicit scheme. As an example let us write the SDE for a 1D random walk in a space x direction of a particle represented by a diffusion coefficient $D(x, t)$. Let us assume also that the particle is advected with a speed $u(x, t)$. The first order forward explicit Ito scheme then writes the increment of the position of the particle within a time step $\Delta t$ as
$$ \Delta x = \left( u(x, t) + \frac{\partial D(x, t)}{\partial x} \right) \Delta t + \xi_x \sqrt{2D(x, t)\Delta t} = \Delta x_{\text{adv}} + \xi_x \Delta x_{\text{diff}}. \quad (38) $$

We note $V(x, t) = u(x, t) + \partial D(x, t)/\partial x$. This equation shows that the particle path has two terms, the first term is deterministic and reproduces a forward Euler increment due to advection. The second term is stochastic and describes a diffusion as $\Delta x \propto (\Delta t)^{1/2}$. The term $\xi_x$ is a random variable usually sampled over a Gaussian distribution with 0 mean and variance 1. The particle distribution can then be reconstructed using a large number of particles. The method is simple; however, it can suffer from noise in parts of the phase space sampled by only a few particles. The latter issue can be partly handled using a particle-splitting scheme (Yamazaki et al. 2015) where the weight attributed to a particle is split over several particles when reaching a region of the phase space with low statistics, as can be the case in the energy space if the particle distribution has an exponential cut-off.

More generally, the diffusive term in Eq. (38) is a Wiener processes $W(t, x)$ which models the Brownian motion of a particle in an homogeneous medium, we write $dW(x, t)/dt = \xi_x$. More complex SDE schemes can be interesting to use if necessary, like schemes backward in time in order to start from a known distribution and reconstruct the particle injection at sources. This way has the advantage to improve statistics if we want to have information at a particular location, corresponding for instance to a satellite. Higher order schemes in space and time are possible by doing a Taylor expansion of both advection and stochastic parts of Eq. (38). Schemes stable in time can be obtained by searching advection and diffusive terms at a time $t' = t + \theta \Delta t$, where $\theta$ is to be taken between 0 and 1 (Smith and Gardiner 1989).

**DSA with SDEs** The study of shock acceleration using SDEs requires some care in fixing the time step $\Delta t$ (Kruells and Achterberg 1994). In shock acceleration studies the shock front is usually obtained from a fluid code, so has a finite width $\Delta x_{sh}$ traced by a few grid cells. The condition over the time step to describe the DSA process properly is then $\Delta X_{adv} < \Delta x_{sh} < \Delta X_{diff}$. The first inequality allows particles to stay around the shock to get accelerated whereas the second inequality allows the particle to sample the up- and downstream media correctly. However, if the diffusion coefficient is an increasing function of the particle energy, i.e., $\partial D(x, E, t)/\partial E > 0$ it is possible to find a threshold energy $E^*$ for which the condition $\Delta x_{sh} = \Delta X_{diff}(E^*)$ is fulfilled (Casse and Marcowith 2005; Schure et al. 2010). Below $E^*$ the shock acceleration process can not be properly treated. One possibility to address this problem is to sharpen artificially the shock (Casse and Marcowith 2003). This method can be easily handled in 1D (Marcowith and Casse 2010) but is difficult to construct in 2 or 3D as the shock front starts to corrugate. Another difficulty is that at an non-parallel shock, the MHD Rankine-Hugoniot conditions induce a discontinuous diffusion coefficient up- and downstream. Quite generally the diffusion transition at the shock can be decomposed into a continuous component $D_c$ and a jump at the shock front $\Delta D = D_u - D_d$ expressed in terms of the up- and downstream diffusion coefficients. The diffusion coefficient can then be written as

$$ D(x) = D_c(x) + \Delta D \delta(x - x_{sh}), \quad (39) $$
where $x_{sh}$ is the shock position (Marcowith and Casse 2010). Zhang (2000) proposes to account for the discontinuous part using a skewed Brownian motion which introduces an asymmetric shock crossing probability. To proceed we introduce a new variable $\tilde{x} = x\zeta(x)$ where

$$ \zeta(x) = \begin{cases} \epsilon & x < x_{sh} \\ 1/2 & x = x_{sh} \\ (1 - \epsilon) & x > x_{sh}, \end{cases} \qquad (40) $$

with $\epsilon = D_u(x_{sh})/(D_d(x_{sh})+D_u(x_{sh}))$. Achterberg and Schure (2011) propose a more general scheme adapted to shock configuration with strong gradients in the diffusion coefficient. This situation occurs especially upstream, in the shock precursor, in case of strong magnetic field amplification. The scheme involves a second-order accuracy predictor-corrector method [see Section 4 in Achterberg and Schure (2011) for details]. The scheme is however much slower than the simple Ito scheme and it is necessary to switch from one scheme to the other in order to save simulation resources.

One also has to account for the particle increment in energy or momentum at each shock crossing. It is also possible to use an explicit Ito scheme for CR energy or momentum similarly to Eq. (38). If stochastic acceleration can be neglected, Marcowith and Kirk (1999) introduce an implicit scheme:

$$ \Delta \ln(p) = - \left( a_{\text{loss}} p + \frac{1}{3} \frac{du}{dx} \right) \Delta t, \qquad (41) $$

where $a_{loss}$ is a loss rate and the second term accounts for the increase in particle momentum from shock acceleration. The implicit scheme rewrites the particle position with time as a linear interpolation $x = (\Delta x/\Delta t)t$. Eq. (41) has the solution

$$ \ln(p(t')/p(t)) = -\ln(F_I + L_s), \qquad (42) $$

where $F_I = \exp((\Delta V/3)\Delta t/\Delta x)$ gives the momentum increment by DSA and

$$ L_s = a_s \frac{\Delta t}{\Delta x} p \int_{x(t)}^{x(t')} \exp \left( \frac{\Delta V}{3} \frac{\Delta t}{\Delta x} \right) dx', \qquad (43) $$

where $L_s$ accounts for the effect of losses. The increment $\Delta x$ is calculated from the SDE in $x$, which is evaluated at $t' = t + \Delta t$. The new momentum is obtained from Eq. (42) calculated using $x(t') = x(t + \Delta t)$. Figure 21 gives the shock solution for electrons including synchrotron losses.

The SDE method has proven to be very efficient in calculating particle acceleration by DSA at non-relativistic shocks in 1D and even in 2.5D in the context of jets (Casse and Marcowith 2003, 2005). The method can in principle be coupled to MHD solutions by sub-cycling the MHD timestep (see Sect. 5.5). To our knowledge no scheme has yet included CR back-reaction over the thermal plasma, but solutions proposed by other Monte-Carlo models (see Sect. 3.6.3) should be applicable to this particular technique.
<figure><img src="image_34.png" /><figcaption>**Fig. 21** Shock electron distribution including synchrotron losses. Particles are injected at momentum $p_0$ for $a_{s,u} = 10$, $a_{s,d} = 1$ and a compression ratio $r = 4$ compared to the analytical solution of Webb et al. (1984). The solid line shows a solution $f(p) \propto p^{-4}$. Image reproduced with permission from Marcowith and Kirk (1999), copyright by ESO</figcaption></figure>

**Relativistic shocks and SDEs** In the context of relativistic shock two difficulties emerge if one wants to apply the SDE method to the problem of particle acceleration because the shock is moving almost as fast as the particles. First, a simple scheme as in Eq. (38) may lead to a violation of the causality principle, because over the diffusive step the particle speed $\Delta x/\Delta t$ can exceed the speed of light. Second, DSA is based on the diffusive approximation which requires the ratio of the particle speed to the shock speed to be small. Achterberg et al. (2001) evaluate particle acceleration in the shock rest frame but simulate the spatial diffusion process from the pitch-angle scattering process and hence reconstruct particle trajectories ([see also Bednarz and Ostrowski 1998]). These works retrieve the shock particle distribution produced by scattering by an isotropic turbulence $f(p) \propto p^{-4.2}$ (see Sect. 2.4).

### 3.6.3 Simulating shock acceleration using a Monte-Carlo method

Ellison and Eichler (1984), Jones and Ellison (1991) proposed a Monte-Carlo method to simulate particle pitch-angle scattering around non-relativistic shocks. The particle mean free path is assumed to scale as a function of the particle rigidity as $\lambda \propto R^a/\rho$, where $R = pc/q$ is the particle rigidity and $\rho$ is the fluid mass density. The momentum vector follows a random walk which produces a variation of the particle pitch-angle $\delta\alpha$. The scattering is assumed to be elastic and isotropic in the fluid rest frame. The particle are injected upstream from the thermal plasma. The CR pressure is reconstructed at different distances from the shock front; particle-splitting technique is used in order to improve statistics at high energies. The CR pressure term is then included in the Rankine-Hugoniot conditions to account for CR backreaction over the shock solutions. Finally, a far escape boundary is adopted to calculate the escaping energy flux carried by the particles. An example of results can be seen in Fig. 43, compared with two other methods discussed elsewhere in this review.
Recently the technique has been used to study the effect of magnetic field amplification in NLDSA at non-relativistic SNR shocks (Vladimirov et al. 2006), as well as acceleration at relativistic shocks in GRB afterglows (Warren et al. 2015).

# 4 Small and meso-scale numerical particle acceleration studies

The understanding of the initial stages of particle acceleration in astrophysical plasma relies on the non-linear interplay between the particle distribution function and electromagnetic fields. The inherent non-linearity of the process prevents the development of robust analytic models, unless numerical simulations provide basic guidelines of the behavior of the system. This is especially true when particles produce the turbulence responsible for their self-confinement and acceleration around the shock front. The improvement in the computation power during last decades allowed for significant progresses in this field of research using computationally expensive, yet considered as *ab-initio*, PIC simulations. Despite the fact that the astrophysical sources space and time scales are out of reach of PIC simulations, a large number of fundamental questions have been addressed and, sometimes, answers provided using this technique. In this section, we review a number of works that investigate the question of particle acceleration efficiency at shocks (Sect. 4.1) and in magnetic reconnection (Sect. 4.2) processes, using PIC simulations. Before we further proceed, we first introduce some vocabulary concerning the different category of shocks investigated with the help of PIC simulations.

**Collisionless shocks classes** Collisionless shocks are mediated by collective plasma effects (Sagdeev 1966). In this sense, more refined classification is required then for hydrodynamical/MHD shocks (weak, strong, radiative, fast, slow, parallel, oblique, perpendicular) where the shock is expected to be mediated by binary collisions between particles. Seminal (Sagdeev 1966) as well as recent studies (e.g., Stockem et al. 2014b, a; Bret et al. 2014; Ruyer et al. 2016) considered and demonstrated the dominant role of small scale plasma instabilities in forming and mediating collisionless shocks. Different types of instabilities are dominant depending on plasma beta, composition, shock Mach number and upstream magnetic field orientation with respect to the shock propagation direction. This translates into a bestiary of different plasma instabilities at play when describing the shock structure. Commonly, the separation into electrostatic, Weibel-mediated and magnetised shocks is done. Non-relativistic, weakly magnetised and low-Mach number shocks are believed to be mediated by electrostatic effects (two-stream, Buneman instabilities). With increasing Mach number ($M_d \gg 1$) or going into relativistic regime, weakly magnetised shocks are mediated by Weibel-filamentation (Bret et al. 2014; Huntington et al. 2015). Strongly magnetised shocks are typically mediated by coherent magnetic reflection of particles on the shock barrier. In this case the shock width is of the order of ion gyroradius.

## 4.1 Shock acceleration numerical experiments: PIC simulations

The short chronological version is the following. In the pioneering studies 1D3V geometry was adopted because of numerical cost (e.g., Biskamp and Welter 1972). Several

<figure><img src="image_35.png" /></figure>
physical processes were evidenced in this way, such as shock front self-reformation (e.g., Lembege and Dawson 1987; Lembege and Savoini 1992) or positron acceleration through resonant absorption of ion-cyclotron waves (Hoshino et al. 1992). However, this configuration was found to be too restrictive to trigger the Fermi process in most of cases, especially for quasi-perpendicular shocks. Next, multidimensional simulations were long enough to form the shock, but no evidence of first-order Fermi acceleration was found (Frederiksen et al. 2004; Hededal and Nishikawa 2005; Kato 2007; Dieckmann et al. 2008). This raised the question whether shocks can accelerate particles self-consistently or some external source of turbulence was necessary. More recent studies were able to form the shock and follow its propagation long enough to allow several Fermi cycles and produce extended power-law particle distribution self-consistently (e.g., Spitkovsky 2008b; Martins et al. 2009; Sironi and Spitkovsky 2011; Plotnikov et al. 2018; Crumley et al. 2019; Lemoine et al. 2019b). Yet, even in the longest simulations the power-law spans no more than two orders of magnitude in particle energy, reflecting the challenging nature for fully kinetic simulations to reach astrophysical space and timescales.

### 4.1.1 Ultra-relativistic shocks

The ultra-relativistic regime is particularly interesting for several reasons.<fnref n="41" /> (i) Energy gain per cycle is large, $\Delta E/E \simeq 2$ instead of $\Delta E/E = \beta_{sh} \ll 1$ in the non-relativistic regime and (ii) Scattering time must be short, otherwise particles get advected within the downstream flow as the shock front recedes rapidly (the shock front moves away with velocity equal $c/3$ in the frame where the downstream plasma is at rest). These two reasons mean that the build-up of the non-thermal power law is faster in the ultra-relativistic regime than in the non-relativistic case. Hence, one can diagnose whether relativistic shock accelerate particles efficiently or not for a given parameter regime on a timescale of several $10^3 \omega_{pe}^{-1}$ (given upstream flow magnetization $\sigma$, magnetic field inclination with respect to the shock propagation direction and plasma composition). On the other hand, fast and efficient particle acceleration prompted early Monte-Carlo and semi-analytical studies to suggest that relativistic shocks are viable candidates for the acceleration of Ultra High Energy CRs (UHECRs).

For the reasons outlined just before, the demonstration of the first order Fermi operability in PIC simulations was firstly done in the regime of relativistic shocks by Spitkovsky (2008b). In the non-relativistic case it was done several years later, as it requires much longer simulation time (e.g., Kato 2015; Park et al. 2015).

We now turn to the discussion of different studies that used PIC simulations to understand the physics of ultra-relativistic shocks. The studies go from 1D to 3D, deal with different plasma compositions, different magnetic field geometries and magnetization. Table 1 provides a (non-exhaustive) list of such studies and presents some relevant numerical parameters that they used. Below we discuss different type of micro-instabilities which are reviewed in Marcowith et al. (2016), Bret (2009).

<footnote n="41">⁴¹ See also the recent review article by Vanthieghem et al. (2020) on the physics of weakly magnetized shocks.</footnote>
<table><caption>**Table 1** Table of different PIC studies of relativistic shocks, referenced by the authors names and publication dates</caption><thead><tr><th>References</th><th>Compos.</th><th>Dim.</th><th>m<sub>i</sub>/m<sub>e</sub></th><th>σ</th><th>θ<sub>B</sub><sup>o</sup></th></tr></thead><tbody><tr><td>Langdon et al. (1988)</td><td>e<sup>-</sup> - e<sup>+</sup></td><td>1D</td><td>1</td><td>[0.1; 13.3]</td><td>90</td></tr><tr><td>Gallant et al. (1992)</td><td>e<sup>-</sup> - e<sup>+</sup></td><td>1D</td><td>1</td><td>[3 · 10<sup>-5</sup>; 5]</td><td>90</td></tr><tr><td>Hoshino et al. (1992)</td><td>e<sup>-</sup> - e<sup>+</sup> - i</td><td>1D</td><td>20</td><td>[5 · 10<sup>-3</sup>; 0.5]</td><td>90</td></tr><tr><td>Nishikawa et al. (2003)</td><td>e<sup>-</sup> - i</td><td>3D</td><td>20</td><td>0</td><td>-</td></tr><tr><td>Frederiksen et al. (2004)</td><td>e<sup>-</sup> - i</td><td>3D</td><td>16</td><td>0</td><td>-</td></tr><tr><td>Hededral et al. (2004)</td><td>e<sup>-</sup> - i</td><td>3D</td><td>16</td><td>0</td><td>-</td></tr><tr><td>Spitkovsky (2005)</td><td>e<sup>-</sup> - e<sup>+</sup></td><td>3D</td><td>1</td><td>[0; 0.1]</td><td>90</td></tr><tr><td>Lyubarsky (2006)</td><td>e<sup>-</sup> - i</td><td>1D</td><td>50</td><td>3 · 10<sup>-3</sup></td><td>90</td></tr><tr><td>Kato (2007)</td><td>e<sup>-</sup> - e<sup>+</sup></td><td>2D</td><td>1</td><td>0</td><td>-</td></tr><tr><td>Hoshino (2008)</td><td>e<sup>-</sup> - i</td><td>1D</td><td>50</td><td>2 · 10<sup>-3</sup></td><td>90</td></tr><tr><td>Dieckmann et al. (2008)</td><td>e<sup>-</sup> - i</td><td>1D,2D</td><td>400</td><td>2.5 · 10<sup>-3</sup></td><td>≈ 10</td></tr><tr><td>Chang et al. (2008)</td><td>e<sup>-</sup> - e<sup>+</sup></td><td>2D</td><td>1</td><td>0</td><td>-</td></tr><tr><td>Spitkovsky (2008a)</td><td>e<sup>-</sup> - e<sup>+</sup></td><td>2D</td><td>1</td><td>0</td><td>-</td></tr><tr><td>Spitkovsky (2008b)</td><td>e<sup>-</sup> - i</td><td>2D</td><td>≤ 10<sup>3</sup></td><td>0</td><td>-</td></tr><tr><td>Keshet et al. (2009)</td><td>e<sup>-</sup> - e<sup>+</sup></td><td>2D</td><td>1</td><td>0</td><td>-</td></tr><tr><td>Martins et al. (2009)</td><td>e<sup>-</sup> - i</td><td>2D</td><td>32</td><td>0</td><td>-</td></tr><tr><td>Sironi and Spitkovsky (2009)</td><td>e<sup>-</sup> - e<sup>+</sup></td><td>2D</td><td>1</td><td>0.1</td><td>[0; 90]</td></tr><tr><td>Sironi and Spitkovsky (2011)</td><td>e<sup>-</sup> - i</td><td>2D,3D</td><td><b>16</b></td><td>[10<sup>-5</sup>; 0.1]</td><td>[0; 90]</td></tr><tr><td>Haugbølle (2011)</td><td>e<sup>-</sup> - i</td><td>2D,3D</td><td>16</td><td>0</td><td>-</td></tr><tr><td>Sironi et al. (2013)</td><td>e<sup>-</sup> - e<sup>+</sup>/i</td><td>2D,3D</td><td>1& <b>25</b></td><td>[0; 0.1]</td><td>90</td></tr><tr><td>Bret et al. (2013)</td><td>e<sup>-</sup> - e<sup>+</sup></td><td>2D</td><td>1</td><td>0</td><td>-</td></tr><tr><td>Ardaneh et al. (2015)</td><td>e<sup>-</sup> - i</td><td>3D</td><td>16</td><td>0</td><td>-</td></tr><tr><td>Plotnikov et al. (2018)</td><td>e<sup>-</sup> - e<sup>+</sup></td><td>2D</td><td>1</td><td>[0; 5]</td><td>90</td></tr><tr><td>Crumley et al. (2019)</td><td>e<sup>-</sup> - i</td><td>2D</td><td><b>64</b></td><td>0.007</td><td>10 &amp; 55</td></tr></tbody></table>

Mass ratio values indicated with bold font indicate that larger values were also explored

**1D studies** The first exploration of relativistic shocks was motivated by the study of the termination shock physics in Pulsar Wind Nebulae. Langdon et al. (1988) performed 1D PIC simulations of perpendicular magnetized shocks in pair plasma with upstream Lorentz factors $\gamma_0 = 20$ and 40 and Alfvénic Mach numbers going from 24 to 154. The size of the simulation box was about 10 Larmor radii for $\sigma = 0.1$ and about 100 Larmor radii for $\sigma = 13.3$. The shock front formed by magnetic reflection between the incoming and wall-reflected plasma. Well-formed shocks exhibited a soliton-like structure where most of dissipation occurs through maser synchrotron instability. Strong electromagnetic precursor emission was observed in these shocks but no particle acceleration. More systematic study was performed by Gallant et al. (1992) where electromagnetic precursor energy was systematically derived. The most efficient precursor emission was observed for $\sigma \sim 0.1$ where it carries about 10% of the incoming kinetic energy. Based on these results, the authors concluded that
magnetized perpendicular relativistic shocks in pair plasma are not efficient particle
accelerators. An interesting particle acceleration mechanism was observed by Hoshino
et al. (1992) in the case where a small fraction of the plasma are ions (electron-positron-
ion composition). In this case the shock structure is modified by ions, even if their
number fraction is small compared to positrons. Magnetosonic waves emitted by ions
are resonantly absorbed by upstream positrons that produces a non-thermal tail in
positron distribution function. Electrons and ions distributions are still Maxwellian.

Concerning the studies of shocks in electron-ion plasma using 1D simula-
tions, Lyubarsky (2006) and Hoshino (2008) explored mildly magnetized regimes.
Lyubarsky explored the effect of the electromagnetic precursor on the incoming elec-
trons and found that the relativistic oscillation of electrons in the field of the wave
results in temperature equipartition between electrons and ions, once they reach down-
stream medium. The theoretical model was confirmed by a 1D PIC simulation with
the mass fraction m<sub>i</sub>/m<sub>e</sub> = 200, upstream Lorentz factor γ<sub>0</sub> = 50, upstream magne-
tization σ = 0.003 and 0.6. No supra-thermal tail in electron distribution was found.
This finding was contrasted by the study of Hoshino (2008) where it was found that the
precursor wave produces a non-thermal tail in ion and electron distribution functions.
The mechanism for particle acceleration is expected to be of wakefield nature.

**Early multi-dimensional studies** The opening of additional degrees of freedom in directions transverse to the shock propagation is essential for relativistic shock physics. For instance, the dominant instability in low-σ regime, Weibel-filamentation, is artificially suppressed in 1D because it is only triggered by a non-zero k⊥, where k⊥ is the transverse wavenumber to the beam propagation direction. This regime is particularly relevant for astrophysics (for instance in the case of the GRB or AGN studies of the propagation of the forward shock) because the magnetization in the ISM is σISM ~ 5 10⁻¹¹ B²/µG/ncm⁻³, where the magnetic field strength is in units of µGauss and the ambient density in units of cm⁻³. Consequently, the observation of filamentation and concomitant trigger of the Fermi process relies on multi-dimensional configuration of the simulation. This restriction is less severe in non-relativistic case for quasi-parallel shocks, where non-resonant streaming or Bell instability can be triggered in 1D and sustains particle scattering (or mirroring) on both sides of the shock front.

Additional step from 2D to 3D is important to correctly deal with particle scattering
properties as the topology of the turbulent magnetic field is different.

The early multi-dimensional exploration of unmagnetized relativistic shocks was
done with full 3D3V but short simulations. Nishikawa et al. (2003) simulated a rel-
ativistic jet with Γ<sub>jet</sub> = 5 propagating into an unmagnetized electron-ion plasma
at rest. These simulations were done using TRISTAN code (Buneman 1993). The
simulation box was small [15 × 8 × 8 × 15](c/ω<sub>pe</sub>)<sup>3</sup> but still large enough to cap-
ture the electron-scale Weibel-filamentation instability and the simulation time was
T<sub>sim</sub> = 23.4 ω<sub>pe</sub><sup>-1</sup> = 5.23 ω<sub>pi</sub><sup>-1</sup>. Capturing only the initial stage (shock not formed)
the authors still demonstrated the importance of the Weibel-filamentation instability.
Frederiksen et al. (2004) performed 3D simulation of unmagnetized colliding plasma
clouds with density ratio n<sub>inj</sub>/n<sub>0</sub> = 3, relative Lorentz factor Γ<sub>jet</sub> = 3 and ion to