# A reduced model for relativistic electron beam transport in solids and dense plasmas

M Touati, J-L Feugeas, Ph Nicolaï, J Santos, L Gremillet, V Tikhonchuk

► To cite this version:

M Touati, J-L Feugeas, Ph Nicolaï, J Santos, L Gremillet, et al.. A reduced model for relativistic electron beam transport in solids and dense plasmas. New Journal of Physics, 2014, 16 (7), pp.073014. 10.1088/1367-2630/16/7/073014. [hal-04361231](https://hal.science/hal-04361231v1)

HAL Id: hal-04361231

https://hal.science/hal-04361231v1

Submitted on 22 Dec 2023

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A reduced model for relativistic electron beam transport in solids and dense plasmas

M Touati<fnref n="1" />, J-L Feugeas<fnref n="1" />, Ph Nicolaï<fnref n="1" />, J J Santos<fnref n="1" />, L Gremillet<fnref n="2" /> and V T Tikhonchuk<fnref n="1" />

E-mail: touati@celia.u-bordeaux1.fr

Received 19 February 2014, revised 25 April 2014

Accepted for publication 16 May 2014

Published 11 July 2014

New Journal of Physics 16 (2014) 073014

doi:10.1088/1367-2630/16/7/073014

## Abstract

A hybrid reduced model for relativistic electron beam transport based on the angular moments of the relativistic kinetic equation with a special closure is presented. It takes into account collective effects with the self-generated electromagnetic fields as well as collisional effects with the slowing down of the relativistic electrons by plasmons, bound and free electrons and their angular scattering on both ions and electrons. This model allows for fast computations of relativistic electron beam transport while describing their energy distribution evolution. Despite the loss of information concerning the angular distribution of the electron beam, the model reproduces analytical estimates in the academic case of a monodirectional and monoenergetic electron beam propagating through a warm and dense plasma and hybrid particle-in-cell simulation results in a realistic laser-generated electron beam transport case.

Keywords: Vlasov-Fokker-Planck, relativistic electron beam, hybrid model

## 1. Introduction

Even if electrons were the first elementary particles discovered in the 19th century thanks to cathode rays [1], their relativistic transport through matter is still an intense field of research

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

citation and DOI.

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

<footnote n="1">¹ University Bordeaux-CNRS-CEA, Centre Lasers Intenses et Applications, UMR 5107,F-33405 Talence, France</footnote>

<footnote n="2">² CEA, DAM, DIF, F-91297 Arpajon, France</footnote>
from the runaway electrons during thunderstorms [2] to the heating by fast electrons in the fast
ignition approach of inertial confinement fusion [3, 4]. A relativistic electron beam of density $n_b$
and of current density $\mathbf{j}_b$ propagating through a plasma of electronic density $n_e \gg n_b$ generates
an electric field which tends to eject the plasma electrons out of the beam volume in order to
equilibrate the total charge in a timescale of the order of the inverse of the plasma electron-ion
collision frequency [5]. Over the same timescale, the electric field accelerates a return current
$\mathbf{j}_p \approx -\mathbf{j}_b$ and decelerates the beam electrons in order to cancel the total current density in
agreement with Lenz's law [5]. This return current allows for the propagation of electron
currents in excess of the Alfven-Lawson limit which defines the maximum relativistic electron
beam current in vacuum [6]. Due to the imperfect current neutralization, a magnetic field is
induced that can deflect the beam electrons. The plasma electron temperature-density crossed
gradients, plasma resistivity gradients and mostly the beam current density curl are the main
sources of this residual magnetic field. Later, this magnetic field begins to diffuse leading to a
separation of the beam and the return current [7, 8]. The collisions of the relativistic electrons
with plasma electrons and ions also contribute to the scattering and the slowing down of the
beam.

The equation which takes into account both collisional and collective transport processes is
the Vlasov-Fokker-Planck (VFP) equation for the beam distribution function $f_b$ [9]. This
equation has been extensively studied for 30 years and several numerical methods have been
developed [10]. Several families of codes can be identified. The first family consists in solving
the VFP equation by a particle-in-cell (PIC) method [11]. Historically, this method was used to
solve physical problems where collisional processes can be neglected. The distribution function
of all electrons (both plasma and beam electrons) is sampled by macroparticles which
consequently leads to accurate but time-consuming computations. Moreover, in order to limit
the non-conservative force associated with the particle-grid mapping which leads to self-heating
and numerical instabilities, the space resolution has to be comparable to the plasma Debye
length. This poses a strong constraint in the case of dense and/or cold plasmas. Collisional
processes are treated by Monte-Carlo methods. The second family of codes employs the same
method but restricted to the beam electrons introducing a low-energy cut off [16]. The plasma
electrons are taken into account via hydrodynamic equations of conservation or simplified ones.
This is the family of hybrid PIC codes. Other codes solve the full VFP equation [12, 13] or use a
decomposition of the distribution function in the momentum space with their corresponding
hybrid versions. It has been shown that a spherical harmonic decomposition of the distribution
function enables modeling of arbitrary local anisotropy for large enough expansion orders [14].
Besides their accuracy and the rapid progress in high performance computing ressources, all
these codes are time consuming because of the Courant-Friedrichs-Lewy condition that
restricts the time step computation combined with the high resolution needed and the large
number of the distribution function variables.

In this context, we have developed a new hybrid model for relativistic electron beam
transport which allows fast computations. The main idea consists in solving the two first angular
moments of the VFP equation complemented with a special closure relation based on the
principle of the Minerbo maximum angular entropy approximation from the radiative transfer
theory [15]. This allows us to close the set of equations by evaluating the 2nd order angular
moment of the distribution function needed in the 1st order angular moment equation. Contrary
to the widely used approximation of the distribution function isotropy with one Legendre
polynomial often called P1, our model M1 allows one to describe the distribution function with an arbitrary local anisotropy.

In the second section, the relativistic VFP is introduced and applied to relativistic electron beam transport through dense matter; the collective effects as well as the collisional effects and the temperature dependence of relativistic electron beam transport are presented, bringing out the expressions of the collisional and collective parameters used in the model. In the third section, the spherical harmonic and the Cartesian tensor scalar product expansions are introduced and the M1 approximation is then presented and discussed. Numerical schemes used to solve the model are then briefly presented. In the fourth section, an academic case of monoenergetic beam propagation through a dense hydrogen plasma is considered in order to demonstrate the major features of the M1 approximation. Analytical expressions are also derived to check the simulation results. In the fifth section, a more realistic simulation of a laser-generated electron beam transport through a thin aluminum target is presented and compared to a simulation made with a hybrid PIC code [16, 17].

## 2. Kinetic description of relativistic electron beam transport in plasma

### 2.1. The relativistic VFP equation for relativistic electron beam transport

We note $f_b(\mathbf{r}, \mathbf{p}, t)$ the distribution function of the relativistic electron beam in the laboratory frame. It is the solution of the VFP equation

$$ \frac{\partial f_b}{\partial t} + \frac{\partial}{\partial \mathbf{r}} \cdot \left\{ \begin{matrix} \frac{\mathbf{p}}{\gamma m_e} f \\ \mathbf{r} \end{matrix} \right\} + \frac{\partial}{\partial \mathbf{p}} \cdot \left[ \begin{matrix} \mathbf{E} + \frac{\mathbf{p}}{\gamma m_e c} \times \mathbf{B} \\ f \end{matrix} \right] = \sum_{\alpha} C \begin{bmatrix} f \\ \alpha \end{bmatrix} \quad (1) $$

where

$$ C \begin{bmatrix} f \\ \alpha \end{bmatrix} = - \frac{\partial}{\partial \mathbf{p}} \cdot \int \mathbf{U} \cdot \left[ f \frac{\partial f_{\alpha}}{\partial f_b} - \frac{\partial f_b}{f} \right] d^3 \mathbf{p} \quad (2) $$

with $\mathbf{U}_\alpha$ the Beliaev-Budker collision operator [18]. The collision operator $C_\alpha[f_b]$ comes from a 2nd order expansion of the Bolztmann collision integral under the small angle scattering assumption

$$ \Delta \mathbf{p} = \Delta p \frac{\mathbf{p}}{\| \mathbf{p} \|} + \Delta p \frac{\mathbf{b}}{\| \mathbf{b} \|} \quad \text{with} \quad p \gg \Delta p \gg \Delta p \cdot \quad (3) $$

$\mathbf{b}$ is the impact parameter vector and $\Delta \mathbf{p}$ is the momentum variation of one relativistic beam electron with momentum $\mathbf{p}$ after colliding with one particle $\alpha$ of momentum $\mathbf{p}_\alpha$. $\alpha$ particles can be ions, bound electrons, free electrons or screened free electrons of the medium where the relativistic electron beam is propagating through. We neglect the binary collisions of the beam electrons with themselves. (3) means that each binary collision occurs in the 2D plane $(\mathbf{p}, \mathbf{p}_\alpha)$ and leads to an exchange of momentum mostly in the perpendicular direction to $\mathbf{p}$. Indeed, even if large angle collisions lead to a large change in momentum after each collision, the probability that they occur is small compared to the probability of small angle collisions. The Beliaev-Budker collision tensor reads
$$ \mathbf{U}_{\alpha} = \frac{1}{2} U_{\alpha,0} \left\{ \mathbf{I} - \frac{1}{(\tilde{\gamma}_{\alpha}^2 - 1)} \left[ \begin{array}{ccc} \mathbf{p} & \mathbf{p} & \mathbf{p} \\ & \ddots & \ddots \\ m_e c & m_e c & m_e c \end{array} \right]^{-\frac{\alpha}{2}} \left( \begin{array}{ccc} \mathbf{p} & -\mathbf{p}_{\alpha} & \mathbf{p} \\ & m_e c & m_e c \end{array} \right) + \tilde{\gamma}_{\alpha} \left( \begin{array}{ccc} \mathbf{p} & -\mathbf{p}_{\alpha} & \mathbf{p} \\ & m_e c & m_e c \end{array} \right) \right\} \quad (4) $$

It is often expressed with the scattering potential $U_{\alpha,0} = 4\pi Z_{\alpha}^2 e^4 \ln \Lambda_{\alpha}/v_{\alpha}$ expanded within the classical Rutherford cross section $(d\sigma_{\alpha}/d\Omega)_{Ruth}$ [9, 18, 19]. But, in a more general case i.e.

without integrating it within a given cross section, the scattering potential reads

$$ U_{\alpha,0} = \frac{1}{2} \int \Delta \mathbf{p}^2 v_{\alpha} d\sigma_{\alpha} \quad (5) $$

$v_{r,\alpha} = c\sqrt{\tilde{\gamma}_{\alpha}^2 - 1}/\tilde{\gamma}_{\alpha}$ is the relative velocity of one particle in the rest frame of the other during their collision, $\tilde{\gamma}_{\alpha} = \gamma_{\alpha}\gamma_{\alpha}/\tilde{\gamma}_{\alpha}^2$ and $\gamma_{\alpha} = 1/\sqrt{(1 - \mathbf{v}_{\alpha} \cdot \mathbf{v}_{\alpha}/c^2)}$. In addition to the small angle scattering assumption (3), we make the assumption that the target particles $\alpha$ remain non-relativistic after each binary collision with the relativistic electrons. That is to say we neglect high energy secondary electrons

$$ |\mathbf{p}_{\alpha}| \ll |\Delta \mathbf{p}| \ll m_{\alpha} c \quad (6) $$

Under the assumption (3), the 2D binary collision problem consists in solving six unknowns (the momenta and the energies of each particle after the collision) while having seven equations (one energy conservation equation, four momenta conservation equations and the two Einstein relationships between energies and momenta). Consequently, there is one relation between two chosen free parameters which are for example the scattering angle $\theta$ and the relative energy loss of the relativistic electron $w = \Delta\gamma/(\gamma - 1)$ in the laboratory frame. So, one is free to work with $d\sigma_d/dw$ instead of $d\sigma_d/d\Omega$. Moreover, under the assumption (6), the energy conservation equation for one collision, $\Delta\gamma_{\alpha} m_{\alpha} c^2 + \Delta\gamma m_e c^2 = 0$, provides

$$ \Delta p_{\perp}^2 = -2m_{\alpha}\Delta\gamma m_e c^2. \quad (7) $$

This naturally leads to a relation between the stopping power of relativistic electrons due to binary collisions and the scattering potential $U_{\alpha,0}$. Knowing the differential cross section $d\sigma_{\alpha}/dw$, one can define the loss of electron kinetic energy $\epsilon = (\gamma - 1)m_e c^2$ per unit of path length $ds$ as follows

$$ \left( \frac{d\epsilon}{ds} \right)_{\alpha} = \epsilon n_{\alpha} \int_{w_{\alpha, \min}}^{w_{\alpha, \max}} w \left( \frac{d\sigma_{\alpha}}{dw} \right) dw \quad (8) $$

with $n_{\alpha}$ the density of the $\alpha$ particles. The integration limits $w_{\alpha, \min}$ and $w_{\alpha, \max}$ in (8) depend on the nature of collisions. For collisions on free electrons, $w_{\text{free,max}} = 1/2$ due to the indiscernibility of electrons, which means that this is the electron with the higher energy which is considered as primary and a cut-off $w_{\text{free,min}} = w_c$ is used to distinguish full binary collisions and collisions for which collective effects take place. That is to say for collisions on bound and screened free electrons. These latter collisions can be understood respectively in terms of energy transfer to
the excitation of bound electrons by the beam electron electric field and in terms of energy transfer to plasma waves in quanta of $\omega_{pe}$ (plasmons) where $\omega_{pe} = \sqrt{4\pi n_e e^2/m_e}$ is the Langmuir plasma frequency. For collisions on plasma ions, $w_{i,min}$ and $w_{i,max}$ correspond respectively to the commonly used impact parameters $b_{max} = \max\{\lambda_{\text{Debye}} (3/4\pi n_i)^{1/3}\}$ and $b_{min} = \lambda_{\text{de Broglie}} = \frac{\omega_{pe}}{2p^*}$. * means the value of the beam electron momentum in the collision center of mass frame. $\lambda_{\text{de Broglie}}$ is the de Broglie wavelength of the electron, $\lambda_{\text{Debye}}$ is the screening Debye length and $(3/4\pi n_i)^{1/3}$ is the mean inter ionic distance. More details are given in the subsection 2.3.

Injecting the expression (7) in (5) and noticing that $v_{r,a} \sim v$ under the assumptions (3) and (6), one finds the following relation between the scattering potential $U_{a,0}$ and the electron

$$U_{a,0} = - \frac{m_a v}{n_a} \left( \frac{d\epsilon}{ds} \right)_a \qquad (9)$$

Secondly, still under the assumptions (3) and (6), one has $\gamma_a \sim 1$ and $\tilde{\gamma}_a \sim \gamma$. Consequently, the Beliaev-Budker tensor (4) can be simplified into

$$\mathbf{U}_{\alpha} = \frac{[\mathbf{r}]_{\alpha,0}}{2} \left[ \mathbf{p} \otimes \mathbf{p} \right] \left[ \mathbf{I} - \mathbf{p}^2 \right] \qquad (10)$$

and its divergence into

$$\frac{\partial}{\partial \mathbf{p}_{\alpha}} \cdot \mathbf{U}_{\alpha} = - \frac{U_{\alpha,0}}{m_{\alpha} v} \mathbf{p}_{\alpha}. \qquad (11)$$

Instead of their momentum, it is more convenient to work with the variable kinetic energy of the relativistic electrons. The structure of the resulting collision tensor motivates to use spherical coordinates $(\Omega, \mathbf{e}_{\theta}, \mathbf{e}_{\phi})$ where $\Omega = \mathbf{p}' \cdot \mathbf{p}$. Instead of $f_b$, we make use of the distribution function $\Psi$ depending on the kinetic energy $\epsilon$ and the propagation direction solid angle $\Omega$

$$\Psi(\mathbf{r}, \epsilon, \Omega, t) = (p^2/v) f_b(\mathbf{r}, \mathbf{p}, t). \qquad (12)$$

where $p^2/v$ comes from the Jacobian of the change of variables from $\mathbf{p}$ to $(\epsilon, \Omega)$. By defining

$$S(\mathbf{r}, \epsilon, t) = - \sum_{\alpha} \left( \frac{d\epsilon}{d\mathbf{s}_{\alpha}} \right) \qquad (13)$$

$$\text{the total stopping power and} \qquad v(\mathbf{r}, \epsilon, t) = - \sum_{\alpha} \frac{m_{\alpha} v}{p^2} \left( \frac{d\epsilon}{ds} \right)_\alpha \qquad (14)$$

the total angular isotropization rate, one may integrate the collision integral (2) by parts and express the VFP equation in the laboratory frame (1) as
$$ \frac{\partial \Psi}{\partial t} + \frac{\partial}{\partial \Omega} \cdot (v \Omega \Psi) + \frac{\partial}{\partial \Omega} \left[ v (e \mathbf{E} \cdot \mathbf{\Omega} - S) \Psi \right] + \frac{\partial}{\partial \Omega} \left\{ \begin{array}{c} \left( \mathbf{I} - \mathbf{\Omega} \otimes \mathbf{\Omega} \right) \cdot \left[ \begin{array}{c} e \\ \left( \mathbf{E} + \frac{v \mathbf{\Omega}}{c} \times \mathbf{B} \right) \Psi - \frac{\partial}{\partial \mathbf{\Omega}} \left( \frac{v}{c} \Psi \right) \end{array} \right] \\ \left( \mathbf{I} - \mathbf{\Omega} \otimes \mathbf{\Omega} \right) \cdot \left[ \begin{array}{c} \frac{\partial f}{\partial \mathbf{\Omega}} \\ \left( \mathbf{E} + \frac{v \mathbf{\Omega}}{c} \times \mathbf{B} \right) \Psi - \frac{\partial}{\partial \mathbf{\Omega}} \left( \frac{v}{c} \Psi \right) \end{array} \right] \end{array} \right\} = 0 \quad (15) $$

where

$$ \frac{\partial}{\partial \Omega} \cdot \left[ (\mathbf{I} - \mathbf{\Omega} \otimes \mathbf{\Omega}) \cdot \mathbf{A} \right] = \frac{1}{\sin \theta} \frac{\partial}{\partial \theta} \left( \sin \theta A_{\theta} \right) + \frac{1}{\sin \theta} \frac{\partial A_{\varphi}}{\partial \varphi} $$

and

$$ \frac{\partial}{\partial \Omega} \cdot \left[ (\mathbf{I} - \mathbf{\Omega} \otimes \mathbf{\Omega}) \cdot \frac{\partial f}{\partial \mathbf{\Omega}} \right] = \frac{1}{\sin \theta} \frac{\partial}{\partial \theta} \left( \sin \theta \frac{\partial f}{\partial \theta} \right) + \frac{1}{\sin^2 \theta} \frac{\partial^2 f}{\partial \varphi^2} $$

is the Laplace-Beltrami operator.

## 2.2. Collective effects of relativistic electron beam transport in plasma

We consider that the beam is neutralized in charge. That is to say, $n_e = Z^*n_i - n_b \approx Z^*n_i$ ($n_b \ll n_e$), and we neglect the displacement current in the Maxwell-Ampere equation

$$ \frac{\partial}{\partial \mathbf{r}} \times \mathbf{B} = \frac{4\pi}{c} (\mathbf{j}_e + \mathbf{j}_b) \quad (16) $$

while considering times greater than the plasma collisional relaxation time (see [5]). The plasma dynamic is taken into account by the generalized Ohm equation [20]

$$ \mathbf{E} = \eta \mathbf{j}_e + \frac{1}{n_e} \frac{\partial}{\partial \mathbf{r}} (n_e T_e) \quad (17) $$

Considering the space scale greater than the plasma skin depth $c/\omega_{pe}$, which is typically less than a fraction of microns, the electron inertia has been neglected in (17). It has been assumed an isotropic resistivity tensor (no magnetization effects) and the ideal gas expression for the electron pressure $P_e = n_e T_e$. Also, the electron viscosity, the thermal force, the magnetic force and the friction force due to collisions with beam electrons have been neglected compared to the friction force by colliding with plasma particles. In order to account for the induced electric field, one has to add the Maxwell-Faraday equation

$$ \frac{\partial}{\partial \mathbf{r}} \times \mathbf{E} = - \frac{1}{c} \frac{\partial \mathbf{B}}{\partial t} \quad (18) $$

Thus, the system of equations (16), (17) and (18) provides the self-generated electromagnetic field equations

$$ \mathbf{E} = -\eta \mathbf{j}_e + \frac{\eta c}{4\pi} \frac{\partial}{\partial \mathbf{r}} \times \mathbf{B} + \frac{1}{n_e e} \frac{\partial}{\partial \mathbf{r}} (n_e T_e) \quad (19) $$

and

$$ \frac{1}{c} \frac{\partial \mathbf{B}}{\partial t} + \frac{\partial}{\partial \mathbf{r}} \times \left( \frac{\eta c}{4\pi} \frac{\partial}{\partial \mathbf{r}} \times \mathbf{B} \right) = \eta \frac{\partial}{\partial \mathbf{r}} \times \mathbf{j}_e + \frac{\partial \eta}{\partial \mathbf{r}} \times \mathbf{j}_e + \frac{1}{n_e e} \frac{\partial n_e}{\partial \mathbf{r}} \times \frac{\partial T_e}{\partial \mathbf{r}} \quad (20) $$
This last equation describes the self-generated magnetic field diffusion with three different source terms depending on the beam current density, the electrical resistivity and the plasma temperature-density cross gradients, while the self-generated electric field **E** can be evaluated knowing the self-generated magnetic field **B** following (19). These self-generated electro-magnetic fields play an important role in relativistic electron beam transport. The magnetic field due to the curl of the beam current density tends to pinch the relativistic electron beam, the magnetic field due to the resistivity gradients tends to move the relativistic electrons from low electrical resistivity zones to higher electrical resistivity ones, while the resistive electric field slows down the relativistic electrons [21]. The plasma electron temperature-density crossed gradients in (20) may modify the beam transport in a time scale of a few picoseconds [22] but in shorter time scales, this effect can be neglected compared to the previous ones.

## 2.3. Collisional effects of relativistic electron beam transport in plasma

The analysis presented above in section 2.1 shows that the VFP equation (15) usually derived for free electrons can be generalized to a more general case by substituting the Coulomb scattering potential $U_{\alpha,0}$ by the realistic stopping power according to (9). We use the general expression for the stopping powers in solids and dense plasmas

$$ \left( \frac{d\epsilon}{ds} \right)_\alpha = -4\pi \frac{n}{m_\alpha} \frac{Z^2 e^4}{\nu^2} \ln \Lambda_{\alpha\alpha}^{\text{rel}}. \quad (21) $$

$Z_\alpha$ equals 1 for electrons and equals the nuclear charge Z for ions. The Coulomb logarithm $\ln \Lambda_{\alpha\alpha}^{\text{rel}}$ is calculated using the Möller cross section [23] for collisions with free plasma electrons. The cut-off used to distinguish the binary part (collision with plasma free electrons) and the collective one (collisions with plasma bound and screened electrons) is evaluated by $w_c = \Lambda_{\text{de Broglie}}^2 / \max\{\Lambda_{\text{Debye}}^2 (3/4\pi n_i)^{2/3}\}$. The Mott cross section [24] is used for collisions with ions. An extension of the Bethe formula [25] is used with a mean excitation potential $I_{\text{ex}}$ provided by [26] to take into account collisions with plasma bound electrons. The Fermi density effect correction [27] is taken into account according to the Pines and Bohm cross section [28] for collisions with screened free electrons (plasmons). These expressions for stopping powers are derived in the first Born approximation for low and intermediate Z plasmas and only for electrons with kinetic energies greater than approximately 10 keV. The expressions for the densities $n_\alpha$ and the Coulomb logarithms $\ln \Lambda_{\alpha\alpha}^{\text{rel}}$ of (21) are summarized in table 1. Even if bremsstrahlung losses of the relativistic electrons can be neglected in the considered range of energies (10 keV - 10 MeV), a radiative stopping power from [29]

$$ \left( \frac{d\epsilon}{ds} \right)_{\text{rad}} = -4\pi \frac{(Z - Z^*)(Z - Z^* + 1) n_i e^4}{m_e c^2 / \gamma} \left[ \ln(2\gamma) - \frac{1}{3} \right] $$

is added into S. Due to the very low mass ratio $m_e/m_i \ll 1$, the contribution of the stopping power on ions $(d\epsilon/ds)_i$ is negligible compared to those on electrons.

Concerning the angular diffusion of the beam electrons, it is worth noting that the isotropization rate
<figure><img src="image_3.png" /><figcaption>Figure 1. Total stopping power S (collisional and bremsstrahlung losses) (a) and total isotropization rate (b) of relativistic electrons in solid-density and room temperature aluminum (black) and copper (red) and in highly compressed ($\rho = 50$ g cm$^{-3}$) and warm ($T = 1$ eV) hydrogen (blue) as a function of the incident beam electron kinetic energy.</figcaption></figure>

<table><caption>Table 1. Expressions for the Coulomb logarithms $\ln \Lambda_{e\alpha}^{\text{rel}}$ and the densities $n_{\alpha}$ of the stopping powers (21). $Z^*$ is the ionization state and $\ln \Lambda_{e\alpha}^{\text{cfas}} = \max\{2, \ln(b_{\text{max}}/b_{\text{min},\alpha})\}$ is the ‘classical’ Coulomb logarithm where $b_{\text{max}} = \max\{\Lambda_{\text{Debye}}, (3/4\pi n_i)^{1/3}\}$ is the upper impact parameter and $b_{\text{min},\alpha} = \min\{(2m_e c\sqrt{\gamma^2-1}), (2m_e c\sqrt{(\gamma-1)/2})\}$ if $\alpha = \text{ions}$ ($b_{\text{min},\alpha} = \min\{(2m_e c\sqrt{\gamma^2-1}), (2m_e c\sqrt{(\gamma-1)/2})\}$ if $\alpha = \text{free electrons}$) is the lower impact parameter.</caption><thead><tr><th></th><th>n<sub>α</sub></th><th>ln Λ<sup>rel</sup><sub>eα</sub></th></tr></thead><tbody><tr><td>Ions<br>Free e<sup>-</sup></td><td>n<sub>i</sub><br>Z<sup>*</sup>n</td><td>ln Λ<sup>clas</sup><sub>e<sub>i</sub></sub> - (2γ<sup>2</sup> - 1) / (2γ<sup>2</sup>)</td></tr><tr><td>Bound e<sup>-</sup><br>Plasmons</td><td>(Z *Z*)n</td><td>ln (e<sub>e</sub> - ln 2 + 2) + (γ - 1) / (2γ<sup>2</sup>) - (1 / (2γ<sup>2</sup>)) - ln 2 + ln(1 + γ / (16 / (γ - 1)))</td></tr><tr><td></td><td>(Z *Z*)n</td><td>ln | (γ - 1) / (I<sub>ex</sub>) * (γ + 1) / 2 | + (2 / (2γ<sup>2</sup>)) - (1 / (2γ<sup>2</sup>)) - (1 / (γ - 1)) - (1 / (γ - 1))<sup>2</sup></td></tr></tbody></table>

$$v = \sum_{\alpha} \frac{4\pi}{\gamma^2 m_e^2 v^3} \frac{Z^2 e^4}{\ln \Lambda_{e\alpha}^{\text{rel}}} \quad (22)$$

deduced from (14) and (21) does not depend on the targeted $\alpha$ particle mass. Moreover, for intermediate and high Z plasmas, the electron beam scattering on ions is dominant compared to their scattering on electrons by a factor Z. For the hydrogen plasmas, the scattering on both plasma ions and electrons provides comparable contributions. The resulting collisional terms S and v are plotted in figure 1 for solid-density and room temperature copper, aluminum and for highly compressed and warm hydrogen.
## 2.4. Temperature dependence of relativistic electron beam transport

The energy deposition by the electron beam produces a heating and a hydrodynamic motion of the plasma. On the picosecond time scale, the dominant effect is the electron heating while the ion motion is much less important. So, in our simplified model, the ions are assumed to be immobile. Both collisional and collective effects contribute to the plasma heating. The electron temperature $T_e$ of the plasma is calculated in our model according to the following heat equation [20]

$$C_v \frac{\partial T_e}{\partial t} - \frac{\partial}{\partial \mathbf{r}} \cdot \left( \kappa_e \frac{\partial T_e}{\partial \mathbf{r}} \right) = W_e - \frac{C_v T_e - T_i}{Z^* \tau_{eq}} \quad (23)$$

where the advective term, the hydrodynamic velocity divergence, the electron viscosity effects and the thermal force have been neglected. Also, the thermal electron conductivity tensor has been assumed to be scalar (no magnetization effects) and it has been noted $C_v = (3/2) Z^* n_i$ the heat capacity.

$$W_e = \int d\epsilon \int_{S_2} d^2\Omega (vS\Psi) + \eta \mathbf{j}^2 \quad (24)$$

is the heating source term evaluated by calculating the direct collisional energy loss of the relativistic electron beam from (15) added to the Ohmic heating. Here, $S = S_e = (d\epsilon/ds)_{free e^-} + (d\epsilon/ds)_{bound e^-} + (d\epsilon/ds)_p$ due to the fact that $m_e \ll m_i$ which implies $(d\epsilon/ds)_i \ll S_e$. The ion temperature is evaluated from the simple equation

$$\frac{\partial T_i}{\partial t} = \frac{T_e - T_i}{\tau_{eq}} \quad (25)$$

still because the direct frictional energy losses of the relativistic electrons on ions are negligible ($m_i \gg m_e$). Also, the thermal ion conductivity is neglected. The part of the energy converted into density of electromagnetic power

$$W_{em} = \frac{\partial}{\partial t} \left( \frac{\mathbf{E}^2 + \mathbf{B}^2}{8\pi} \right) + c \frac{\partial}{\partial \mathbf{r}} \cdot \mathbf{E} \times \mathbf{B} \quad (26)$$

is negligible compared to $W_e$. The electron-ion equilibration time $\tau_{eq}$ is taken from [30], allowing one to express the thermalization of heated electrons on ions with an arbitrary degree of electron degeneracy. Even if the ideal gas expression that has been used for $C_v = (3/2) Z^* n_i$, that does not take into account the electron degeneracy at low temperatures, the results are not too greatly affected. While the collisional energy losses and the collisional angular diffusion of relativistic electrons weakly (logarithmically) depend on the plasma temperature, the self-generated electromagnetic fields are very sensitive to the different temperature regimes. Indeed, the heat capacity $C_v$, the transport coefficients $\eta$ and $\kappa_e$ as well as the ionization state $Z^*$ depend on the temperature. The latter is evaluated using the expression from [26], except for hydrogen plasmas for which it predicts too high values for $Z^*$ at temperatures below the Fermi temperature $T_F$. In this special case, assuming that the temperature $T_e = T_i = T$ remains lower than $T_F$ during and after the heating by the beam, we take $C_v = (3/2) n_i$ while the electron-ion equilibration term vanishes in (23).
<figure><img src="image_4.png" /><figcaption>Figure 2. Electrical resistivity for aluminum (black) and copper (red) as a function of the electron temperature plotted in two particular cases: $T_i = 300$ K (dashed curves) and $T_i = T_e$ (solid curves) and for hydrogen ($\rho = 50$ g cm$^{-3}$) with $T = T_e$ (solid blue curve).</figcaption></figure>

The electrical resistivity can be written as $\eta = m_e v_e / n_e e^2$ and the thermal conductivity $\kappa_e = \gamma_L T_e / \eta$ where $v_e$ is the electron collision frequency and $\gamma_L$ the Lorenz constant. In the following, the electrical resistivity $\eta$ and the thermal conductivity $\kappa_e$ for the copper and the aluminum are evaluated using the Eidmann-Chimier model [31, 32] which allows one to calculate the transport coefficients with different ion and electron temperatures. In this model, the collision frequency of the electrons $v_e$ is evaluated by taking the average of the plasma electron collision frequency expressions in different temperature regimes $v_e^{-2} = (v_{e-ph} + v_{e-e})^{-2} + v_c^{-2} + v_{\text{Spitzer}}^{-2}$. In the low temperature regime, the mean free path is evaluated by $v_e / (v_{e-ph} + v_{e-e})$ where $v_{e-ph}$ is the collision rate of conduction plasma electrons on phonons, $v_{e-e}$ is the collision rate of conduction plasma electrons on themselves and $v_e = \sqrt{(2T_F + T_e)/m_e}$ is the electron velocity. In the high temperature regime, the mean free path is evaluated by $v_e / N_{\text{Spitzer}}$ where $N_{\text{Spitzer}}$ is the Spitzer collision frequency [34]. In the intermediate range of temperature, the mean free path is written as $v_e / N_c$ where $N_c = v_e / (3/4\pi n_i)^{1/\beta}$. For the hydrogen plasmas, as suggested in [33], the transport coefficients are calculated as $\eta^{-2} = \eta_{\text{Hubbard}}^{-2} + \eta_{\text{Spitzer}}^{-2}$ and $\kappa_e^{-2} = \kappa^2 + \kappa^2$ where $\eta_{\text{Hubbard}}$ and $\kappa_{\text{Hubbard}}$ are the hot temperature Spitzer transport coefficient [34] while $\eta_{\text{Hubbard}}$ and $\kappa_{\text{Hubbard}}$ are the low temperature Hubbard transport coefficients [35]. The transport coefficients $\eta$ used in the model are plotted for aluminum, copper and highly compressed hydrogen in figure 2. Concerning the Hubbard transport coefficients, a least square fit from the ion-ion coupling parameter $\Gamma$ table of [35] is used to avoid discontinuities.
### 3. M1 model for relativistic electron beam transport

#### 3.1. Introduction-spherical harmonic and Cartesian tensor scalar product expansions

A standard method of solution of the VFP equation (15) consists in using a spherical harmonic decomposition. This takes advantage of the fact that the spherical harmonics constitute a full set of orthogonal functions on the unity sphere $S_2$ and are the eigenvectors of the Laplace-Beltrami operator [14].

$$ \Psi = \sum_{l=0}^{\infty} \int_{S_2} \Psi^{\ell} \odot \Omega^{(\ell)} \Omega^{(\ell)} \quad (27) $$

where $\Psi_{\ell}$ is the $\ell$th order angular moment of the distribution function $\Psi$

$$ \Psi_{\ell}(\mathbf{r}, \varepsilon, t) = \int_{S_2} \Psi \Omega^{(\ell)} d^2\Omega, \quad (28) $$

$\Omega^{(\ell)}$ equals the $(\ell - 1)$ tensor products of $\Omega$ with itself $\Omega \otimes \dots \otimes \Omega$ and $\odot_{\ell}$ is the $\ell$ times contracted product. It has been shown that the $N$th order Lagrange polynomial approximation (Cartesian tensor scalar product expansion with $\ell_{max} = N$) is strictly equivalent to the spherical harmonic expansion approximation with $\ell_{max} = N$ [36]. The $l$th order angular moment of the VFP equation can be obtained by integrating (15) multiplied by $\Omega^{(\ell)}$ over the unity sphere $S_2$.

The resulting expression involves the $(\ell + 1)$th order angular moment. Consequently, the $(N + 1)$th order angular moment has to be approximated to close the set of equations. This is done by imposing the PN closure relation

$$ \Psi_{N+1}(\mathbf{r}, \varepsilon, t) = \int_{S_2} \Psi \Omega^{(N+1)} d^2\Omega \approx \int_{S_2} \left| \sum_{l=0}^{N} \frac{1}{4\pi} \ell! \Psi_{\ell} \odot_{\ell} \Omega^{(\ell)} \right| \Omega^{(N+1)} d^2\Omega. $$

In practice, the models are often limited to the 1st order approximation (P1) in order to make fast computations. In this approach, the first two equations read

$$ \frac{\partial \Psi_0}{\partial t} + \frac{\partial}{\partial \mathbf{r}} \cdot (\nu \Psi_1) + \frac{\partial}{\partial \varepsilon} \cdot (\nu (e\Psi_1 \cdot \mathbf{E} - S\Psi_0)) = 0 \quad (29) $$

$$ \frac{\partial \Psi_1}{\partial t} + \frac{\partial}{\partial \mathbf{r}} \cdot (\nu \Psi_2) + \frac{\partial}{\partial \varepsilon} \cdot \left[ \nu (e\Psi_2 \cdot \mathbf{E} - S\Psi_1) \right] = \frac{e}{p} (\Psi_1 - \Psi_2) \cdot \mathbf{E} + \frac{e}{\gamma m_e c} \Psi_1 \times \mathbf{B} - \nu \Psi_2 \cdot \mathbf{I} \quad (30) $$

The second order angular moment $\Psi_2$ is evaluated by using the P1 approximation $\Psi \approx \Psi_{P1} = \Psi_0 \frac{\partial \mathbf{r}}{\partial \ell} + \Psi_1 \cdot \Omega \frac{\partial \mathbf{r}}{\partial \ell}$ of (27)

$$ \Psi_2 = \int \Psi \Omega \otimes \Omega d^2\Omega \approx \int \Psi_{P1} \Omega \otimes \Omega d^2\Omega = \frac{1}{3} \Psi_{P1} \mathbf{I} \quad (31) $$

which is exactly the second order angular moment of an isotropic angular distribution function. Consequently, the P1 approximation (31) is limited to weakly anisotropic distributions and it does not allow one to evaluate the anisotropic part of $\Psi_2$. Moreover, $\Psi_{P1}$ may become negative [39].
## 3.2. The M1 closure

Another approach consists in solving the set of equations (29) and (30) with a closure relation

$$ \Psi \approx \int_{S_2} \Psi \Omega \otimes \Omega d^2\Omega = \frac{1}{\mu} \Psi I + \mu \Psi \left| \begin{array}{c} \frac{\Psi \otimes \Psi}{I} \\ 1 \end{array} - I \right| \quad (32) $$

where $\mu = \mu(\Psi_0, \Psi_1) \in [0, 1]$ is a closure parameter depending on the two first angular moments in such a way that the closure relation (32) is exact for both totally isotropic angular distributions ($\mu = 0$) and totally anisotropic angular distributions ($\mu = 1$). According to the method derived by Minerbo [15], Levermore [37] and Dubroca [38, 39], the approximate distribution function $\Psi_{M1}$ is found by maximizing the local angular entropy of beam electrons with a given kinetic energy under the constraints of the definition of the angular moments $\Psi_0$ and $\Psi_1$. This local angular entropy per unit energy of the beam is defined by

$$ H_{\varepsilon}[\Psi] = - \int_{S_2} \Psi (\ln \Psi - 1) d^2\Omega. \quad (33) $$

The distribution function maximizing $H_{\varepsilon}[\Psi]$ is found by the standard procedure of the Lagrange multipliers $\alpha_0$ and $\alpha_1$ constrained by the definitions (28) of $\Psi_0$ and $\Psi_1$. The Lagrangian of this maximization problem is

$$ L_{\varepsilon}[\Psi] = H_{\varepsilon}[\Psi] - \alpha_0 \left( \Psi_0 - \int_{S_2} \Psi d^2\Omega \right) - \alpha_1 \left( \Psi_1 - \int_{S_2} \Psi \Omega d^2\Omega \right) \quad (34) $$

and the maximizing distribution function $\Psi_{M1}$ is then defined by

$$ \frac{dL_{\varepsilon}}{d\Psi} (\Psi = \Psi) = 0 \quad (35) $$

where $d/d\Psi$ is the functional derivative. The solution of (35) is $\Psi_{M1} = A \exp(\alpha_0 + \alpha_1 \cdot \Omega)$ where the parameter A (which has been introduced for dimensional reasons) and the Lagrange multipliers $\alpha_0$ and $\alpha_1$ have to be evaluated in terms of physical quantities by using the definitions (28) of the moments $\Psi_0$ and $\Psi_1$. One has

$$ \Psi_0 = \int_{S_2} \Psi_{M1} d^2\Omega = A \exp(\alpha_0) \frac{4\pi \sinh |\alpha_1|}{|\alpha_1|} \quad (36) $$

and

$$ \Psi_1 = \int_{S_2} \Psi_{M1} \Omega d\Omega = A \exp(\alpha_0) \frac{4\pi \sinh |\alpha_1|}{|\alpha_1|} \left( \coth |\alpha_1| - \frac{1}{|\alpha_1|} \right) \frac{\alpha_1}{|\alpha_1|} \quad (37) $$

From (36) and (37), one deduces an expression for the anisotropy vector $\Psi_1 = \frac{1}{|\alpha_1|} \alpha_1$

$$ \Omega_{\varepsilon} = \frac{1}{\Psi_0} = \left( \coth |\alpha_1| - \frac{|\alpha_1|}{|\alpha_1|} \right) \quad (38) $$

It is defined as the mean propagation direction of the electrons having the energy $\varepsilon$ at the position $\mathbf{r}$ and time $t$. Due to the triangular inequality applied to (28) with $\ell = 1$, one has necessarily $0 \le |\Omega_{\varepsilon}| \le 1$. Although the bijective relation (38) cannot be inverted analytically,
one can fit the real values of $\alpha_1(\Omega_\epsilon)$ by

$$ \alpha_1(\Omega_\epsilon) \approx \frac{3\Omega_\epsilon}{1 - \frac{|\Omega_\epsilon|^2}{2} \left(1 + |\Omega_\epsilon|\right)^2} \quad (39) $$

Finally, one deduces an explicit expression

$$ \Psi_{M1} = \Psi_0 \frac{|\alpha_1|}{4\pi \sinh |\alpha_1|} \exp (\alpha_1 \cdot \Omega) \quad (40) $$

In the isotropic case where $|\alpha_1| \ll 1$ ($|\Omega_\epsilon| \ll 1$), the M1 model reduces to the one-polynomial approximation P1. But, in the opposite case of a strong anisotropy $|\alpha_1| \to \infty$ ($|\Omega_\epsilon| \to 1$), the function $\Psi_{M1} \to \Psi_0 \delta[\Omega - \Omega_\epsilon]$ where $\delta$ is the Dirac distribution. By substituting (40) in the definition (28) of $\Psi_2$, one obtains

$$ \mu = 1 - 3 \left( \frac{\coth \frac{a_1}{|\Omega_\epsilon|} - \frac{1}{|\Omega_\epsilon|^2}}{a_1} \right) \approx \frac{|\Omega_\epsilon|^2}{|\Omega_\epsilon|^2} \left( 1 + \frac{|\Omega_\epsilon|^2}{|\Omega_\epsilon|^2} \right) \quad (41) $$

The closure (32) is also known as the Eddington approximation which has been developed previously for photons [40]. According to studies on radiative transfer [15, 37, 38, 41], the closure parameter can be written as $\mu = (3\chi - 1)/2$ where $\chi$ is the Eddington factor. This M1 closure relation (32) is a linear combination between the local beam-like case where all electrons at the position $\mathbf{r}$ with the energy $\epsilon$ move in the same direction $\Omega_\epsilon$ and the local isotropic case where all electrons at the position $\mathbf{r}$ with the energy $\epsilon$ can move in all directions with the same probability. It preserves consequently the advantage of the P1 model in fast computing while angular distributions are described with a much better precision. The dependence of the anisotropy vector $\Omega_\epsilon$ on the closure parameter $\mu$ is plotted in figure 3.

It is not evident that the local angular entropy maximization under the constraints of the definition of the angular moments $\Psi_0$ and $\Psi_1$ provides the best possible closure relation. We refer here to Minerbo who justified this closure concerning monoenergetic photons by saying: ‘in communication theory, it is shown that the information content is the negative of the entropy of the distribution. Thus, by using the maximum entropy criterion, one avoids introducing information that is not available. This approach is conceptually superior to the use of an ad hoc model for the intensity.’ Concerning the relativistic electrons considered here, even if their angular scattering tends to isotropize their angular distribution and increase their angular entropy, self-generated electromagnetic fields may not follow the same trend. One can deduce the local angular entropy dissipation rate $(\partial H_\epsilon / \partial t)_{\text{col}}$ starting from the VFP equation (15) and by integrating it over the unity sphere in the momentum space $S_2$. By noting $\phi = \Psi_1 \ln \Psi_1 - \Psi_2$, the equation reads

$$ \frac{\partial H_{\epsilon}}{\partial t} + \frac{\partial}{\partial \mathbf{r}} \cdot \int_{S_2} v \Omega_{\epsilon} \phi d\Omega_{\epsilon} + \left| \int_{S_2} v (eE \cdot \Omega_{\epsilon} - S) \phi d^2\Omega_{\epsilon} \right| = \left( \frac{\partial H_{\epsilon}}{\partial t} \right)_{\text{col}} \quad (42) $$
<figure><img src="image_5.png" /><figcaption>Figure 3. Closure parameter as a function of the anisotropic parameter $|\Omega_{\epsilon}| = |\Psi_1|/\Psi_0$ (solid blue curve) plotted within the approximation (39). The dots are the exact values of $\mu$ for some values of $|\Omega_{\epsilon}|$.</figcaption></figure>

with

$$ \left( \frac{\partial H}{\partial t} \right)_{\text{col}} = \frac{\nu}{2} \int_{S_2} \frac{1}{\phi} \left[ \left( \frac{\partial \phi}{\partial \theta} \right)^2 + \left( \frac{\partial \phi}{\partial \varphi} \right)^2 \right] d^2 \Omega - \Psi \frac{\nu}{p^2} \frac{\partial}{\partial \varepsilon} (p^2 S). \quad (43) $$

The local angular entropy dissipation rate $(\partial H_{\epsilon}/\partial t)_{\text{col}}$ contains two terms. The term depending on the isotropization rate $\nu$ is positive and consequently increases the angular entropy $H_{\epsilon}$. But, the sign of $(\partial H_{\epsilon}/\partial t)_{\text{col}}$ depending also on the total stopping power $S$ is not defined. It means that energy exchange between energy groups of beam electrons does not allow us to consider each group as a closed system. Moreover, the third term in the left hand side of (42) shows that the self-generated electric field affects the angular entropy time evolution. Nevertheless, our experience of using the M1 model confirms that the closure relation (32) allows for a sufficiently accurate and fast computation of the VFP equation (15) with an arbitrary degree of anisotropy. Moreover, it continuously relates the anisotropic and isotropic regimes while satisfying the physical constraints $\Psi_0 \ge 0$ and $0 \le |\Omega_{\epsilon}| \le 1$ thanks to the exponential form (40) [42]. The numerical experiments carried out for the comparison of M1 computations with the full VFP code [13] have shown good agreements with a much reduced computation time [43].

### 3.3. Numerical schemes

The set of equations (29), (30) and (32) is resolved explicitly using a second order HLL Riemann solver [44] (from the name of its founders Harteen, Lax and Van Leer) for the advective terms in space and energy. The collisional slowing down term is solved using the first order downwind scheme. The electron beam density
$$n_b = \int_{\varepsilon_{\min}}^{\infty} \Psi_0 d\varepsilon \qquad (44)$$

and the electron beam current density

$$\mathbf{j}_b = e n_b \mathbf{v}_b = e \int_{\varepsilon_{\min}}^{\infty} \mathbf{\Psi}_{1\nu} d\varepsilon. \qquad (45)$$

are calculated with the two first angular moments $\Psi_0$ and $\Psi_1$, respectively. The self-generated electromagnetic fields (19), (20) are solved explicitly with second order schemes except for the magnetic diffusion term, which is solved implicitly with a second order scheme using a conjugated gradients algorithm. The plasma heat equations (23), (25) are solved explicitly with second order schemes. The energy deposition into the plasma electrons by direct collisional heating by the beam electrons $W_{e,col}$ and the indirect heating $W_{e,res}$ due to the return current are taken into account by writing (24)

$$W_e = W_{e,col} + W_{e,res} \qquad (46)$$

with $W_{e,col} = \int_{\varepsilon_{\min}}^{\infty} S\nu\Psi_0 d\varepsilon$ and $W_{e,res} = \eta \left[ \frac{c \partial}{4\pi \partial r} \times \mathbf{B} - e \int_{\varepsilon_{\min}}^{\infty} \mathbf{\Psi}_{1\nu} d\varepsilon \right]^2$

## 4. 2D-3 V academic case

### 4.1. Introduction

As a first illustration of the M1 model, we consider the simple case of a quasi-monoenergetic and monodirectional anisotropic ($\Omega_{\varepsilon} (z = 0) = e_z$) relativistic electron beam injected at $z = 0$ in a 2D box ($100 \ \mu\text{m} \times 100 \ \mu\text{m}$) of a dense hydrogen plasma with a density $\rho = 50 \ \text{g cm}^{-3}$ and an initial temperature $T_0 = 1 \text{ eV}$. A Gaussian distribution centered at $\varepsilon_0 = (\gamma_0 - 1) m_e c^2 = 1 \text{ MeV}$ with a 50 keV standard deviation is used for the beam energy spectrum. A Gaussian temporal shape centered at $t_0 = 1750 \text{ fs}$ with a standard deviation of $\sigma_t = 500 \text{ fs}$ and a Gaussian spatial shape with a standard deviation of $\sigma_x = 10 \ \mu\text{m}$ have also been used. The electron beam has a total energy of $U = 10 \text{ J}$. As already mentioned in section 2.4, we neglect the electron-ion energy exchange, we assume $T_i = T_e = T$ and $C_v = (3/2) n_i$ in (23). The spatial resolution has been chosen $\Delta x = \Delta z = 1 \ \mu\text{m}$ while the energy resolution has been chosen $\Delta\varepsilon = 5 \text{ keV}$ from $\varepsilon_{\min} = 20 \text{ keV}$ to 1.2 MeV, so that the computation time needed is about 20 hours on 8 CPU. In this academic case, due to the high hydrogen density—which induces a lot of collisions of the beam electrons with plasma particles, and a low plasma electrical resistivity because of plasma electrons degeneracy—for the value of the initialized beam current density

$$j_b(x, z, t) = j_{b0} \exp \left\{ -\frac{x^2}{2\sigma_x^2} - \frac{\left[ z - v(t-t) \right]^2}{2(v_0 \sigma_t)^2} \right\} \qquad (47)$$

with $j_{b0} = Ue/(2\pi)^{3/2} \varepsilon_0 \sigma_x^2 \sigma_t = -1.27 \ 10^{12} \text{ A cm}^{-2}$ and $v_0 = c(1 - 1/\gamma_0^2)^{1/2}$, the collisional effects are predominant compared to the collective ones.
<figure><img src="image_6.png" /><figcaption>Figure 4. Hydrogen temperature $T$ [eV] at $t = 3.5$ ps from the simulation (a) and a comparison of the temperature $T$ profiles at $z = 0.5$ µm from the simulation (solid black) and the estimate (48) (dashed red) at $t = 3.5$ ps (b).</figcaption></figure>

## 4.2. Plasma heating and self-generated electromagnetic fields

By assuming a Dirac distribution in energy centered at $\epsilon_0$ for $\Psi_0$ and by neglecting the Ohmic heating as well as the electron thermal conduction, one can evaluate from the heat equation (23) the hydrogen plasma temperature due to the electron beam energy deposition close to $z = 0$ (to ensure the rigid beam approximation):

$$ T(x, z, t) \approx T_0 + T_1 \exp\left(-\frac{x^2}{2\sigma_x^2}\right) F(z, t) \quad (48) $$

where $F(z, t) = 1 - \text{erf}[(t_0 + z/v_0 - t)/\sigma_\sqrt{2}]$, erf is the error function and $T_1 = S(\epsilon_0)U/4\pi^2 G_\xi \sigma^2 \approx 17.6$ eV. A comparison of the simulation profile at $z = 0.5$ µm and $t = 3.5$ ps with the estimate (48) shows a good agreement (see figure 4(b)). That confirms that neglecting the temperature diffusion as well as the indirect electron beam energy deposition via Ohmic heating $W_{e,\text{res}}$ is a good approximation. Indeed, the diffusion time of the temperature is about $C \sigma_x^2 / \kappa \approx 100 - 1000$ ns, which is large compared to the few ps time interval considered here. Along the z-axis, the temperature rises from $z = 0$ to $z \approx 30$ µm reaching a maximum value of $T = 21$ eV and then it decreases to the initial value $T_0$ due to the beam electron energy losses discussed in the next section.

By neglecting the resistive diffusion of the magnetic field $(\eta c^2 (t - z/v_0 - t_0)/2\pi\sigma_x^2 \ll 1$ at the considered times of a few ps) and by approximating the temperature dependence of the resistivity in the self-generated magnetic diffusion equation (20) as $\eta \approx \eta_0 (T/T_0)^\alpha$ where $\eta_0 = 9 \times 10^{-9} \Omega \text{ m}$ and $\alpha = 1.4$ in agreement with the Hubbard theory in this regime (see figure 2), the estimate (48) of the temperature allows one to evaluate the self-generated magnetic field close to $z = 0$:
<figure><img src="image_7.png" /><figcaption>Figure 5. Magnetic field $B_y$ [T] at 3.5 ps from the simulation (a) and comparison between the magnetic field $B_y = B_{j,y} + B_{\eta,y}$ profile at $z = 0.5$ µm and $t = 3.5$ ps (solid black) where $B_{y,j}$ is the contribution due to the beam current density gradients (solid red) and $B_{j,\eta}$ the contribution due to the resistivity gradient (solid blue) from the simulation and their respective estimates (49) (dashed black), (50) (dashed red) and (51) (dashed blue) (b).</figcaption></figure>

$$B_y = B_{y,j} + B_{y,\eta} \qquad (49)$$

where

$$B_{y,j} \approx B^0 \frac{T}{\alpha + 1} \left[ \left( \frac{T}{T} \right)^{\alpha+1} \right]^{-1} \qquad (50)$$

is the contribution due to the beam current density and

$$B_{y,\eta} \approx B^0 \frac{T}{\alpha + 1} \left[ \frac{\alpha (x/\sigma_x) \exp\left(-\frac{x^2}{2\sigma_x^2}\right)}{1 + \exp\left(-\frac{x^2}{2\sigma_x^2}\right)} - F(z, t) \right] \qquad (51)$$

is the contribution due to the resistivity gradients. Here, $B_0 = j_{b0} (\eta c \sigma_x) \sigma_t \sqrt{\pi/2} \approx -7.95$ T and $F(z, t)$ of (48) has been approximated by the Heavyside function $H(t - z, \zeta_0 - t_0)$. These analytic estimates are plotted and compared with the simulation results in figure 5(b). It confirms that the resistivity gradient makes a significant contribution to the self-generated magnetic field even for hydrogen temperatures below 20 eV. The temperature-density crossed gradients do not contribute to the magnetic field generation due to the fact that the plasma electron density is constant. By neglecting the plasma pressure gradients and the self-generated magnetic field in (19), one can also evaluate the self-generated electric field

$$E_z \approx -\eta j_{b0} \exp \left\{ -\frac{x^2}{2\sigma_x^2} - \frac{[z - v(t - t)]^2}{2(v_0 \sigma_t)^2} \right\}. \qquad (52)$$
<figure><img src="image_8.png" /><figcaption>Figure 6. Electron beam density $n_b[\text{cm}^{-3}]$ at $t = 1750$ fs from the simulation. The electron trajectory (oscillating red curve) is arbitrary.</figcaption></figure>

The maximum value of the slowing down electric force $-e\mathbf{E} \cdot \mathbf{\Omega}_\varepsilon$ is $e\eta j_{b0} \approx 0.1 \text{ keV } \mu\text{m}^{-1}$. It is very small compared to $S$ (see figure 1(a)). This confirms again that the resistive heating $W_{e,\text{res}}$ in the heating equation (23) is negligible. According to the estimates (50) and (51), the maximum of the beam cyclotron frequency $\omega_c = eB_y/\gamma m_e c$ is about $eB_0/m_e c \approx 1 \text{ ps}^{-1}$ at $x = \sigma_x$, $z = 0$ and $t = 3.5$ ps. Consequently, the inequality $\omega_c \ll \nu$ is verified in this particular case (see figure 1(b)). Thus, the effects of the self-generated electromagnetic fields are negligible and the evolution of the electron beam is purely collisional.

### 4.3. Kinetic evolution of the electron beam

One can evaluate the mean position on the z-axis of a beam electron with an initial energy $\varepsilon_0$ and an initial velocity $v_0 = v_0 e_z$ at $z = 0$ by

$$ \langle z \rangle (\boldsymbol{\varepsilon}) = \int_{0}^{s} \langle \cos \theta \rangle (s) \, ds = \int_{\varepsilon_0}^{\varepsilon} \langle \cos \theta \rangle (\varepsilon) \left( \frac{d\boldsymbol{\varepsilon}}{ds} \right)^{-1} d\varepsilon. \quad (53) $$

Notations are illustrated in figure 6. $\langle \cos \theta \rangle$ is the mean cosine of the angle between the z-axis and the position of the beam electron and $s$ is the path length following the electron trajectory in the $(\mathbf{r}, \boldsymbol{v})$ space. In this academic case, the total stopping power can be written

$$ \frac{d\boldsymbol{\varepsilon}}{ds} = -S = -4\pi \frac{n r^2 m c^2}{\varepsilon^2} \ln \Lambda_e \quad (54) $$

where $\ln \Lambda_e$ is the sum of the Coulomb logarithms $\ln \Lambda_{rel}$ of the beam electrons scattering on the bound, free and screened plasma electrons (21) (see table 1), $r_e$ is the classical electron radius and the stopping power on plasma ions has been neglected since $m_e/m_i \ll 1$. For the case of a plasma with generated electrons ($I = 1 - 20 \text{ eV} \ll I_F = 351 \text{ eV}$ for $\rho = 50 \text{ g cm}^{-3}$), the drag number $\ln \Lambda_e$ can be evaluated in the limit $\beta \to 1$ as
$$ \ln \Lambda_e = \ln \frac{m_e c^2}{\Box \omega_{pe}} + \frac{9}{16} - \frac{\ln 2}{2} + f(\gamma) \quad (55) $$

where $f(\gamma) = \ln(\beta\sqrt{\gamma-1}) - [(1/8) + \ln 2]/\gamma + [(1/16) + (1/2)\ln 2]/\gamma^2$ [45, 48].

The mean cosine $\langle \cos \theta \rangle$ can be evaluated in the M1 model (29), (30) and (32) by noticing that $\langle \cos \theta \rangle = \Omega_{\epsilon} \cdot e_z$. Then, neglecting the self-generated electromagnetic field $E$ and $B$, one finds

$$ \frac{d}{ds} \langle \cos \theta \rangle = -k \langle \cos \theta \rangle - \frac{1}{\epsilon} \frac{\partial}{\partial \epsilon} \cdot (\Pi \cdot e) \quad (56) $$

where $\Pi = \frac{1-\mu}{\epsilon} \Psi I + \Psi \mu - \frac{\Omega}{\epsilon} \Psi_1 \otimes \Psi_1 \frac{\partial \mathbf{r}}{|\Psi_1|^2}$ and $\frac{\partial}{\partial \epsilon} = \frac{1}{\epsilon} [\frac{\partial}{\partial \epsilon} + \nu \Omega \cdot \frac{\partial}{\partial \epsilon} - S \nu \frac{\partial}{\partial \epsilon}]$. $\nu = \frac{1}{\epsilon}$, $\frac{\partial}{\partial t} = \frac{\partial}{\partial \tau}$, $\frac{\partial}{\partial \mathbf{r}} = \frac{\partial}{\partial \mathbf{r}}$, $\frac{\partial}{\partial \mathbf{\epsilon}} = \frac{\partial}{\partial \mathbf{\epsilon}}$.

It has also been noted

$$ k = \frac{\nu}{\epsilon} = 4\pi \frac{n r^2}{\nu^2 \beta^4} (\ln \Lambda_e + \ln \Lambda_{ei}^{\text{rel}}) \quad (57) $$

the inverse of the beam electrons mean free path where $\ln \Lambda_{ei}^{\text{rel}} = \ln [2(3 \pi n)^{1/3} / (\Box / m_e)] - 1 + \ln(\sqrt{\gamma^2 - 1}) + 1/\gamma^2$ is the Coulomb logarithm from the stopping power of the beam electrons on ions (21). Assuming that $|\Omega_{\epsilon}| \approx 1$ which induces $\Pi_{\epsilon} \approx 0$, one may neglect the second term in the right hand side of (56) and one obtains, in agreement with the multiple scattering theory of Lewis [46], that

$$ \langle \cos \theta \rangle (\epsilon) \approx \exp \left( - \int_{0}^{s} k_{1}(s) ds \right) = \exp \left( - \int_{\epsilon_0}^{\epsilon} k_{1}(\epsilon) \left( \frac{d\epsilon}{ds} \right)^{-1} d\epsilon \right) \quad (58) $$

As was suggested in [45] and [47], the ratio $\ln \Lambda_{ei}^{\text{rel}} / \ln \Lambda_e$ can be considered as a constant and one obtains

$$ \langle \cos \theta \rangle \approx \left[ \frac{(\gamma - 1)(\gamma + 1)}{(\gamma - 1)(\gamma + 1)} \right]^{\frac{\ln \Lambda_e + \ln \Lambda_{ei}^{\text{rel}}}{2}} \quad (59) $$

Following the arguments of [48], we neglect $f(\gamma)$ in $(d\epsilon/ds)^{-1}$ of (53) compared to $\ln \Lambda_e^* = \ln(m_e c^2 / \Box \omega_{pe}) + (9/16) - (\ln 2/2) \approx 7.98$ for $\rho = 50$ g cm$^{-3}$. Moreover, by considering $(\ln \Lambda_e + \ln \Lambda_{ei}^{\text{rel}}) / 2 \ln \Lambda_e \approx 1$ in $\langle \cos \theta \rangle$ of (53), the mean electron propagation distance can be estimated as

$$ \langle z \rangle \approx 4\pi n_e r_e^2 \ln \Lambda_e^* \frac{\nu_0^0}{\nu_0} - 1 \left[ \frac{\nu_0^0}{\nu_0} - \frac{\nu_0^0}{\nu_0} - 2 \ln \left( \frac{\nu_0^0}{\nu_0} \right) \right] \quad (60) $$

Then, the penetration depth of the beam electrons with an initial kinetic energy $\epsilon_0$ and an initial velocity $v_0 = v_0 e_z$ at $z = 0$ can be written as
$$L_p = \langle z \rangle (\varepsilon \to 0) = \xi R \qquad (61)$$

where

$$R = \int_{\varepsilon_0}^{0} \left( \frac{d\varepsilon}{ds} \right)^{-1} d\varepsilon = \frac{1}{4\pi m_e r_e^2 \ln \Lambda_e^*} \frac{\left( \varepsilon_0 / m_e c^2 \right)^2}{1 + \varepsilon_0 / m_e c^2} \qquad (62)$$

is the range of the beam electrons with an initial kinetic energy $\varepsilon_0$ and

$$\xi = \left( \frac{\gamma + 1}{\gamma_0 - 1} \right)^2 \left( \frac{1}{\beta_0^2} \left( \beta^2 - 2 \frac{\ln \gamma_0}{\gamma_0} \right) \right) \qquad (63)$$

is the correction due to their angular scattering. It is equal to $2/3$ when $\gamma_0 \to 1$ and it increases to 1 when $\gamma_0 \to \infty$. These values are in agreement with the approximation of the penetration depth (30) in [48], which assumes a factor $\xi \approx 0.8$ to recover the Monte-Carlo simulation. For $\rho = 50$ g cm$^{-3}$ and $\varepsilon_0 = 1$ MeV, one finds $R \approx 54$ µm and $\xi \approx 0.7$, which corresponds to $L_p \approx 38$ µm. This is in agreement with our simulation results shown in figure 4(a) and figure 6. The estimates (59) and (60) predict the mean position $\langle z \rangle$ and the mean diffusion angle $\arccos \langle \cos \theta \rangle$ of $\approx 10$ µm and $\approx 20^\circ$ respectively for $p_z/m_ec = 2.5$, $\approx 20$ µm and $\approx 30^\circ$ for $p_z/m_ec = 2$ and $\approx 30$ µm and $\approx 45^\circ$ for $p_z/m_ec = 1.7$. This is in agreement with the numerical results obtained for the electron beam distribution function in the M1 approximation (40) as illustrated in figures 7(b), (c) and (d), respectively. Close to $z = L_p$, the analytic estimates $\langle z \rangle$ and $\langle \cos \theta \rangle$ differ from the numerical results as shown in figure 8. This is due to the singularity at the penetration depth $L_p$ contained in (59) and (60) and to the fact that the last term in the right hand side of (56) cannot be neglected anymore when $\gamma \to 1$.

## 5. Realistic relativistic electron beam transport simulation

The M1 model is also compared to hybrid PIC simulations of a relativistic electron beam propagation in a thin solid target, motivated by an experimental campaign [49]. The target is composed of three successive layers of 1 µm of aluminum, 3 µm of copper and 1 µm of aluminum. A linearly polarized laser pulse with a wavelength $\lambda = 800$ nm, a total energy $E_L = 0.7$ J and a 26 fs full width half maximum time duration is focused with a peak intensity of $I_L = 3 \times 10^{19}$ W cm$^{-2}$ at a $45^\circ$ incident angle. Plasma mirrors have been used during these experiments to avoid prepulse/preplasma so that the main electron acceleration mechanism during the laser-target interaction is the $j \times B$ heating. Thus, the accelerated electrons propagate mainly in the laser pulse propagation direction. The beam's initial properties are obtained from a 2D fully PIC simulation of the laser-plasma interaction using the code CALDER [50-52]. The beam transport in the dense part of the target is simulated by the hybrid PIC code PARIS [16, 17] and by the M1 code described in the previous sections. A low cut-off energy $\varepsilon \gg \varepsilon_{min} = 20$ keV has been used to distinguish between bulk electrons and beam electrons. The energy distribution of the laser-generated electron beam, its angular and spatial distribution and the instantaneous conversion efficiency from the laser to the beam have been interpolated and adapted to initialize the angular moments $\Psi_0$ and $\Psi_1$ at $z = 0$. For this transport simulation, the
<figure><img src="image_9.png" /><figcaption>Figure 7. Electron beam distribution function $f_b [\text{cm}^{-3} \cdot \text{m}_e \text{c}^{-3}]$ from the simulation on the z-axis at $t = 1750$ fs, $p_x = 0$ and at different depths $z = 1$ µm (a), $z = 10$ µm (b), $z = 20$ µm (c) and $z = 30$ µm (d). The dashed red curves represent the analytical estimates of arccos $\langle \cos \theta \rangle (\varepsilon)$ evaluated at the kinetic energies $\varepsilon$ corresponding to $\langle z \rangle (\varepsilon) = 10$ µm (b), $20$ µm (c) and $30$ µm (d).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 8. Comparison between the simulation results at $t = t$: $\arccos(\Omega_{\varepsilon,z}^{\text{max}})$ versus $z$ where $\Omega_{\varepsilon,z}^{\text{max}} = \Omega_{\varepsilon,z}(x=0, z, \varepsilon, t) = \max_{m,0} \{\Omega_{\varepsilon,z}(x=0, z, \varepsilon, t)\}$ (black curve), $\arccos(\Omega_{\varepsilon,z}^{\text{mean}})$ versus $z$ where $\Omega_{\varepsilon,z}^{\text{mean}} = \langle \Omega_{\varepsilon,z}(x=0, z, \varepsilon, t)_0 \rangle$ over $|\varepsilon - \varepsilon_m| < 200$ keV (blue curve) and the analytical estimates $\arccos(\cos \theta)(\varepsilon)$ (from (59)) versus $\langle z \rangle(\varepsilon)$ (from (60)) (red curve).</figcaption></figure>

spatial resolution has been chosen $\Delta x = \Delta z = 0.25$ µm while the energy resolution has been chosen $\Delta \varepsilon = 10$ keV, in the range from $\varepsilon_{\min} = 20$ keV to 3 MeV. The computation time needed is about 20 hours on 8 CPU. Absorbing conditions at the target borders have been used so that the refluxing of the beam electrons at both the rear and the irradiated side of the target was suppressed.

Almost all beam electrons reach the rear side of the target after approximately 100 fs (see figure 9). The maximum beam density is close to the critical density $n_{b,max} \approx 10^{21} \text{ cm}^{-3}$ while the maximum value of the beam current density $j_{b,max}$ is above $10^{12} \text{ A.cm}^{-2}$. Thus, the collective effects play an important role in the beam transport [53] (aluminum and copper electrical resistivity is much higher than that of hydrogen studied in the previous section; see figure 2). The self-generated magnetic field reaches its maximum value of approximately 200 T in the first aluminum layer at the end of the beam propagation and then decreases away down to 100 T at 500 fs (see figure 10(c)). Indeed, the main contribution to the magnetic fields is due to the curl of the beam current density but the electrical resistivity gradients also play an important role: at 27 fs, the plasma electrons have been heated up to the Fermi temperature in the first aluminum layer. Consequently, the electrical resistivity in this area goes from the cold solid-liquid phase to the hot plasma one and decreases with the temperature. The same scenario appears in the copper layer at 40 fs and later in the rear aluminum layer. The consequence is that it creates electrical resistivity gradients in the target which tend to hollow the beam (see figures 10(a) and (b)), as explained in [54].

The time evolution of the energy in the target is plotted in figure 11. The energy is conserved within a 0.5 % error in this simulation. The total injected energy at $z=0$ is $\approx 50$ mJ, which represents a conversion efficiency from the laser to the electron beam of $\approx 7\%$. The electromagnetic energy is negligible compared to the beam energy by a factor $\approx 1000$. The
<figure><img src="image_11.png" /><figcaption>Figure 9. Electron beam density $n_b$ [cm⁻³] from the M1 simulation at $t = 25$ fs (a), $t = 50$ fs (b), $t = 75$ fs (c), and $t = 100$ fs (d).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 10. Plasma electron temperature $T_e$ [eV] with a logarithm scale (a), plasma electrical resistivity $\eta$ [Ω m] with a logarithm scale (b), self-generated magnetic field $B_y$ [T] (c) and time integrated density of Kα photons emitted per steradian $n_{K\alpha}$ [cm⁻³ sr⁻¹] (d) from the M1 simulation at $t = 500$ fs.</figcaption></figure>

heating of the target due to the return current (≈ 10 mJ) exceeds by roughly two times the direct collisional heating by the beam electrons.

In spite of the relative complexity of this laser-generated electron beam transport, the results of this M1 simulation are close to those of the hybrid PIC simulation. The same behavior
<figure><img src="image_13.png" /><figcaption>Figure 11. Instantaneous beam energy (solid black), integrated beam energy balanced between injected and escaping electrons at $z = 0$ (solid red), escaped energy at $z = 5$ µm (dashed red), total collisional energy loss (solid magenta), total ‘collective’ (resistive) energy loss (solid blue), instantaneous electric energy ×1000 (solid green) and instantaneous magnetic energy ×1000 (solid cyan) from the M1 simulation.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 12. Comparison at $t = 500$ fs between the mean electron temperature $\langle T_e \rangle$ over $|x_{max} - x| < 5$ µm profile obtained with M1 (blue) and that obtained with the hybrid PIC code Paris (green). $x_{max}$ is defined as the position where $T_e(x_{max}, z)$ is the maximum electron temperature at a given depth $z$ (red).</figcaption></figure>

is recovered for the temperature profile of the target as shown in figure 12. One can observe the signature of a predominant Ohmic heating in the copper layer ($z = 1 \rightarrow 4$ µm) temperature profile of the target. Indeed, copper is less resistive than aluminum; if direct collisional losses were predominant, due to the fact that copper is denser than aluminum, one would see the opposite case with a more important heating in the copper layer. The time integrated density of Kα emission of radiations per steradian have been computed within the empirical expression for K-shell ionization cross section by electron impact [55] and the K-shell fluorescence yield probability provided by [56]. It is plotted in figure 10(d) at $t = 500$ fs. There is an important
discrepancy between the simulation result and the experimental data concerning the size of the $K\alpha$ spot size. Indeed, the refluxing of the beam electrons across the target which enhances their lateral expansion, and consequently increases the induced $K\alpha$ emission of radiations spot size, is not taken into account here.

## 6. Conclusion

The relativistic kinetic equation has been adapted for the description of relativistic electron beam transport in solids and dense plasmas by expressing the Beliaev-Budker scattering potential with the stopping powers of electrons. The collective effects as well as the collisional effects have been highlighted to bring out the key parameters of relativistic electron beam transport such as the total stopping power and the isotropization rate of beam electrons, the evolution of the electrical resistivity, the thermal conductivity and the electron-ion equilibration time of the dense plasma. We presented a reduced model suitable to numerically solve the VFP equation with less computation time than existing models. The M1 model consists of approaching the solutions of the relativistic VFP equation with its two first angular moments and by approximating the second order angular moment of the distribution function using the Minerbo maximum entropy closure from the radiative transfer theory. An academic simulation of a monoenergetic, monodirectional and anisotropic MeV electron beam propagating through a highly compressed and cold hydrogen plasma is presented. In this ideal situation, the collisional effects are predominant compared to collective ones. It allows one to check the model by deriving analytical expressions for the temperature and the self-generated magnetic-fields induced by the beam as well as the mean scattering angle and the penetration depth of beam electrons. The comparison between these estimates and the simulation shows good agreement. Despite the local angular average performed to each energy group of beam electrons, the model also reproduces the results obtained with a hybrid PIC simulation concerning a realistic laser-generated electron beam transport. The electron beam is injected in a thin aluminum target with an embedded copper layer which has been used during an experimental campaign conducted on the UHI100 laser. The distribution of the beam is deduced from a fully 2D PIC simulation of the laser-plasma interaction. The obtained electron beam density and its current density as well as the self-generated electromagnetic fields shows the same comportment. The final temperature in the target recovers the hybrid PIC results in this transport case where the collective effects are predominant compared to collisional ones. Besides, in order to reproduce the experimental data, the refluxing of beam electrons at both the irradiated and the rear side of the target have to be taken into account. The future development of the M1 model may include a second set of M1 equations to describe the counter propagating electron beam.

## Acknowledgments

The work presented here was performed under the auspices of the French Agency for Research and the joint competitiveness cluster Alpha Route des lasers through project TERRE ANR-2011-BS04-014, and the EURATOM within the ‘Keep-in-Touch’ activities. The authors thank B Dubroca and his PhD student S Guisset for fruitful discussions concerning the numerical aspects of the paper.
References

[1] Thomson J J 1901 On bodies smaller than atoms Pop. Sci. 59 325-36

[2] Gurevich A V, Mesyats G A, Zybin K P, Yalandin M I, Reutova A G, Shpak V G and Shunailov S A 2012 Observation of the avalanche of runaway electrons in air in a strong electric field *Phys. Rev. Lett.* **109** 085002

[3] Tabak M, Hammer J, Glinsky M E, Kruer W L, Wilks S C, Woodworth J, Campbell E M, Perry M D and Mason R J 1994 Ignition and high gain with ultrapowerful lasers *Phys. Plasmas* **1** 1626

[4] Strozzi D J, Tabak M, Larson D J, Divol L, Kemp A J, Bellei C, Marinak M M and Key M H 2012 Fast-ignition transport studies: realistic electron source, integrated particle-in-cell and hydrodynamic modeling, imposed magnetic fields *Phys. Plasmas* **19** 072711

[5] Hammer D A and Rostoker N 1970 Propagation of high current relativistic electron beams *Phys. Fluids* **13** 1831

[6] Alfvén H 1939 On the motion of cosmic rays in interstellar space *Phys. Rev.* **5** 425

[7] Lee R and Sudan R N 1971 Return current induced by a relativistic beam propagating in a magnetized plasma
*Phys. Fluids* **14** 1213

[8] Bret A, Gremillet L and Dieckmann M E 2010 Multidimensional electron beam-plasma instabilities in the relativistic regime *Phys. Plasmas* **17** 120501

[9] Landau L D and Lifshitz E M 1981 A course of theoretical physics *Phys. Kinetics* **10** 168-216

[10] Thomas A G R, Tzoufras M, Robinson A P L, Kingham R J, Ridgers C P, Sherlock M and Bell A R 2012 A review of Vlasov-Fokker-Planck numerical modeling of inertial confinement fusion plasma *J. Comput. Phys.* **231** 1051

[11] Takizuka T and Abe H 1977 A binary collision model for plasma simulation with a particle code *J. Comput. Phys.* **25** 205

[12] Yokota T, Nakao Y, Johzaki T and Mima K 2006 Two-dimensional relativistic Fokker-Planck model for core plasma heating in fast ignition targets *Phys. Plasmas* **13** 022702

[13] Duclous R, Dubroca B, Filbet F and Tikhonchuk V 2009 High order resolution of the Maxwell-Fokker-Planck-Landau model intended for ICF applications *J. Comput. Phys.* **228** 5072

[14] Tzoufras M, Bell A R, Norreys P A and Tsung F S 2011 A Vlasov-Fokker-Planck code for high energy density physics *J. Comput. Phys.* **230** 6475

[15] Minerbo G N 1977 Maximum entropy Eddington factors LA Report UC-32 LA-6882-MS

[16] Gremillet L, Bonnaud G and Amiranoff F 2002 Filamented transport of laser-generated relativistic electrons penetrating a solid target *Phys. Plasmas* **9** 941

[17] Martinolli E et al 2006 Fast-electron transport and heating of solid targets in high-intensity laser interactions measured by Kα fluorescence *Phys. Rev. E* **73** 046402

[18] Belyaev S T and Budker G I 1956 The relativistic kinetic equation *Sov. Phys. Doklady* **107** 807

[19] Braams B J and Karney C F 1987 Differential form of the collision integral for a relativistic plasma *Phys. Rev. Lett.* **59** 1817

[20] Braginskii S 1965 *Reviews of Plasma Physics* ed M Leontovitch (New York: Consultants Bureau) p 205

[21] Davies J R, Bell A R, Haines M G and Guérin S M 1997 Short-pulse high-intensity laser-generated fast electron transport into thick solid targets *Phys. Rev. E* **56** 7193

[22] Nicolaï Ph, Feugeas J-L, Regan C, Olazabal-Loumé M, Breil J, Dubroca B, Morreeuw J-P and Tikhonchuk V 2011 Effect of the plasma-generated magnetic field on relativistic electron transport *Phys. Rev. E* **84** 016402

[23] Möller C 1932 Zur Theorie des durchgangs schneller electronen durch materie *Ann. Phys.* **14** 531

[24] Mott N F 1932 The polarization of electrons by double scattering *Proc. R. Soc. London A* **135** 429

[25] Bethe H 1932 Bremsformel fr elektronen relativistischer geschwindigkeit Z. Phys. **76** 293

[26] More R M 1985 Pressure ionization, resonances and the continuity of bound and free states *Adv. At. Mol. Phys.* **21** 305
[27] Fermi E 1940 The ionization loss of energy in gases and in condensed materials *Phys. Rev.* **57** 485

[28] Pines D and Bohm D 1952 A collective description of electron interactions: II. Collective vs individual particle aspects of the interactions *Phys. Rev.* **85** 338

[29] Heitler W and Sauter F 1933 Stopping of fast particles with emission of radiation and the birth of positive electrons *Nature* **132** 892

[30] Brysk H, Campbell P M and Hammerling P 1975 Thermal conduction in laser fusion *Plasma Phys.* **17** 473

[31] Eidmann K, Meyer-ter-Vehn J, Schlegel T and Hüller S 2000 Hydrodynamic simulation of subpicosecond laser interaction with solid-density matter *Phys. Rev.* E **62** 1202

[32] Chimier B, Tikhonchuk V T and Hallo L 2007 Heating model for metals irradiated by a subpicosecond laser pulse *Phys. Rev.* B **75** 195124

[33] Lambert F, Recoules V, Decoster A, Clérouin J and Desjarlais M 2011 On the transport coefficients of hydrogen in the inertial confinement fusion regime a) *Phys. Plasmas* **18** 056306

[34] Cohen R S, Spitzer L and McR R P 1950 The electrical conductivity of an ionized gas *Phys. Rev.* **80** 230

[35] Hubbard W B 1966 Studies in stellar evolution. V. Transport coefficients of degenerate stellar matter *Astrophys. J.* **146** 858

[36] Johnston T W 1960 Cartesian tensor scalar product and spherical harmonic expansions in Boltzmann's equation *Phys. Rev.* **120** 2277

[37] Levermore C D 1979 A Chapman-Enskog approach to flux-limited diffusion theory *LLNL Report* UCID-18229

[38] Dubroca B and Feugeas J-L 1999 Etude théorique et numérique d'une hiérarchie de modèles aux moments pour le transfert radiatif *C. R. Sci. Paris* t. **329** 915

[39] Dubroca B, Feugeas J-L and Franck M 2010 Angular moment model for the Fokker-Planck equation *Eur. Phys. J.* D **60** 301

[40] Eddington A 1926 *The Internal Constitution of the Stars* (New York: Dover)

[41] Pomeraning G C 1981 Flux limiters and Eddington factors *LLNL Report* UCRL-86561

[42] Wright T P and Hadley G R 1975 Relativistic distribution functions and applications to electron beams *Phys. Rev.* A **12** 686

[43] Mallet J, Brull S and Dubroca B 2014 An entropic scheme for an angular moment model for the classical Fokker-Planck equation of electrons *Commun. Comput. Phys.* **15** 422-50

[44] Berthon C, Dubois J, Dubroca B, Nguyen-Bui T-H and Turpault R 2010 A free streaming contact preserving scheme for the M1 model *Adv. Appl. Math. Mech.* **3** 259-85

[45] Solodov A A and Betti R 2008 Stopping power and range of energetic electrons in dense plasmas of fastignition fusion targets *Phys. Plasmas* **15** 042707

[46] Lewis H W 1950 Multiple scattering in an infinite medium *Phys. Rev.* **78** 526

[47] Robiche J, Rax J-M, Bonnaud G and Gremillet L 2010 Fast electron energy deposition in a magnetized plasma: kinetic theory and particle-in-cell simulation *Phys. Plasmas* **17** 033106

[48] Atzeni S, Schiavi A and Davies J R 2009 Stopping and scattering of relativistic electron beams in dense plasmas and requirements for fast ignition *Plasma Phys. Control. Fusion* **51** 015016

[49] Santos J J *et al* 2013 Supra-thermal electron beam stopping power and guiding in dense plasmas *J. Plasma Phys.* **79** 429

[50] Lefebvre *et al* 2003 Electron and photon production from relativistic laser-plasma interactions *Nucl. Fusion* **43** 629

[51] Nuter R, Gremillet L, Lefebvre E, Lévy A, Ceccotti T and Martin P 2011 Field ionization model implemented in particle in cell code and applied to laser-accelerated carbon ions *Phys. Plasmas* **18** 033107

[52] Pérez F, Gremillet L, Decoster L, Drouin M and Lefebvre E 2012 Improved modeling of relativistic collisions and collisional ionization in particle-in-cell codes *Phys. Plasmas* **19** 083104

[53] Vauzour B, Santos J J, Debayle A and Hulin S 2012 Relativistic high-current electron-beam stopping-power characterization in solids and plasmas: collisional versus resistive effects *Phys. Rev. Lett.* **109** 255002
[54] Davies J R, Green J S and Norreys P A 2006 Electron beam hollowing in laser-solid interactions *Plasma Phys. Control. Fusion* **48** 1181

[55] Hombourger C 1998 An empirical expression for K-shell ionization cross section by electron impact *J. Phys. B: At. Mol. Opt. Phys* **31** 3693

[56] Bambynek W 1984 *X-Ray and Inner-Shell Processes in Atoms, Molecules and Solids* ed A Meisel and J Finster (Leipzig: Karl-Marx-Universitat)