# A unified description of atomic physics for electron Fokker-Planck calculations

Y. Savoye-Peysson, D. Mazon, J. Bielecki, D. Dworak, K. Król, A. Jardin, M. Scholz, J. Walkowiak, J. Decker

► To cite this version:

Y. Savoye-Peysson, D. Mazon, J. Bielecki, D. Dworak, K. Król, et al.. A unified description of atomic physics for electron Fokker-Planck calculations. Nuclear Fusion, 2023, 63 (12), pp.126041. 10.1088/1741-4326/acffd9.cea-04274020

HAL Id: cea-04274020

https://cea.hal.science/cea-04274020

Submitted on 7 Nov 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
A unified description of the atomic physics for electron
Fokker-Planck calculations

Y. Savoye-Peysson¹, D. Mazon¹, J. Bielecki², D. Dworak², K. Król², A. Jardin², M. Scholz², J. Walkowiak² and J. Decker³

¹CEA, IRFM, F-13108, Saint-Paul-lez-Durance, France

²Institute of Nuclear Physics, Polish Academy of Sciences, PL-31342, Krakow, Poland

³Ecole Polytechnique Fédérale de Lausanne (EPFL), Swiss Plasma Center (SCP), CH-1015 Lausanne, Switzerland

**Abstract.** Most of the realistic kinetic calculations for tokamak plasmas require now to incorporate the effect of partially ionized high-Z elements arising either from uncontrolled influxes of metallic impurities like tungsten in high input power regimes or from mitigation of runaway electrons generated after possible major disruptions by massive gas injection. The usual electron-ion Fokker-Planck collision operator must be therefore modified, since all atoms in the plasma are not fully ionized, as it can be considered for light elements. This represents a challenge, in order to perform fast but also accurate calculations, regardless the types of elements present in the plasma, but also their local levels of ionization, while covering a wide range of electron energies in a consistent way, from few keV to tens of MeV in plasmas whose electron temperature may itself vary from ten eV to several keV. In this context, a unified description of the atomic models is proposed, based on a multi-Yukawa representation of the electrostatic potential calibrated against results obtained by advanced quantum calculations. Besides the possibility to improve the description of inner and outer atomic shells in the determination of the atomic form factor, this model allows to derive analytical formulations for both elastic and inelastic scattering which can be then easily incorporated in kinetic calculations. The impact of the number of exponentials in the description of the atomic potential is discussed, and the comparison with simple or advanced atomic models is also performed.

E-mail: yves.savoye-peysson@cea.fr

PACS numbers: 32.10.-f, 34.50.-s, 52.65.Ff
# 1. Introduction

The use of tungsten (W) as the plasma-facing material in present-day experimental fusion devices such as WEST[1], EAST[2] and the International Thermonuclear Experimental Reactor (ITER) currently being constructed [3], has raised the question of the impact of partially ionized high-Z impurities on the performances of hot plasmas. For example, the ability to drive efficiently the toroidal plasma current by Radio-Frequency (RF) electromagnetic waves may be reduced by an enhanced electron pitch-angle scattering and electron-ion slowing-down due to tungsten, thus limiting the capability of control for improved plasma performances [4]. In standard tokamak fusion plasmas, the usual electron temperature $T_e$ is supposed to be in the range between 1 to 10 keV, such that most of low-Z impurities are fully ionized over a large volume, except possibly in the outermost regions near the separatrix. Conversely, high-Z elements remain partially ionized everywhere, even in the core of the plasma, and the nucleus charge $Z_s$ of the species $s$ may still be partially screened by many bound electrons. For the tungsten element whose atomic number is $Z_s = 74$, the mean screened ion charge is $\bar{Z}_{0,s} \simeq 42$ in a plasma whose electron temperature is $T_e = 3$ keV according to the OPEN-ADAS database, as shown in Fig. 1, such that $\bar{N}_s \simeq 32$ electrons are still bound [5]. Even at $T_e = 10$ keV, as expected in ITER plasmas, $\bar{N}_s$ is still large for the tungsten element, of the order of twenty$^\dagger$.

The role played by the screening of partially ionized high-Z elements has been first considered to describe accurately the dynamics of runaway electrons in very cold post-disruptive plasmas but also to investigate the possibility to mitigate them by massive gas injection of high-Z elements up to argon. It is shown that the dynamics in momentum space of the non-thermal electrons can be notably modified as compared to the traditional picture because of the partial screening, with a significant impact on the critical electric field (Hastie-Connors) beyond which electrons may run away [6–8]. This original work, implemented in the *CODE* code dedicated to runaway electrons physics in almost zero-temperature post-disruptive plasmas [9], has been extended later to standard tokamak regimes in the *LUKE* solver of the 3-D linearized bounce-averaged relativistic electron Fokker-Planck equation [10], allowing to describe the consequences of uncontrolled impurity influxes of high-Z metallic elements like tungsten on RF current drive for example [4]. More recently, kinetic calculations have been carried out, showing that RF current driven by the Lower Hybrid wave is moderately lowered despite a strong thermal collapse ascribed to an uncontrolled accumulation of tungsten in the plasma core of WEST tokamak [11]. In both studies, a standard Yukawa potential (single exponential) was used in the *LUKE* code.

Even if the atomic processes which must be described in *CODE* and *LUKE* kinetics codes are rather similar, some differences specific to hot plasmas must be investigated. Indeed, while the atomic physics of argon and elements with lower $Z_s$ values has been thoroughly studied by quantum non-relativistic codes describing ground-state and mean excitation energies for different ionization states [12–15], the knowledge of atomic properties for metallic elements with higher $Z_s$ values, and in particular for tungsten, is much more sparse. This is the consequence of the relativistic effects and the resulting complex orbitals coupling, which must be fully incorporated in quantum calculations, making them considerably more difficult. Indeed, by combining the virial theorem with the quantum uncertainty principle, relativistic effects become significant when the relativistic Lorentz factor $\gamma_s$ exceeds significantly unity, where $\gamma_s^2 = (\alpha Z_s)^2 + 1$ and $\alpha$ is the usual fine structure

<footnote>$^\dagger$The mean screened ion charge is defined as $\bar{Z}_{0,s} = \sum_i f_{0,s,i} Z_{0,s,i}$ where $f_{0,s,i}$ is the local fraction of all ionization states $Z_{0,s,i}$. By definition, $\sum_i f_{0,s,i} = 1$, $Z_{0,s,i} \in \{0, Z_s\}$. The number of bound electrons is $N_{s,i} = Z_s - Z_{0,s,i}$ and the mean value is $\bar{N}_s = \sum_i f_{0,s,i} N_{s,i}$.</footnote>
constant. While for argon, relativistic corrections are negligible since $\gamma_s^{Ar} \approx 1.0086$, they become more significant for tungsten, as $\gamma_s^W \approx 1.136$. If the ground-state may be obtained for the field-free tungsten element using the density functional theory (DFT) or the Multi-Configuration Dirac-Hartree-Fock (MCDHF) approaches‡ implemented in GAUSSIAN and GRASP codes respectively [12, 13] §, the mean excitation energies for all ionization states, which play an important role in the inelastic electron-ion scattering processes, are still not available. Much in the same way, while screening effects on bremsstrahlung by runaway electrons may be reasonably described using a standard angular-averaged formula in the first Born approximation [18], such an approximation cannot be considered for less energetic electrons resonantly accelerated by RF waves for example, since the angular cone of emission is much larger [19]. Therefore, in order to cover continuously the range of kinetics energies from few keV to several tens MeV photons, a fully numerical integration over the electron emission angle of the cross-section differential in photon energy and in photon and electron emission angles must be carried out, which represents also a considerable numerical task.

In this context, the accurate incorporation of the atomic physics in kinetic codes, while keeping computational effort at a reasonable level, is a serious challenge, especially for describing inelastic scattering. A similar effort should concern the screening effects on the bremsstrahlung, a major moment of the non-thermal distribution function for diagnosing fast electron dynamics. The use of analytical formulas based of simple parametrized models is consequently the more suitable approach, with absolute calibration against results obtained by advanced numerical quantum codes. However, simple atomic models have usually a limited applicability, which prevents a systematic and consistent use for all quantities that must be modified for taking into account of the atomic physics. While the well-known Thomas-Fermi model [20], and its approximate formulations [21, 22] are well suited for neutral or weakly ionized atoms, they give usually less accurate results when the number $N_s$ of bound electrons is small as compared to $Z_s$. The charge density of inner shells is usually better described by a Yukawa electrostatic potential [23], which itself is not relevant for neutral or weakly ionized atoms, since the charge density fall-off is generally too sharp at large distance from the nucleus as compared to DFT or MCDHF calculations. However, the Yukawa atomic model is widely used for bremsstrahlung studies even if the target atom is neutral or weakly ionized, since this physics process involves usually deep electronic shells to calculate radiation emission [24–28]. Consequently, a unified and accurate description of the atomic electrostatic potential that can be used either in kinetic calculations or for bremsstrahlung without a significant degradation of the numerical performances of the kinetic codes, whatever plasma conditions (cold or hot) and the type of element, is of a great interest.

The purpose of this paper is therefore to propose a general and global approach for incorporating the atomic physics in electron Fokker-Planck solvers, allowing existing codes to be easily and robustly updated for realistic simulations, whatever the consequences on the fast electron dynamics, which will be the object of a separate study.

This objective is addressed by expressing the atomic electrostatic potential as a series of Yukawa potentials. It turns out to be a trade-off to keep codes fast and accurate over a wide range of electron kinetic energies, regardless the ionization state of the elements, with a reasonably

<footnote>‡In tokamak plasmas, The mean distance $\bar{d} = n^{-1/3}$ between particles, where n is the plasma density, is always much greater than the atomic radius ranging approximately between Bohr radius $a_0$ and $4 \times a_0$ for neutral atoms and less for corresponding ions. Therefore, for all elements and regardless their states, they can always be considered as in vacuum or field-free, which simplifies considerably ground-state calculations. This is not the case in inertial fusion plasmas [16].</footnote>

<footnote>§The FAC code based on the modified multiconfigurational Dirac-Hartree-Fock-Slater (MC-DHFS) procedure may be also an interesting alternative tool [17]</footnote>
correct description of the physics involving both outer and inner electron shells. This approach has already been successfully considered for modelling results of Hartree-Fock-Slater calculations for neutral elements only, from Hydrogen to Uranium, using up to three exponentials [29], and was also successfully applied for deriving an analytical formula of the bremsstrahlung with screening effects, valid from the classical to the fully relativistic limits [30, 31]. The modelling here considered, so-called multi-Yukawa (MY), is basically an extension of the Moliere's approach initially used for describing the Thomas-Fermi potential of neutral elements as a linear combination of three exponentials [32]. The great advantage of this method is the possibility to obtain easily analytical derivations for many physical quantities of interest in the first Born approximation, owing to the simple analytical expression of the Fourier transform of an exponential function in the calculation of the atomic form factor, assuming a spherical symmetry for the density of bound electrons in the ground state [29].

In Sec. 2, the multi-Yukawa atomic model is introduced, and the procedure of calibration is explained in detail for an arbitrary number of exponentials. Comparison with DFT or MCDHF calculations is presented. In Fokker-Planck calculations, the electron-ion collision operator is described by a friction vector and a diffusion tensor resulting from elastic and inelastic collisions. In presence of partially-ionized high-Z elements, inelastic electron-ion collisions must be also taken into account, since free electrons of the plasma may lose part of their kinetic energy by either atomic excitation or ionization. The latter process may have a critical impact on the early build-up of electron avalanches which play a major role in the dynamics of a runaway electron population in post disruptive plasmas [33]. The incorporation of the atomic physics in the Fokker-Planck collision operator is first detailed in Sec. 3, either for elastic and inelastic processes. The screening function describing the impact of partially ionized atoms on elastic Coulomb collisions is derived in Sec. 4, using the multi-Yukawa atomic model. The inelastic electron-ion collisions are then considered in Sec. 5, with the approximate approach based on the Bethe's formula of the electron energy losses per unit length [34, 35]. In this case, the atomic physics is described by the mean excitation energy from the ground-state which can be calculated by several methods. This quantity is also derived from a non-relativistic variational quantum approach but also with the classical Local Plasma Approximation (LPA), both using the multi-Yukawa atomic model [36, 37]. Conclusions are given in Sec. 6.

Though the impact of the screening on the bremsstrahlung can be described with the same atomic form factor, as that used for the Mott relativistic cross-section, this problem will be addressed in a separate paper.

## 2. Atomic model

### 2.1. Radial distribution of charge in the ground-state and form factor

Kinetic calculations with partially ionized high-Z atoms requires an atomic model that describes accurately the spatial distribution of the bound electrons $\rho_{Z_{0,s}}(\mathbf{r})$ in the ground state, regardless the type of atom and its level of ionization, where $\mathbf{r}$ is the distance to the nucleus. Indeed, excited states are transient and their lifetimes are generally much shorter than the mean time between two collisions in standard tokamak plasmas¶. By definition, $\int_V \rho_{Z_{0,s}} d\mathbf{r} = N_s = Z_s - Z_{0,s}$ where $V$ is a volume of reference characterizing the ion size, $N_s$ is the number of bound electrons, $Z_{0,s}$ is the screened ion charge and $Z_s$ is the atomic number. The screening effects are determined

<footnote>¶This condition may be marginally fullfiled in very cold post-disruptive plasmas.</footnote>
by evaluating the form factor $F_{Z_{0,s}}(\mathbf{q}) \equiv \int_V \exp(-i\mathbf{q} \cdot \mathbf{r}/\hbar) \rho_{Z_{0,s}} d\mathbf{r}$, where $\mathbf{q}$ is the usual recoil momentum by Coulomb collisions. Since the kinetic energy is conserved in the elastic scattering process and assuming that small-angle scattering predominates, $|\mathbf{q}| \simeq 2|\mathbf{p}|\sin(\theta/2)$, $\mathbf{p}$ being the incoming electron momentum and $\theta$ the deflection angle. Here, $\mathbf{q} = \mathbf{p}_f - \mathbf{p}_i$, where $\mathbf{p}_{i[f]} = \hbar\mathbf{k}_{i[f]}$, while $\mathbf{k}_{i[f]}$ are the wave vectors associated to the spinless wavefunctions $|i[f]\rangle = \exp(i\mathbf{k}_{i[f]} \cdot \mathbf{r})/\sqrt{V}$ of the incoming $|i\rangle$ and outgoing electron $|f\rangle$ respectively, both being considered as plane waves (1st Born approximation)$. From the Fermi's Golden rules, the relativistic Mott cross-section that describes Coulomb collisions in kinetic calculations must be modified according to the simple rule $Z_s \to Z_s - F_{Z_{0,s}}(\mathbf{q})$ in order to account for the partial atomic screening. For low energy electrons, since $\lim_{\|\mathbf{q}\| \to 0} F_{Z_{0,s}}(\mathbf{q}) = N_s$, the ion is fully screened i.e. $Z_s \to Z_{0,s} = Z_s - N_s$, while conversely, for very energetic electrons, it is fully stripped, since $\lim_{\|\mathbf{q}\| \to \infty} F_{Z_{0,s}}(\mathbf{q}) = 0$ **. Within this framework, the form factor $F_{Z_{0,s}}(\mathbf{q})$ is simply the Fourier transform of the spatial distribution of the bound electrons, whose determination is the starting point for investigating the effect of atomic screening in Fokker-Planck calculations.

## 2.2. Description of the multi-Yukawa electrostatic potential description

For the approximate formulation of the Thomas-Fermi atomic model as derived by Kirillov, et al., [22], but also for a Yukawa electrostatic potential (single exponential), the form factor may be expressed analytically in the same way according to the formula $F_{Z_{0,s}}(\bar{q}) = N_s / (1 + (\bar{q}\bar{a}_{Z_{0,s}}/2)^m)$, where $\bar{a}_{Z_{0,s}} \equiv 2a_{Z_{0,s}}/\alpha$ following the notation used in Ref. [7], and $\bar{q} = q/(m_e c)$ with $q = \|\mathbf{q}\|$. Here, $a_{Z_{0,s}}$ may be considered as an effective radius of the ion of charge $Z_{0,s}$ which depends on the chosen atomic model, $\alpha$ is the fine structure constant, $c$ the speed of light, $m_e$ is the electron rest mass. While $m = 3/2$ with $a_{Z_{0,s}} = 3N_s^{2/3}/(4Z_s)$ for the approximate Thomas-Fermi model [22], $m = 2$ and $a_{Z_{0,s}} = \lambda_{Z_{0,s}}^{-1}$ for the Yukawa electrostatic potential, where $\lambda_{Z_{0,s}}^{-1}$ is a characteristic screening length, its value being usually determined by a best fit of results obtained with advanced atomic calculations using DHFS (Dirac-Hartree-Fock-Slater) codes [24–28]††.

The cloud of bound electrons is assumed to be spherically symmetric around the nucleus, an approximation which turns out to be reasonably well satisfied for most ground states here considered. Indeed, the level of spherical symmetry can be evaluated from the matrix elements of the quadrupole moment rank-two tensor, directly obtained from DFT calculations [39], which measures essentially the deviation of the charge distribution $\rho_{Z_{0,s}}(\mathbf{r})$ from spherical symmetry. It is evaluated by a global parameter, $\Delta\Theta = |(\max(XX, YY, ZZ) - \min(XX, YY, ZZ)) / \max(XX, YY, ZZ)|$ where the diagonal elements of the tensor are $XX, YY, ZZ$. According to this simple definition, $\Delta\Theta = 0$ is corresponding to a perfect spherical symmetry for which all diagonals elements are identical. It turns out that this parameter is progressively increasing with the ionization level as shown in Fig 2. It is always much lower than 0.15 for $Z_{0,s} \le 40$, and very small for all noble gas-like electronic configurations, regardless the $Z_{0,s}$ value. Above $Z_{0,s} = 40$, some electronic configurations exhibit larger departure from spherical symmetry, but they concern primarily few values for $Z_{0,s} \ge 56$, which will be almost never found in tokamak plasmas.

$\|\text{The volume } V \text{ is chosen such that } \langle i | i \rangle = \langle f | f \rangle = 1$, a condition to have a probabilistic interpretation of the wave functions.

**The procedure is general and may be applied to all cross-sections derived within the 1st Born approximation like for bremsstrahlung.

††For a neutral atom of atomic number $Z_s$ and if the Thomas-Fermi model is used, $\lambda_{0,s} \equiv \lambda_{Z_{0,s}=0} = b_s^{-1} = 4[9\pi^2/2]^{-1/3} Z_s^{1/3} a_0^{-1} \simeq 1.13 Z_s^{1/3} a_0^{-1}$, a value frequently found in the literature [38]. Here, $b_s$ is the atomic radius in the Thomas-Fermi model.
For both approximate atomic models, calculations of the screening effects on elastic scattering in kinetic calculations can be fully performed analytically [7]. However, the derivation of the bremsstrahlung cross-section differential in photon energy and angle with partial screening effects, which requires an angular integration over the deflection angle of the scattered electron [18], cannot be carried out fully analytically with $m = 3/2$. An explicit analytical formulation can be obtained with $m = 2$ only, as demonstrated for the case of a neutral atom [30, 40]. Therefore, with the constraint to perform fast and accurate kinetic, but also bremsstrahlung calculations, based on analytical formulas with a unified atomic model, the use of the Yukawa electrostatic potential is unambiguously more appropriate. In order to keep its technical advantages without the intrinsic limitations for neutral or weakly ionized atoms, the simplest approach is to consider, instead, a generalized Yukawa potential, here named multi-Yukawa, which can be expressed as a series of exponentials,

$$4\pi\epsilon_0 r U_{Z_{0,s}}(r) = -Z_{0,s} - \sum_i A_{Z_{0,s,i}} \exp(-\lambda_{Z_{0,s,i}} r), \quad (1)$$

each of them describing accurately the charge distribution around the nucleus, either close or far from it. From the Poisson's equation $\Delta U_{Z_{0,s}} = \nabla^2 U_{Z_{0,s}} = -\rho_{Z_{0,s}}/\epsilon_0$, the radial normalized distribution of bound electrons is

$$\bar{\rho}_{Z_{0,s}}(\bar{r}) = \frac{Z_s - Z_{0,s}}{4\pi\bar{r}} \sum_i \bar{\lambda}_{Z_{0,s,i}}^2 \bar{A}_{Z_{0,s,i}} \exp(-\bar{\lambda}_{Z_{0,s,i}} \bar{r}), \quad (2)$$

where the density $\bar{\rho}_{Z_{0,s}}(\bar{r}) = \rho_{Z_{0,s}}(\bar{r}) a_0^3$ is in atomic units, $\bar{\lambda}_{Z_{0,s,i}} \equiv \lambda_{Z_{0,s,i}} a_0$ and $\bar{A}_{Z_{0,s,i}} = A_{Z_{0,s,i}} / (Z_s - Z_{0,s})$, $\bar{\lambda}_{Z_{0,s,i}}$ being the inverse of the normalized characteristic length and $\bar{A}_{Z_{0,s,i}}$ the weight of the $i^{th}$ Yukawa exponential respectively. Here, $\bar{r} \equiv r/a_0$, where $a_0$ is the classical Bohr radius. The corresponding form factor is therefore

$$F_{Z_{0,s}}(\bar{q}) = (Z_s - Z_{0,s}) \sum_i \frac{\bar{A}_{Z_{0,s,i}}}{1 + (\bar{q} \bar{a}_{Z_{0,s,i}}/2)^2}, \quad (3)$$

where $\bar{a}_{Z_{0,s,i}} \equiv 2\bar{\lambda}_{Z_{0,s,i}}^{-1}/\alpha$. By definition $\sum_i \bar{A}_{Z_{0,s,i}} = 1$, which guarantees that $F_{Z_{0,s}}(0) = Z_s - Z_{0,s} = N_s$.

Such an approach has been considered long time ago to describe the Thomas-Fermi atomic potential by Moliere using three exponentials [32]. The correspondence between coefficients $(B_i, \beta_i)$ found in the literature and $(\bar{A}_{Z_{0,s,i}}, \bar{\lambda}_{Z_{0,s,i}})$ is given in Appendix A. This method has been also used to fit the density of bound electrons calculated by a DHFS code for neutral atoms only, whose $Z_s$ value is ranging from 1 (hydrogen) to 92 (uranium) [29]. For most elements above argon approximately, three exponentials are necessary to reproduce accurately the radial distribution of charges when ionization is weak. Naturally, the analytical density given by Eq. (2) can only partially reproduce the oscillations of the DHFS linear density $4\pi\bar{r}^2\bar{\rho}_{Z_{0,s}}(\bar{r})$ associated with different inner shell contributions. However, the approximate form factor $F_{Z_{0,s}}(\bar{q})$ given by Eq. 3 remains very close to the numerical value determined from DHFS calculations, as core oscillations of the linear density have a small spatial weight, which validates the overall procedure.

### 2.3. Absolute model calibration

In the present work, the method used in Ref. [29] is generalized to all ionization states of any type of element. In this case, the effective number of exponentials used in (2) and (3) is determined by
the possibility to find a full set of positive $\bar{\lambda}_{Z_{0,s,i}}$ values. For some elements with an atomic number larger than tungsten, like gold ($Z_s = 79$), up to four exponentials can be found by the calibration procedure, but for lower $Z_s$ values, the number of exponentials never exceeds three usually, as found for tungsten. The determination of $\bar{A}_{Z_{0,s,i}}$ and $\bar{\lambda}_{Z_{0,s,i}}$ cannot be performed using a conventional least-squares fit method, because of the non-linearity of the problem and the existence of many local minima in the function to be minimized [29].

The method is consequently based on a technique of moments which guarantees the uniqueness of the solution with strict conditions, if it exists. However, the solution may not correspond to the best adjustment of the numerical atomic density. Nevertheless, as shown by the rather good agreement with quantum calculations, it is likely very close to it, by construction. The approach considered here ensures that the elastic Born cross sections practically coincide with that derived from the DFT or MCDHF calculations, because the error on the form factor is rather small, as this term is an integral of the bound electron density.

The coefficients $(\bar{A}_{Z_{0,s,i}}, \bar{\lambda}_{Z_{0,s,i}})$ of the multi-Yukawa description are determined from the condition $\langle \bar{r}^l \rangle = \langle \bar{r}^l \rangle^{num}$, $\langle \bar{r}^l \rangle^{num}$ being the moment of order $l$ calculated numerically from the density of bound electrons $\rho_{Z_{0,s}}^{num}(\bar{r})$ obtained by advanced atomic quantum codes. Here, from the multi-Yukawa density given by Eq. (2),

$$ \langle \bar{r}^l \rangle = \sum_i \bar{\lambda}_{Z_{0,s,i}}^{-l} \bar{A}_{Z_{0,s,i}} \Gamma(l+2) = (l+1)! \sum_i \bar{\lambda}_{Z_{0,s,i}}^{-l} \bar{A}_{Z_{0,s,i}} \quad (4) $$

where $\Gamma(z)$ is the Gamma function. Defining $\mathcal{R}_l \equiv \langle \bar{r}^l \rangle / (l+1)!$, a set of $2l$ equations depending upon the number of parameters $(\bar{A}_{Z_{0,s,i}}, \bar{\lambda}_{Z_{0,s,i}})$ to be determined is obtained,

$$ \mathcal{R}_l = \sum_i \bar{A}_{Z_{0,s,i}} \bar{\lambda}_{Z_{0,s,i}}^{-l} \quad (5) $$

and $(\bar{A}_{Z_{0,s,i}}, \bar{\lambda}_{Z_{0,s,i}})$ are calculated by solving the equation $\mathcal{R}_l = \mathcal{R}_l^{num}$, where

$$ \mathcal{R}_l^{num} = \frac{\langle \bar{r}^l \rangle^{num}}{(l+1)!} = \frac{1}{(l+1)!} \frac{4\pi}{(Z_s - Z_{0,s})} \int_0^\infty \bar{r}^{l+2} \bar{\rho}_{Z_{0,s}}^N (\bar{r}) d\bar{r} \quad (6) $$

the number $l$ being an integer greater than $-1$.

For a fit with a single Yukawa potential, only two terms remain, and since $\sum_i \bar{A}_{Z_{0,s,i}} = 1$, it can be deduced that $\bar{A}_{Z_{0,s,1}} = 1$ and $\bar{\lambda}_{Z_{0,s,1}} = \mathcal{R}_{-1}^N$. For two exponentials, four equations with four unknowns must be considered. By grouping the equations,

$$ \bar{\lambda}_{Z_{0,s,1}} + \bar{\lambda}_{Z_{0,s,2}} = \frac{\mathcal{R}_1^{num} - \mathcal{R}_{-1}^{num} \mathcal{R}_2^{num}}{(\mathcal{R}_1^{num})^2 - \mathcal{R}_2^{num}} \quad (7) $$

$$ \bar{\lambda}_{Z_{0,s,1}} \bar{\lambda}_{Z_{0,s,2}} = \frac{1 - \mathcal{R}_{-1}^{num} \mathcal{R}_1^{num}}{(\mathcal{R}_1^{num})^2 - \mathcal{R}_2^{num}} \quad (8) $$

and the values $\bar{\lambda}$ are therefore solutions of the quadratic equation

$$ ((\mathcal{R}_1^{num})^2 - \mathcal{R}_2^{num}) \bar{\lambda}^2 - (\mathcal{R}_1^{num} - \mathcal{R}_{-1}^{num} \mathcal{R}_2^{num}) \bar{\lambda} + (1 - \mathcal{R}_{-1}^{num} \mathcal{R}_1^{num}) = 0 \quad (9) $$

If both roots are real and positive, they correspond to $(\bar{\lambda}_{Z_{0,s,1}}, \bar{\lambda}_{Z_{0,s,2}})$ respectively, such that

$$ \bar{A}_{Z_{0,s,1}} = \frac{\mathcal{R}_{-1}^{num} - \bar{\lambda}_{Z_{0,s,2}}}{\bar{\lambda}_{Z_{0,s,1}} - \bar{\lambda}_{Z_{0,s,2}}} \quad (10) $$
and $\overline{A}_{Z_{0,s,2}} = 1 - \overline{A}_{Z_{0,s,1}}$, otherwise, a single exponential must be considered for the modelization. In this case, there is some loss of accuracy on the modeling of the bound electron density, but it has a moderate impact on the form factor, regarding its definition given by Eq. 3. This is an intrinsic limitation of this method, despite its robustness, highlighting that not all bound electron density profiles may be described by a series of multiple exponentials. It arises principally for highly ionized atoms, because highest moments $\mathcal{R}_l^{num}$ are too small as compare to the lowest ones. The weight of $\overline{\rho}_{Z_{0,s}}^N(\bar{r})$ at large $\bar{r}$ is therefore useless to identify a single exponential from this method. In this case, it is worth noting that standard non-linear techniques usually do not converge.

The procedure may be extended to three exponentials, and all $\bar{\lambda}$ values must be real and positive solutions of the polynomial equation $\bar{\lambda}^3 - X_1\bar{\lambda}^2 + X_2\bar{\lambda} - X_3 = 0$ to ensure that the atomic potential may be well described by a multi-Yukawa potential with the use of three exponentials. Otherwise, two exponentials must be considered in the modeling procedure, thus removing useless moments associated to three exponentials. Defining the vector $\mathbf{X}$ from coefficients $\{X_1, X_2, X_3\}$ of the polynomial equation in $\bar{\lambda}$,

$$ \mathbf{X} = \begin{bmatrix} X_1 \\ X_2 \\ X_3 \end{bmatrix} = \begin{bmatrix} \bar{\lambda}_{Z_{0,s,1}} + \bar{\lambda}_{Z_{0,s,2}} + \bar{\lambda}_{Z_{0,s,3}} \\ \bar{\lambda}_{Z_{0,s,1}}\bar{\lambda}_{Z_{0,s,2}} + \bar{\lambda}_{Z_{0,s,2}}\bar{\lambda}_{Z_{0,s,3}} + \bar{\lambda}_{Z_{0,s,1}}\bar{\lambda}_{Z_{0,s,3}} \\ \bar{\lambda}_{Z_{0,s,1}}\bar{\lambda}_{Z_{0,s,2}}\bar{\lambda}_{Z_{0,s,3}} \end{bmatrix} \quad (11) $$

the coefficients $\overline{A}_{Z_{0,s,i}}$ are determined from the matrix relation $\mathbf{A} = \mathbf{N}^{-1}\mathbf{M}\mathbf{X}$, where parameters $\overline{A}_{Z_{0,s,i}}$ are components of the vector

$$ \mathbf{A} = \begin{bmatrix} \overline{A}_{Z_{0,s,1}} \\ \overline{A}_{Z_{0,s,2}} \\ \overline{A}_{Z_{0,s,3}} \end{bmatrix} \quad (12) $$

with

$$ \mathbf{M} = \begin{bmatrix} 1 & -\mathcal{R}_1^{num} & \mathcal{R}_2^{num} \\ \mathcal{R}_1^{num} & -\mathcal{R}_2^{num} & \mathcal{R}_3^{num} \\ \mathcal{R}_2^{num} & -\mathcal{R}_3^{num} & \mathcal{R}_4^{num} \end{bmatrix} \quad (13) $$

and

$$ \mathbf{N} = \begin{bmatrix} \bar{\lambda}_{Z_{0,s,1}} & \bar{\lambda}_{Z_{0,s,2}} & \bar{\lambda}_{Z_{0,s,3}} \\ \frac{1}{\bar{\lambda}_{Z_{0,s,1}}} & \frac{1}{\bar{\lambda}_{Z_{0,s,2}}} & \frac{1}{\bar{\lambda}_{Z_{0,s,3}}} \\ \end{bmatrix} \quad (14) $$

Formally, it is possible to extend easily this method by recurrence to any number of exponentials. The dimensions of $\mathbf{X}$, $\mathbf{N}$ and $\mathbf{M}$ must be adjusted according to the number of exponentials, as well as the degree of the polynomial equation in $\bar{\lambda}$ to be solved, its coefficients being determined by expanding the product $\prod_i (\bar{\lambda} - \bar{\lambda}_{Z_{0,s,i}})$. In the numerical implementation of the method of moments, the possibility to find up to four exponentials has been considered. However, the larger number of exponentials is rarely found, only for few low ionization states of elements heavier than tungsten, like gold. For tungsten, the maximum number of exponentials never exceeds three, regardless its ionization state. As discussed previously, if no solution is found for a given set of multiple exponentials, a solution is searched for a number of exponentials decremented by one unity, and the procedure is repeated until a set of positive and real $\bar{\lambda}$ values is found. The case with a single exponential corresponding to the standard Yukawa potential is the ultimate solution, if a multi-Yukawa potential may not be found numerically.
## 2.4. Comparisons between multi-Yukawa model and quantum relativistic calculations

Numerical calculations of the radial profiles of the bound electrons have been performed for most of the elements that can be found in a plasma, regardless their ionization states using GAUSSIAN and GRASP codes respectively [12, 13]. They are all implemented in the LUKE suite of codes for studies of the atomic physics on fast electron dynamics in magnetized plasmas [10]. Details about the parameters used for the simulations with GAUSSIAN and GRASP codes are given in Appendix B. In Fig. 3, the radial profiles of the density of bound electrons for all ionizations states of tungsten calculated by the GAUSSIAN code are displayed. For low ionization states, the density exhibits clearly several bumps which correspond approximately to the principal quantum numbers *n* of the atomic orbitals. Excellent agreement is found between the results of the two codes for all ionization states, as shown for the neutral atom of tungsten and the ions $W^{10+}$, $W^{42+}$, $W^{56+}$ in Fig. 4. Consequently, numerical densities of reference $\rho_{Z_{0,s}}^N(\bar{r})$ of one of the two codes can be used indifferently for determining the coefficients $(\bar{A}_{Z_{0,s,i}}, \bar{\lambda}_{Z_{0,s,i}})$ of the multi-Yukawa description.

In the literature, a comparison between the DFT and simple atomic models frequently used in the publications has been carried out. Here, the radial dependencies of the density of bound electrons are evaluated for the neutral tungsten and the ion $W^{42+}$ using the Thomas-Fermi model but also the standard Yukawa one (single exponential). For the latter, two inverse screening lengths have been considered : $\lambda_{0,s}^B \approx 0.9Z_s^{0.42}a_0^{-1}$ from a fit of the Herman-Skillman potentials determined by solving the Dirac-Hartree-Fock-Slater (DHFS) equations [38] and $\lambda_{0,s}^{TF} = b_s^{-1} \approx 1.13Z_s^{1/3}a_0^{-1}$, where $b_s$ is the reference length in the Thomas-Fermi model. As shown in Fig. 5, a good quantitative agreement is observed between the DFT and the Thomas-Fermi model for the neutral tungsten as expected from the theory, while the agreement is poor with the standard Yukawa model, regardless the inverse screening length. Conversely, the agreement for $W^{42+}$ between the DFT and standard Yukawa model is better than with the Thomas-Fermi model. This highlights the fact that none of the simple models have a wide range of application for describing accurately the atomic physics in kinetic and radiation calculations, since the types of elements in the plasma may change with operating conditions, while their states of ionization can also vary considerably with the temperature of the plasma.

Using results obtained with the DFT model, the set of coefficients $(\bar{A}_{Z_{0,s,i}}, \bar{\lambda}_{Z_{0,s,i}})$ has been determined for up to four exponentials, but for elements lighter than gold, the maximum number of exponentials never exceed three, like for tungsten. The full list for all ionization states of tungsten is given as a reference in the Table 1. For very weakly ionized states, $Z_{0,s} \le 5$, three exponentials are found by the numerical procedure, because of the different slopes in the radial density, as shown in Fig. 3, while the number of exponentials is usually lower for larger $Z_{0,s}$ values, since the decrease of the radial density from the nucleus is becoming more regular and generally steeper. A comparison of the impact of the number of exponentials on the density profile for neutral tungsten is shown in Figs. 6. With the use of three exponentials, an excellent quantitative agreement is found between $\bar{\rho}_{Z_{0,s}}^N(\bar{r})$ determined by the DFT and the approximate multi-Yukawa description $\bar{\rho}_{Z_{0,s}}(\bar{r})$ at almost all radii. With a reduced number of exponentials, the agreement tends to deteriorate and is poor for a single exponential corresponding to the standard Yukawa description. It is interesting to note that the Molière's description of the Thomas-Fermi model is also in very good agreement with the results of the DFT, which is consistent with results shown in Fig. 5. Even if the radial dependence of the density from the nucleus is well reproduced, in the core of the atom, oscillations of the linear density $4\pi\bar{r}^2\bar{\rho}_{Z_{0,s}}^{num}(\bar{r})$ cannot be well reproduced by a series of exponentials, as shown in Fig. 7. Such an intrinsic limitation has already been observed in Ref. [29], where the same approach
is considered but with another reference atomic model (DHFS). Since the discrepancy occurs for $\bar{r} < 0.5$, its impact on the form factor remains however small, as displayed in Fig. 8, as far as the normalized recoil momentum $\bar{q}$ is less than 0.1. Knowing that most of the Coulomb collisions occur principally for very low $\bar{q}$ values corresponding to the first Born approximation, the multi-Yukawa description is therefore remarkably robust for very weakly ionized atoms.

For higher ionization states, the advantage of the multi-Yukawa description is that it remains accurate for describing both the density and the form factor. An example is given in Fig. 9 for $W^{42+}$. In this case, the agreement is very good either for two or three exponentials, while the standard Yukawa description with a single exponential has rather a poor agreement with the radial dependence of the density determined with the DFT. However, as expected, the discrepancy is less pronounced as compared to the case of the neutral atom due to the small remaining bumps of the radial density. As for the neutral atom case, the small oscillations of the numerical linear density $4\pi\bar{r}^2\bar{\rho}_{Z_{0,s}}^{num}(\bar{r})$ are not well reproduced by the series of exponentials, as shown in Fig. 10, but the departure from $\bar{\rho}_{Z_{0,s}}^{num}(\bar{r})$ has again a very small impact on the form factor (see Fig. 11).

In Table 2, the set of values $(\bar{A}_{Z_{0,s,i}}, \bar{\lambda}_{Z_{0,s,i}})$ obtained for the neutral tungsten from DFT and DHFS methods are given for comparison [29]. The coefficients of the Molière's method are also reported [32]. Even if the methodology is similar to the one detailed in Ref. [29], the coefficients for the three exponentials case have a quite significant difference. Nevertheless, as shown in Fig. 13, the approximate multi-Yukawa linear densities remain fairly close to the value obtained by the DFT, even if all the oscillations cannot be well reproduced. The differences between the coefficients are resulting from their large sensitivity due to small changes, illustrating the ill-conditionned nature of the problem here addressed. This justifies a posteriori the chosen method, as compared to a standard least-squares fit procedure which cannot converge when multiple close solutions exist [41]. In Fig. 12, the usual coefficient of determination $R^2$ is displayed for all the ionization states of tungsten to illustrate how well the multi-Yukawa reproduces the results of DFT for the atomic density$^{‡‡}$. The fact that $R^2 > 0.44$ regardless the number of exponentials indicates that the MY is an appropriate simplified atomic model for describing quantum code outputs. As expected, with a single exponential, $R^2$ gradually increases up to unity for the hydrogen-like atom, indicating that the standard Yukawa model is more appropriate for highly ionized atoms. It illustrates also the limits of this model for very weakly ionized atoms, which justified the need of the multi-Yukawa description. Conversely, for up to three exponentials, $R^2$ remains always upper than 0.8, and often close to unity even for the neutral or weakly ionized atom. Between $Z_{0,s} = 45$ and $Z_{0,s} = 55$, but also in the interval $Z_{0,s} = 64 - 68$, the method of moments is not able to identify a set of two or three exponentials, because the density $\bar{\rho}_{Z_{0,s}}(\bar{r})$ fall-off with the distance from the nucleus has a nearly single exponential dependence. In this case, the $R^2$ coefficient is lower, as shown in Fig. 12, and the density of bound electrons is therefore less accurately described as compared to other ionization states. Nevertheless, detailed calculations have shown that the impact on the atomic form factor remains still moderate. The fact that $R^2 \ge 0.7$ indicates that the MY model still remains well consistent with the results of DFT or MCDHF codes, even for these more difficult cases.

$^{‡‡}$The coefficient of determination is calculated according to the standard formula $R^2 = 1 - SS_{res}/SS_{tot}$ where $SS_{res} = \sum_j (\bar{\rho}_{Z_{0,s}}^{num}(\bar{r}_j) - \bar{\rho}_{Z_{0,s}}(\bar{r}_j))^2$ is the residual sum of squares and $SS_{tot} = \sum_j (\bar{\rho}_{Z_{0,s}}^{num}(\bar{r}_j) - \langle \bar{\rho}_{Z_{0,s}}^{num} \rangle)^2$ is the total sum of squares. The sum is performed on all radial locations $\bar{r}_j$. Here $\langle \bar{\rho}_{Z_{0,s}}^{num} \rangle = (1/N_{gridpoint}) \sum_j \bar{\rho}_{Z_{0,s}}^{num}(\bar{r}_j)$ is the mean of the calculated bound electron density. With this definition, $0 \le R^2 \le 1$ if the model is consistent with numerical data. Zero indicates a very poor agreement, and if $R^2 = 1$, the agreement is perfect.
In order to illustrate the capability of the method to identify a best fit with more than three exponentials, the case of the weakly ionized gold atom $Au^{1+}$ is shown in Fig. 13. Some differences can be seen in the inner part of the bound electron density between three and four exponentials, especially when $\bar{r} < 0.5$, but they have no impact on the atomic form factor. The coefficient of determination increases from $R^2 = 0.4223$ for a single exponential to $R^2 = 0.6597$ for two exponentials, $R^2 = 0.8546$ for three and finally to $R^2 = 0.9646$ for four exponentials, indicating that the multi-Yukawa solution with the highest possible number of exponentials gives a better agreement with DFT or MCDHF calculations.

Though the Yukawa model fails to provide an accurate description of quantum relativistic calculations, except for rather highly stripped atoms, it is interesting to evaluate, with the procedure here detailed, the ratio $\lambda^2_{Z_{0,s}}/\lambda^2_{0,s} = \varphi_s(x)$ describing the relative change of the screening length with the normalized ionization state, i.e., $x = Z_{0,s}/Z_s$. For the approximate Thomas-Fermi model derived by Kirillov et al. [22], $\varphi_s(x) = (1-x)^{-4/3}$, while from a fit of DHFS calculations, an heuristic dependency of the form $\varphi_s(x) = (1-x^{n_s+1})/(1-x)$ was guessed with $n_s = Z_s(1/3 - 0.0020 \times Z_s)$ [26, 28]. As shown in Fig. 14, $\lambda^2_{Z_{0,s}}/\lambda^2_{0,s}$ is an increasing function of $Z_{0,s}/Z_s$ whose order of magnitude is reasonably well reproduced by the approximate Thomas-Fermi model. If the heuristic dependence in Ref. [26, 28] is far from the numerical calculations based on DFT, it is found that a modified formulation $\varphi_s(x) = (1-x^{n_s+1})/(1-x)^{3/2}$ gives a much better agreement. This improvement is valid for all elements, whatever $Z_s$, in particular for light elements. Nevertheless, it should not hide the fact that the Yukawa model with a single exponential is not appropriate for describing the density of bound electrons for weakly ionized elements, even if the variation of the relative quantity $\lambda^2_{Z_{0,s}}/\lambda^2_{0,s}$ with $Z_{0,s}/Z_s$ can be reasonably well reproduced.

### 3. Generalized electron-ion collision operator

#### 3.1. Elastic scattering

The Coulomb collision operator in kinetic calculations may be expressed formally as $d\bar{f}_e/dt|_{coll} \equiv \sum_s \sum_{Z_{0,s}} \mathbf{C}_{e,Z_{0,s}}(f_e, f_{Z_{0,s}}) + \mathbf{C}_{e,e}(f_e, f_e)$ where $\mathbf{C}_{e,Z_{0,s}}(f_e, f_{Z_{0,s}})$ describes the interactions between the momentum distribution function $f_e(t, \mathbf{x}, \mathbf{p})$ of test electrons and the momentum distribution function $f_{Z_{0,s}}(t, \mathbf{x}, \mathbf{p}_s)$ of atoms of species $s$ with an ionization state $Z_{0,s}$, while $\mathbf{C}_{e,e}(f_e, f_e)$ is the linearized electron-electron collision operator§§ [10]. Here, all ionization states present in the plasma must be considered, with $Z_{0,s}$ ranging from zero for the neutral atom to $Z_s$ for the fully stripped one. The density $n_{Z_{0,s}}(t, \mathbf{x}) = \int f_{Z_{0,s}}(t, \mathbf{x}, \mathbf{p}_s) d^3\mathbf{p}_s$ of ions with a net charge $Z_{0,s}$ at the spatial location $\mathbf{x}$ is resulting from the local balance between ionization and recombination

§§The Fokker-Planck collision operator here used describes electron dynamics in plasmas whose temperatures are ranging from a few eV to several keV. Test electrons may be classical or relativistic when the Belaiev-Budker $e-e$ collision operator is considered. When the distortion of the distribution function from a Maxwellian represents a small fraction of the electron population, it is possible to linearize the electron-electron collision operator and by construction the Maxwellian is an eigenfunction of it. To account of self-collisions between fast electrons and the thermal bulk is particularly important for an accurate quantitative estimates of the rf-driven or Ohmic current source [42, 43]. It is an integral term usually determined to conserve particles and momentum, but not energy. Therefore, the electron temperature of the plasma must be a given parameter (by transport code for example), which is assumed to change slowly at the scale of the collision time. The full self-consistency between the bulk electron temperature and the fast electron energy losses may be obtained, but requires always the use of an external transport code. Therefore, the possible radiative cooling of the bulk electrons must be a part of the transport code, but not of the Fokker-Planck calculations themselves, thanks to the linearization. If time scales of radiative cooling of the bulk electrons and collisions are similar, the time ordering will fail, and the whole approach should be revisited.
ten times less for tungsten. Møller corrections to $\bar{u}_s$ are therefore always negligible, since the energy of these electrons cannot exceed 30 MeV, because of synchrotron radiation losses [45].

The fully screened relativistic Mott cross-section of collision between an electron and an ion of charge $Z_{s,0}^2$ is

$$ \frac{d\sigma_{e,Z_{0,s}}}{d\Omega} = Z_{0,s}^2 \frac{r_e^2 (1-x^2) \bar{p}^2 + 1}{4 \bar{p}^4 x^4} \qquad (21) $$

where $x = \sin(\theta/2)$, $\theta$ being the usual deflection angle of the electron and $r_e$ is the classical electron radius. Spin and relativistic corrections are negligible in the non-relativistic limit, and when $\bar{p}^2 \ll 1$, the Mott cross-section merges with the usual Rutherford expression [46, 47]. Since $\bar{p} = \gamma\beta$, where $\beta = \bar{v} = \sqrt{1-1/\gamma^2}$ and $\gamma$ is the Lorentz factor, for 200 keV slide-away electrons, $\bar{p}^2 \simeq 0.93$, while for 20 MeV run-away electrons, $\bar{p}^2 \simeq 1610$.

Knowing that $|\Delta\mathbf{p}| = 2|\mathbf{p}|\sin(\theta/2)$ and that the angular integral is taken over $\int d\Omega = \int_{x_{\min}}^{x_{\max}} \sin\theta d\theta \int_0^{2\pi} d\phi$, where $\phi$ is the azimuthal angle in the center of mass frame, the fully screened friction vector $\mathbf{A}_{Z_{0,s}}$ and diffusion tensor $\mathbb{D}_{Z_{0,s}}$ are

$$ \mathbf{A}_{Z_{0,s}} = -\Upsilon \int d^3\bar{\mathbf{v}}_s f_{Z_s}(t, \mathbf{x}, \bar{\mathbf{v}}_s) \frac{\hat{\mathbf{u}}_s}{\bar{u}_s^2} \int_{x_{\min}}^{x_{\max}} dx \frac{(1-x^2)\bar{p}^2 + 1}{x} Z_{0,s}^2 \qquad (22) $$

and

$$ \mathbb{D}_{Z_{0,s}} = \frac{1}{2}\Upsilon \int d^3\bar{\mathbf{v}}_s f_{Z_{0,s}}(t, \mathbf{x}, \bar{\mathbf{v}}_s) \frac{1}{\bar{u}_s} (\mathbf{I} - \hat{\mathbf{u}}_s\hat{\mathbf{u}}_s) \times \int_{x_{\min}}^{x_{\max}} dx \frac{(1-x^2)\bar{p}^2 + 1}{x} (1-x^2) Z_{0,s}^2 \qquad (23) $$

where $\Upsilon = 4\pi r_e^2 c$, while $x_{[\min] \max} = (1 + \bar{b}_{[\max] \min})^{-1/2}$, $\bar{b} = b/b_{90}$ being the normalized impact parameter of the Coulomb collision, with respect to the perpendicular deflection impact parameter $b_{90} = r_e Z_{0,s}^2 \beta^{-2}$. Here, the electron velocity normalized to the speed of light is linked to the Lorentz factor by the relation $\gamma = (1 - \beta^2)^{-1/2}$. The values of $\bar{b}_{[\max] \min}$ are discussed in the Appendix C. In Eq. (23), the term $(\mathbf{I} - \hat{\mathbf{u}}_s\hat{\mathbf{u}}_s)$ is the usual perpendicular collision operator.

Since both definite integrals $\int_{x_{\min}}^{x_{\max}} \dots dx/x$ in Eqs (22) and (23) give the same value, one obtains,

$$ \int_{x_{\min}}^{x_{\max}} \dots \frac{dx}{x} \approx \gamma^2 \int_{x_{\min}}^{x_{\max}} \frac{dx}{x} = (\bar{p}^2 + 1) \ln \frac{\bar{b}_{\max}}{\sqrt{1 + \bar{b}_{\min}^2}} \qquad (24) $$

where $\ln \Lambda_{e,Z_{0,s}} = \ln (\bar{b}_{\max}/\sqrt{1 + \bar{b}_{\min}^2})$ is the Coulomb logarithm. Therefore, taking $x_{\min} = 1/\Lambda_{e,Z_{0,s}}$ and $x_{\max} = 1$, Eqs. (22) and (23) may be approximated by

$$ \mathbf{A}_{Z_{0,s}} \approx -\Upsilon \int d^3\bar{\mathbf{v}}_s f_{Z_{0,s}}(t, \mathbf{x}, \bar{\mathbf{v}}_s) \frac{\hat{\mathbf{u}}_s}{\bar{u}_s^2} (\bar{p}^2 + 1) \int_{1/\Lambda_{e,Z_{0,s}}}^{1} Z_{0,s}^2 \frac{dx}{x} \qquad (25) $$

and

$$ \mathbb{D}_{Z_{0,s}} \approx \frac{1}{2}\Upsilon \int d^3\bar{\mathbf{v}}_s f_{Z_{0,s}}(t, \mathbf{x}, \bar{\mathbf{v}}_s) \frac{1}{\bar{u}_s} (\mathbf{I} - \hat{\mathbf{u}}_s\hat{\mathbf{u}}_s) (\bar{p}^2 + 1) \int_{1/\Lambda_{e,Z_{0,s}}}^{1} Z_{0,s}^2 \frac{dx}{x} \qquad (26) $$
The calculation of the Coulomb logarithm ln $\Lambda_{e,Z_{0,s}}$ is discussed in the Appendix C.

The partial screening is taken into account by replacing $Z_{0,s}^2 \rightarrow |Z_s - F_{Z_{0,s}}(\bar{q})|^2$ in Eqs. (25) and (26), where $\bar{q} = 2\bar{p}\sin(\theta/2) = 2\bar{p}x$, and following the definition in Ref. [7], the Fokker-Planck screening function $g_{Z_{0,s}}(\bar{p})$ is defined as

$$ \int_{1/\Lambda_{e,Z_{0,s}}}^{1} |Z_s - F_{Z_{0,s}}(\bar{q})|^2 \frac{dx}{x} = Z_{0,s}^2 \ln \Lambda_{e,Z_{0,s}} + g_{Z_{0,s}}(\bar{p}) \quad (27) $$

or

$$ g_{Z_{0,s}}(\bar{p}) = \int_{1/\Lambda_{e,Z_{0,s}}}^{1} \left( |Z_s - F_{Z_{0,s}}(\bar{q})|^2 - Z_{0,s}^2 \right) \frac{dx}{x} \quad (28) $$

since

$$ Z_{0,s}^2 \int_{1/\Lambda_{e,Z_{0,s}}}^{1} \frac{dx}{x} = Z_{0,s}^2 \ln \Lambda_{e,Z_{0,s}} \quad (29) $$

The formulation (28) guarantees that $\int_{1/\Lambda_{e,Z_{0,s}}}^{1} |Z_s - F_{Z_{0,s}}(\bar{q})|^2 dx/x = Z_{0,s}^2 \ln \Lambda_{e,Z_{0,s}}$ for weakly energetic electrons, while conversely, for very energetic ones, it is $Z_s^2 \ln \Lambda_{e,Z_s}$. However, as pointed out in Ref. [7], partial screening cannot be described in a strict Fokker-Planck sense other than in the complete and no screening limits. In order to keep dominant screening terms and avoid unphysical behaviour for partial screening, only terms to the lowest order in $x$ must be considered, which allows $\bar{q}$ to be significant for large electron energies, and take consequently the full form of $F_{Z_{0,s}}(\bar{q})$. The corresponding Fokker-Planck operator is then equivalent to the first Legendre mode of the Boltzmann operator at non-relativistic energies, and differs by a factor of order $1/\ln \Lambda_{e,Z_{0,s}}$ in the ultra-relativistic limit.

In the limit of an almost zero ion temperature, as in a post-disruptive regime in tokamaks, $f_{Z_{0,s}}(t, \mathbf{x}, \bar{v}_s) \approx n_{Z_{0,s}}(t, \mathbf{x}) \delta(\bar{v}_s) / (4\pi\bar{v}_s^2)$ where $\delta(\bar{v}_s)$ is the Dirac function, thus assuming that ions are at rest. In this case, the integration over $\mathbf{v}_s$ may be performed analytically, and expressions in Ref. [7] may be retrieved. However, the implementation of the screening effects in kinetic codes for studying standard regimes like in the LUKE code is slightly different, because of the finite ion temperature which requires to express $\mathbf{A}_{Z_{0,s}}$ and $\mathbb{D}_{Z_{0,s}}$ in term of Rosenbluth potentials, allowing a convenient conservative formulation of the collision operator. Knowing that $\nabla_{\bar{v}}(1/\bar{u}_s) = \partial(1/\bar{u}_s)/\partial\bar{v} = -\hat{\mathbf{u}}_s/\bar{u}_s^2$, $\nabla_{\bar{v}}\bar{u}_s = \hat{\mathbf{u}}_s$ and $\nabla_{\bar{v}}\nabla_{\bar{v}}u = (\mathbb{I} - \hat{\mathbf{u}}_s\hat{\mathbf{u}}_s)/\bar{u}_s$, the integral $\int d^3\mathbf{v}_s$ may be permuted with the derivatives $\partial/\partial\bar{v}$, and the term $(\bar{p}^2 + 1) \left[ \int_{1/\Lambda_{e,Z_{0,s}}}^{1} |Z_s - F_{Z_{0,s}}(\bar{q})|^2 dx/x \right]$ itself, which is independent of $\bar{u}_s$, so

$$ \mathbf{A}_{Z_{0,s}} \approx -\Upsilon (\bar{p}^2 + 1) \left[ \int_{1/\Lambda_{e,Z_{0,s}}}^{1} |Z_s - F_{Z_{0,s}}(\bar{q})|^2 \frac{dx}{x} \right] \frac{\partial}{\partial \bar{v}} H_{Z_{0,s}}(t, \mathbf{x}, \bar{v}) \quad (30) $$

where

$$ H_{Z_{0,s}}(t, \mathbf{x}, \bar{v}) = \int d^3 \bar{v}_s f_{Z_{0,s}}(t, \mathbf{x}, \bar{v}_s) \frac{1}{\bar{u}_s} \quad (31) $$

while

$$ \mathbb{D}_{Z_{0,s}} \approx \frac{1}{2} \Upsilon (\bar{p}^2 + 1) \left[ \int_{1/\Lambda_{e,Z_{0,s}}}^{1} |Z_s - F_{Z_{0,s}}(\bar{q})|^2 \frac{dx}{x} \right] \frac{\partial}{\partial \bar{v}} \frac{\partial}{\partial \bar{v}} G_{Z_{0,s}}(t, \mathbf{x}, \bar{v}) \quad (32) $$
with

$$G_{Z_{0,s}}(t, \mathbf{x}, \bar{\mathbf{v}}) = \int d^3\bar{\mathbf{v}}_s f_{Z_{0,s}}(t, \mathbf{x}, \bar{\mathbf{v}}_s) \bar{u}_s \quad (33)$$

In the coordinate system $(p, \xi, \varphi)$ used by the LUKE Fokker-Planck solver in momentum space [48], where $\xi$ is the cosine of the pitch-angle, the expression of the collision operator in terms of a divergence of the electron flux in momentum space $\nabla_{\mathbf{p}} \cdot \mathbf{S}_{\mathbf{p}}^{\text{coll}}(f_e)$ is

$$S_{\mathbf{p}}^{\text{coll}} = -D_{\mathbf{p}\mathbf{p}}^{\text{coll}} \frac{\partial f_e}{\partial p} + \frac{\sqrt{1-\xi^2}}{p} D_{\mathbf{p}\xi}^{\text{coll}} \frac{\partial f_e}{\partial \xi} + F_{\mathbf{p}}^{\text{coll}} f_e \quad (34)$$

$$S_{\xi}^{\text{coll}} = -D_{\xi p}^{\text{coll}} \frac{\partial f_e}{\partial p} + \frac{\sqrt{1-\xi^2}}{p} D_{\xi \xi}^{\text{coll}} \frac{\partial f_e}{\partial \xi} + F_{\xi}^{\text{coll}} f_e \quad (35)$$

assuming a local axisymmetric plasma. By symmetry, $D_{p\xi}^{\text{coll}} = D_{\xi p}^{\text{coll}} = F_{\xi}^{\text{coll}} = 0$. The contribution of elastic electron-ion collisions in the non-zero diffusion and friction terms $D_{pp}^{\text{coll}}$ and $F_{p}^{\text{coll}}$ which describes momentum exchange between particles is always very small as compared to the one of the electron-electron collisions, because of the very large difference of mass between electrons and ions, and may be therefore neglected. The single large non-zero term arising from electron-ion collision is $D_{\xi\xi}^{\text{coll}}$, which is proportional to $Z_{0,s}^2 \ln \Lambda_{e, Z_{0,s}}$. Consequently, introducing the partial screening in kinetic calculations requires simply to make the transformation $Z_{0,s}^2 \ln \Lambda_{e, Z_{0,s}} \to Z_{0,s}^2 \ln \Lambda_{e, Z_{0,s}} + gZ_{0,s} (\bar{p})$ for the pitch-angle diffusion $D_{\xi\xi}^{\text{coll}}$, with a careful account for the Coulomb logarithm $\ln \Lambda_{e, Z_{0,s}}$ with $Z_{0,s}$, as discussed in the Appendix C. Here, $D_{\xi\xi}^{\text{coll}}$ incorporates the contribution of all ion species present in the plasma and their respective ionization states.

### 3.2. Inelastic scattering

In presence of partially ionized high-Z atoms in the plasma, energetic electrons may lose a part of their kinetic energy by interacting with the bound electrons of a partially ionized atom which jumps consequently in a transient excited state. The slowing-down process which is taken into account by Fokker-Planck calculations is therefore the sum of multiple terms, the usual one from $e-e$ collisions, described in the LUKE kinetic code by the relativistic Belaiev-Budker collision operator [10, 49], the Abraham-Lorentz-Dirac reaction force for very energetic electrons arising from synchrotron radiation losses [45] and the new one from $e-i$ excitation. The latter can be deduced from the Bethe's stopping-power formula describing the losses of energy $dE$ per unit length $dx$ [35, 50]

$$-\left. \frac{dE}{dx} \right|_{Z_{0,s}} = 4\pi r_e^2 n_{Z_{0,s}} (Z_s - Z_{0,s}) \frac{m_e c^2}{\beta^2} [\ln B_{Z_{0,s}} - \beta^2] \quad (36)$$

with

$$B_{Z_{0,s}} = \frac{\sqrt{2\gamma\beta\sqrt{\gamma-1}}}{\langle \hbar\omega_{Z_{0,s}} \rangle / m_e c^2} \quad (37)$$

where $\langle \hbar\omega_{Z_{0,s}} \rangle$ is the mean excitation energy for the ion of net charge $Z_{0,s}$. Since the energy loss over a distance $\Delta x$ is equivalent of the work of an effective drag force $F_{p}^{\text{excitation}}(p)$ over that distance, its expression in Fokker-Planck calculations in simply

$$F_{p}^{\text{excitation}}(p) = -\sum_{s} \sum_{Z_{0,s}} \left. \frac{dE}{dx} \right|_{Z_{0,s}} \quad (38)$$
The validity of the Bethe slowing down formula holds principally for fast electrons, whose kinetic energy $E$ is much larger than the mean excitation energy $\langle \hbar\omega_{Z_{0,s}} \rangle$. In that case, the logarithm term always predominates over the small spin corrections given the $-\beta^2$ term or possible very small additional terms. For electrons whose kinetic energy is becoming low as compared to the mean excitation energy, the Bethe formula indicates that the stopping power tends to decrease. In this limit, the dominant inelastic term is coming from $e-e$ collisions, as Bethe formula goes to zero. Unfortunately, the Bethe formula may reverse sign, a non physical effect which is an intrinsic limitation of the Bethe's approach. This problem was identified in Ref. [7], and bypassed by performing an interpolation. In this case, $\ln B_{Z_{0,s}}$ in Eq. 36 is replaced by $\ln(1 + B_{Z_{0,s}}^{n_B})/n_B$, with $n_B$ an integer that is chosen heuristically to be 5. When $B_{Z_{0,s}} \gg 1$, $\ln(1 + B_{Z_{0,s}}^{n_B})/n_B \simeq \ln B_{Z_{0,s}}$, and the Bethe's formula is well retrieved. The Bethe-like expression guarantees that in the limit $\bar{p} = \gamma\beta \to 0$, $\ln(1 + B_{Z_{0,s}}^{n_B})/n_B - \beta^2$ is always positive and is smoothly becoming very small. The exact value of $\ln(1 + B_{Z_{0,s}}^{n_B})/n_B - \beta^2$ is not critical, since the dominant inelastic term is from $e-e$ collisions. Another approach, is also to enforce inelastic collisions from excitation of high-Z elements to zero, when $\ln B_{Z_{0,s}} - \beta^2$ becomes negative***. Both methods are equivalent numerically.

Regarding the formulation of the Fokker-Planck solver in the LUKE code as shown Eqs. (34, 35) in [10, 45, 48], the drag force $F_p^{excitation}(p)$, as given by Eq. (38), may be readily incorporated in $F_p^{coll}$.

### 4. Fokker-Planck screening function

From the definition of the Fokker-Planck screening function given by Eq. (28), and making the change of variable $y = \hat{q}/x = 2\bar{p}/\alpha$, where $\bar{p} = p/(mc)$, $g_{Z_{0,s}}(\bar{p})$ may be expressed as the sum of two terms $g_{Z_{0,s,1}}(\bar{p})$ and $g_{Z_{0,s,2}}(\bar{p})$ where

$$g_{Z_{0,s,1}}(\bar{p}) = 2Z_s \int_{y/\Lambda}^{y} (N_s - F_{Z_{0,s}}(\hat{q})) \frac{d\hat{q}}{\hat{q}} \quad (39)$$

and

$$g_{Z_{0,s,2}}(\bar{p}) = \int_{y/\Lambda}^{y} \left( F_{Z_{0,s}}^{2}(\hat{q}) - N_{s}^{2} \right) \frac{d\hat{q}}{\hat{q}} \quad (40)$$

Here, $\Lambda \equiv \Lambda_{Z_{0,s}}$, in order to simplify notations.

The form factor given by Eq. (3) may be recast in the simple form

$$F_{Z_{0,s}}(\hat{q}) = N_s \sum_i \frac{\bar{A}_{0,s,i}}{1 + (\hat{q}a_{Z_{0,s,i}})^2} \quad (41)$$

since $a_{Z_{0,s,i}} = \alpha\bar{a}_{Z_{0,s,i}}/2$ and $\alpha\hat{q} = \bar{q}$. Therefore

$$g_{Z_{0,s,1}}(\bar{p}) = 2Z_s \int_{y/\Lambda}^{y} \left( N_s - N_s \sum_i \frac{\bar{A}_{0,s,i}}{1 + (\hat{q}a_{Z_{0,s,i}})^2} \right) \frac{d\hat{q}}{\hat{q}} \quad (42)$$

***For low energy electrons, $\gamma \simeq 1 + \beta^2/2$, such that $\sqrt{\gamma-1} \simeq \beta/\sqrt{2}$. In this limit $B_{Z_{0,s}} \simeq 2E/\langle \hbar\omega_{Z_{0,s}} \rangle$. The term $\ln B_{Z_{0,s}} - \beta^2$ becomes negative if $\ln B_{Z_{0,s}} < \beta^2$ which leads to a transcendental equation in $E$ for determining this threshold.
or

$$g_{Z_{0,s},1}(\bar{p}) = Z_s(Z_s - Z_{0,s}) \sum_i \bar{A}_{0,s,i} \ln \left(1 + (\bar{p}\bar{a}_{Z_{0,s,i}})^2\right) \quad (43)$$

assuming that the condition $2\bar{p}a_{Z_{0,s,i}}/(\alpha\Lambda) = \bar{p}a_{Z_{0,s,i}}/\Lambda \ll 1$ holds. It is always valid in tokamak plasmas since the Debye sphere has a large number of particles.

Much in the same way,

$$g_{Z_{0,s},2}(\bar{p}) = + N_s^2 \sum_i \bar{A}_{0,s,i}^2 \tilde{g}_{Z_{0,s},2,i}(\bar{p}) + N_s^2 \sum_i \sum_{j \neq i} \bar{A}_{0,s,i} \bar{A}_{0,s,j} \tilde{g}_{Z_{0,s},2,i,j}(\bar{p}) \quad (44)$$

where

$$\tilde{g}_{Z_{0,s},2,i}(\bar{p}) = \int_{y/\Lambda}^{y} \frac{1 - (1 + (\hat{q}a_{Z_{0,s,i}})^2)^2}{(1 + (\hat{q}a_{Z_{0,s,i}})^2)^2} \frac{d\hat{q}}{\hat{q}} \quad (45)$$

and

$$\tilde{g}_{Z_{0,s},2,i,j}(\bar{p}) = \int_{y/\Lambda}^{y} \frac{1 - (1 + (\hat{q}a_{Z_{0,s,i}})^2) (1 + (\hat{q}a_{Z_{0,s,j}})^2)}{(1 + (\hat{q}a_{Z_{0,s,i}})^2) (1 + (\hat{q}a_{Z_{0,s,j}})^2)} \frac{d\hat{q}}{\hat{q}} \quad (46)$$

which can be integrated analytically such that

$$\tilde{g}_{Z_{0,s},2,i}(\bar{p}) = -\frac{1}{2} \frac{(\bar{p}\bar{a}_{Z_{0,s,i}})^2}{1 + (\bar{p}\bar{a}_{Z_{0,s,i}})^2} - \frac{1}{2} \ln \left(1 + (\bar{p}\bar{a}_{Z_{0,s,i}})^2\right) \quad (47)$$

while

$$\tilde{g}_{Z_{0,s},2,i,j,2}(\bar{p}) = -\frac{1}{4} \ln \left| 1 + \bar{p}^2 (\bar{a}_{Z_{0,s,i}}^2 + \bar{a}_{Z_{0,s,j}}^2) + \bar{p}^4 \bar{a}_{Z_{0,s,i}}^2 \bar{a}_{Z_{0,s,j}}^2 \right| \\ + \frac{\bar{a}_{Z_{0,s,i}}^2 + \bar{a}_{Z_{0,s,j}}^2}{2\sqrt{\Delta}} \left[ \ln \left| \frac{4\bar{a}_{Z_{0,s,i}}^2 \bar{a}_{Z_{0,s,j}}^2 \bar{p}^2 + 2(\bar{a}_{Z_{0,s,i}}^2 + \bar{a}_{Z_{0,s,j}}^2) - \sqrt{\Delta}}{4\bar{a}_{Z_{0,s,i}}^2 \bar{a}_{Z_{0,s,j}}^2 \bar{p}^2 + 2(\bar{a}_{Z_{0,s,i}}^2 + \bar{a}_{Z_{0,s,j}}^2) + \sqrt{\Delta}} \right| \\ - \ln \left| \frac{2(\bar{a}_{Z_{0,s,i}}^2 + \bar{a}_{Z_{0,s,j}}^2) - \sqrt{\Delta}}{2(\bar{a}_{Z_{0,s,i}}^2 + \bar{a}_{Z_{0,s,j}}^2) + \sqrt{\Delta}} \right| \right] \quad (48)$$

with

$$\sqrt{\Delta} = 2 \left| \bar{a}_{Z_{0,s,i}}^2 - \bar{a}_{Z_{0,s,j}}^2 \right| \quad (49)$$

Gathering all terms,

$$g_{Z_{0,s}}(\bar{p}) = Z_s(Z_s - Z_{0,s}) \sum_i \bar{A}_{0,s,i} \ln (1 + (\bar{p}\bar{a}_{Z_{0,s,i}})^2) \\ - \frac{(Z_s - Z_{0,s})^2}{2} \sum_i \bar{A}_{0,s,i}^2 \left( \frac{(\bar{p}\bar{a}_{Z_{0,s,i}})^2}{1 + (\bar{p}\bar{a}_{Z_{0,s,i}})^2} + \ln (1 + (\bar{p}\bar{a}_{Z_{0,s,i}})^2) \right) \\ + 2(Z_s - Z_{0,s})^2 \sum_{i>j>i} \sum_i \bar{A}_{0,s,i} \bar{A}_{0,s,j} \tilde{g}_{Z_{0,s},2,i,j}(\bar{p}) \quad (50)$$
with

$$
\begin{equation}
\begin{aligned}
\tilde{g}_{Z_{0,s},2,i,j}(\bar{p}) = & -\frac{1}{2} \frac{a_{Z_{0,s,i}}^2 + a_{Z_{0,s,j}}^2}{a_{Z_{0,s,j}}^2 - a_{Z_{0,s,i}}^2} \ln \left| \frac{1 + \bar{p}^2 \bar{a}_{Z_{0,s,j}}^2}{1 + \bar{p}^2 \bar{a}_{Z_{0,s,i}}^2} \right| \\
& -\frac{1}{4} \ln \left| 1 + \bar{p}^2 \left( \bar{a}_{Z_{0,s,i}}^2 + \bar{a}_{Z_{0,s,j}}^2 \right) + \bar{p}^4 \bar{a}_{Z_{0,s,i}}^2 \bar{a}_{Z_{0,s,j}}^2 \right| \\
& + \frac{\bar{a}_{Z_{0,s,i}}^2 + \bar{a}_{Z_{0,s,j}}^2}{4 \left| \bar{a}_{Z_{0,s,i}}^2 - \bar{a}_{Z_{0,s,j}}^2 \right|} \\
& \times \left[ \ln \left| \frac{2 \bar{a}_{Z_{0,s,i}}^2 \bar{a}_{Z_{0,s,j}}^2 \bar{p}^2 + \left( \bar{a}_{Z_{0,s,i}}^2 + \bar{a}_{Z_{0,s,j}}^2 \right) - \left| \bar{a}_{Z_{0,s,i}}^2 - \bar{a}_{Z_{0,s,j}}^2 \right|}{2 \bar{a}_{Z_{0,s,i}}^2 \bar{a}_{Z_{0,s,j}}^2 \bar{p}^2 + \left( \bar{a}_{Z_{0,s,i}}^2 + \bar{a}_{Z_{0,s,j}}^2 \right) + \left| \bar{a}_{Z_{0,s,i}}^2 - \bar{a}_{Z_{0,s,j}}^2 \right|} \right| \\
& - \ln \left| \frac{\left( \bar{a}_{Z_{0,s,i}}^2 + \bar{a}_{Z_{0,s,j}}^2 \right) - \left| \bar{a}_{Z_{0,s,i}}^2 - \bar{a}_{Z_{0,s,j}}^2 \right|}{\left( \bar{a}_{Z_{0,s,i}}^2 + \bar{a}_{Z_{0,s,j}}^2 \right) + \left| \bar{a}_{Z_{0,s,i}}^2 - \bar{a}_{Z_{0,s,j}}^2 \right|} \right| \right]
\end{aligned}
\tag{51}
\end{equation}
$$

and, as expected, $\lim_{\bar{p} \to 0} g_{Z_{0,s}}(\bar{p}) = 0$ is verified regardless the element and its ionization state $Z_{0,s}$.

In the case of a single exponential corresponding to the standard Yukawa atomic potential, Eq. (50) simplifies to the usual form

$$
g_{Z_{0,s}} = \frac{1}{n} \left[ (Z_s^2 - Z_{0,s}^2) \ln \left(1 + (\bar{p}\bar{a}_{Z_{0,s,1}})^n\right) - (Z_s - Z_{0,s})^2 \frac{(\bar{p}\bar{a}_{Z_{0,s,1}})^n}{1 + (\bar{p}\bar{a}_{Z_{0,s,1}})^n} \right] \quad (52)
$$

where $n = 2$, since $\bar{A}_{0,s,1} = 1$, while $\tilde{g}_{Z_{0,s,2,i,j}}(\bar{p}) = 0$ by definition. Using the Thomas-Fermi-Kirillov model [22], $n = 3/2$, and Eq. 6 in Ref. [6] is well retrieved. The difference is generally small for tungsten, few percent, between the MY and Thomas-Fermi-Kirillov models.

The analytical expression of $g_{Z_{0,s}}(\bar{p})$ may be easily implemented in Fokker-Planck solvers, allowing fast and accurate kinetic calculations, whatever $\bar{p}$ and the type of elements and its level of ionization. As shown in Fig. 15, $g_{Z_{0,s}}(\bar{p})$ from Eq. (50) for neutral tungsten with the use of three exponentials is very close to the numerical estimate $g_{Z_{0,s}}^{num}(\bar{p})$ directly obtained in the limit $\bar{p} \gg 1$, from DFT calculations using results of the GAUSSIAN code. Its formula is

$$
g_{Z_{0,s}}^{num}(\bar{p}) = (Z_s^2 - Z_{0,s}^2) (\ln(2\bar{p}/\alpha) + \gamma_{EM} - 1) + 2Z_s N_s \hat{I}_{1,Z_{0,s}} + N_s^2 \left(\frac{1}{2} - \hat{I}_{2,Z_{0,s}}\right) \quad (53)
$$

where

$$
\hat{I}_{1,Z_{0,s}} \equiv \frac{4\pi}{N_s} \int_0^\infty \bar{\rho}_{Z_{0,s}}^{num}(\bar{r}_1) \bar{r}_1^2 \ln \bar{r}_1 d\bar{r}_1 \quad (54)
$$

with

$$
\hat{I}_{2,Z_{0,s}} = \frac{4\pi^2}{N_s^2} \int_0^\infty \bar{\rho}_{Z_{0,s}}^{\text{num}}(\bar{r}_1) \hat{J}_{2,Z_{0,s}}(\bar{r}_1) \bar{r}_1 d\bar{r}_1 \quad (55)
$$

and

$$
\hat{J}_{2,Z_{0,s}}(\bar{r}_1) = \int_0^\infty \bar{\rho}_{Z_{0,s}}^{\text{num}}(\bar{r}_2) \bar{r}_2 d\bar{r}_2 \left( (\bar{r}_1 + \bar{r}_2)^2 \ln(\bar{r}_1 + \bar{r}_2) - (\bar{r}_1 - \bar{r}_2)^2 \ln|\bar{r}_1 - \bar{r}_2| \right) \quad (6)
$$

from Ref. [7], $\gamma_{EM}$ being the Euler-Mascheroni constant [51].
Since Eq. (53) is derived in the limit $\bar{p} \gg 1$, it is consequently not valid at low $\bar{p}$, and $g_{Z_{0,s}}^{num}(\bar{p})$ does not converge towards zero when $\bar{p} \le 0.5$, as shown in Fig. 16 for the neutral atom of tungsten. When the number of exponentials is reduced, $g_{Z_{0,s}}(\bar{p})$ is always lower than $g_{Z_{0,s}}^{N}(\bar{p})$. The relative error is about 13% at $\bar{p} = 1$ for a single exponential. This tendency is similar for an ionized atom as shown for $W^{42+}$ in Fig. 17. In this case, only two exponentials are necessary to accurately reproduce the atomic potential. For ionization states ranging between $W^{45+}$ and $W^{55+}$, where only a single exponential can be found by the MY procedure described in Sec. 2, the relative deviation of $g_{Z_{0,s}}(\bar{p})$ from $g_{Z_{0,s}}^{num}(\bar{p})$ remains small whatever the $\bar{p}$ value, of the order of few percent, even if the coefficient of determination $R^2$, shown in Fig. 12, is lower. This results from the fact that $g_{Z_{0,s}}(\bar{p})$ is itself an integral, which smoothes out possible errors.

## 5. Mean excitation energy

The mean excitation energy $\langle \hbar\omega_{Z_{0,s}} \rangle$ is the key parameter to describe enhanced slowing down of the electrons by transferring energy to partially ionized high-Z elements in a hot plasma. It is formally defined as $\langle \hbar\omega_{Z_{0,s}} \rangle = (1/Z_{0,s}) \sum_{ik} f_{ik} \ln(\hbar\omega_{ik})$, where $f_{ik}$ is the dipole oscillator strength of the transition $\omega_{ik}$ for the atomic system between quantum states $|i\rangle$ and $|k\rangle$, according to the Bethe's theory [34, 52]. Its determination from first principles calculations is a considerable challenge, so except for elements that do not require relativistic corrections, $\langle \hbar\omega_{Z_{0,s}} \rangle$ is generally obtained from empirical laws constrained by measurements for neutral atoms only [53–56]. Recent advanced calculations carried out by a non-relativistic Multi-Configurational Self-Consistent Field (MCSCF) code have allowed to estimate $\langle \hbar\omega_{Z_{0,s}} \rangle$ for all ionization states of the elements lighter than argon $Z_s \le 18$ [14, 15, 57]. Though this result represents a considerable progress, the accurate determination of $\langle \hbar\omega_{Z_{0,s}} \rangle$ for much higher-Z elements like tungsten is still missing, which represents a difficulty for studying the impact of inelastic processes by electron-ion interaction in a hot plasma. In this context, several simple models have been introduced to compare their impact on kinetic calculations.

In general, MCSCF calculations show that $\langle \hbar\omega_{Z_{0,s}} \rangle$ has an exponential-like dependence with $Z_{0,s}$, which can be easily determined at the two limits, i.e. for the neutral atom and for the hydrogen-like atom characterized by a single valence electron. For low-Z neutral elements, $\langle \hbar\omega_{Z_{0,s}} \rangle$ has a rather complex structure, while it is becoming almost proportional to $Z_s$ for $Z_s > 18$, i.e. $\hbar\langle \omega_{Z_{0,s}} \rangle \approx 10Z_s eV$, this relation being known as the Bloch relation [53]. For $Z_s < 18$, $\langle \hbar\omega_{Z_{0,s}} \rangle$ oscillates with $Z_s$ and tends to increase up to 50% approximately as its value decreases. The $Z_s$ dependence of $\langle \hbar\omega_{Z_{0,s}} \rangle$ for neutral atoms can be well described by a statistical approach of the energy loss process, known as the Local Plasma Approximation (LPA) [37, 53]. On the other limit corresponding to a single bound electron, $\langle \hbar\omega_{Z_{0,s}} \rangle = Z_s^2 I_H eV$, where $I_H = 14.9916 eV$ is obtained from non-relativistic quantum calculations [14, 58]. Consequently, $\langle \hbar\omega_{Z_{0,s}} \rangle$ may be approximated by the simple heuristic relation

$$ \ln \langle \hbar\omega_{Z_{0,s}} \rangle \approx \ln \left( \frac{I_H}{10} Z_s \right) \frac{Z_{0,s}}{Z_s - 1} + \ln(10Z_s) \quad (57) $$

where $Z_{0,s}$ is the charge of the fully screened ion, and $\hbar\omega_{Z_{0,s}}$ is expressed in eV units. This relation can be considered as an upper bound of $\langle \hbar\omega_{Z_{0,s}} \rangle$, since electron-electron correlations tend to reduce the mean excitation energy [14]. It has been used to quantify the impact of tungsten on the toroidal plasma current driven by the RF wave at the Lower Hybrid frequency in tokamaks [4].
As shown for argon in Fig. 18, the exponential interpolation given by Eq. (57) is in good agreement for both weakly and highly ionized states as compared to MCSCF calculations [14]. In between, results obtained with the MCSCF code are lower, especially in the interval $Z_{0,s} = [10 - 15]$. However, the difference never exceeds a factor two.

More refined approaches may be considered, taking into account for the density of bound charges calculated in the ground-state. Indeed, since excited states are transient with a characteristic time generally much smaller than the mean collision time, elements in the plasma are principally in a ground-state from which atom transitions must be considered to evaluate $\langle \hbar\omega_{Z_{0,s}} \rangle$. Two models are interesting for this purpose, the LPA approach [37], dedicated principally for very weakly ionized atoms, and a variational quantum description [36]. Though restricted to non-relativistic elements, the latter may be an interesting alternative even for high-Z elements, as it is expected to be valid within a larger range of $Z_{0,s}$ values. For both approaches, the multi-Yukawa description of the density of bound electrons may be used, allowing a unified description of the atomic physics in kinetic calculations, not only for elastic Coulomb collisions but also for inelastic processes.

### 5.0.1. Local Plasma Approximation
The Local Plasma Approximation (LPA) has been widely used for calculating the mean excitation energies of neutral atoms [53]. It can be extended to any ionization state, according to the relation

$$ \ln \langle \hbar\omega_{Z_{0,s}} \rangle = \frac{4\pi}{N_s} \int_0^\infty \bar{r}^2 \bar{\rho}_{Z_{0,s}}(\bar{r}) \ln \left( \gamma_{LPA} \sqrt{4\pi\alpha^2 m_e c^2} \sqrt{\bar{\rho}_{Z_{0,s}}(\bar{r})} \right) d\bar{r} \quad (58) $$

The LPA formula gives generally poor results when the ionization state is high and a few number of electrons remain bound, because of the basic difficulties encountered when one tries to derive this scheme from first principles, i.e., starting with the standard definition of the oscillator strength in terms of dipole matrix elements and carrying out a systematic deduction [59]. In particular, the choice of $\gamma_{LPA}$ is rather arbitrary, and its value, from heuristic arguments, is generally set to $\sqrt{2}$. Incorporating Eq.(2) into Eq.(58) with $\gamma_{LPA} = \sqrt{2}$,

$$ \ln \langle \hbar\omega_{Z_{0,s}} \rangle = \sum_i \bar{\lambda}_{Z_{0,s,i}}^2 \bar{A}_{Z_{0,s,i}} \int_0^\infty \bar{r} \exp(-\bar{\lambda}_{Z_{0,s,i}} \bar{r}) \times \ln \left( \sqrt{2}\alpha^2 m_e c^2 \sqrt{\frac{Z_s - Z_{0,s}}{\bar{r}}} \sum_j \bar{\lambda}_{Z_{0,s,j}}^2 \bar{A}_{Z_{0,s,j}} \exp(-\bar{\lambda}_{Z_{0,s,j}} \bar{r}) \right) d\bar{r} \quad (59) $$

As shown for argon in Fig. 18, which is the highest-Z elements for which advanced numerical quantum calculations are available [14], a good quantitative agreement is found between LPA calculations using Eq. (59) and MCSCF ones for the neutral atom. The value from NIST database is also consistent with LPA level [60]. As the ion charge $Z_{0,s}$ is increasing, the departure from the results of the numerical quantum calculations is more and more pronounced, and with the LPA, the limit for the hydrogen-like ion is never recovered, indicating that the model fails completely in this regime. For the case of tungsten, the lack of quantum calculations prevent an accurate comparison as for argon. Nevertheless, the departure from the hydrogen-like limit is also very large, while for the neutral atom, the agreement between the value given by the NIST database and the estimate from the Bloch relation is very good, as shown in Fig. 19.
5.0.2. Variational quantum description The non-relativistic variational quantum model to calculate the mean excitation energy, initially derived for inertial fusion experiments, is an interesting approach to get a more accurate estimate $\langle \hbar\omega_{Z_{0,s}} \rangle$ as a function of $Z_{0,s}$. According to Ref. [36], the mean excitation energy from the ground state is given by the relation $\ln \langle \hbar\omega_{Z_{0,s}} \rangle = \frac{1}{2} \ln S(1)/S(-1)$ where $S(-1) = 2m_e a_0^2 \langle \bar{r}^2 \rangle / (3\hbar^2)$ and $S(1) = 4K_0/3$. The functions $S$ are moments of the strength distribution of oscillators [52], which may be expressed as a function of $K_0$, the averaged kinetic energy of the cloud of bound electrons, and $\langle \bar{r}^2 \rangle = (4\pi/N_s) \int_0^\infty \bar{r}^4 \bar{\rho}_{Z_{0,s}}(\bar{r}) d\bar{r}$ or $\langle \bar{r}^2 \rangle = 6R_2$ from Eq. (6). Therefore

$$ \langle \hbar\omega_{Z_{0,s}} \rangle^2 = 2 \frac{K_0}{\langle \bar{r}^2 \rangle} \alpha^2 m_e c^2 \qquad (60) $$

where $\alpha^2 m_e c^2 \simeq 27.21 \text{ eV}$, which is about twice the Rydberg unit of energy. The calculation of $K_0$ is performed using the virial theorem, $2K_0 \simeq -\langle U_{Z_{0,s}} \rangle$, where $U_{Z_{0,s}}$ is the atomic potential related to the density of bound charge $\rho_{Z_{0,s}}$ by the Poisson's equation, as given by Eq. (1), using the multi-Yukawa description. Therefore, since

$$ \langle U_{Z_{0,s}} \rangle = - \frac{Z_s - N_s}{4\pi\epsilon_0 a_0} \sum_i \bar{\lambda}_{Z_{0,s,i}} \bar{A}_{Z_{0,s,i}} \\ - \frac{N_s}{4\pi\epsilon_0 a_0} \left( \frac{1}{2} \sum_i \bar{A}_{Z_{0,s,i}}^2 \bar{\lambda}_{Z_{0,s,i}} \right. \\ \left. + \sum_i \sum_{j \neq i} \bar{A}_{Z_{0,s,i}} \bar{A}_{Z_{0,s,j}} \frac{\bar{\lambda}_{Z_{0,s,i}}^2}{\bar{\lambda}_{Z_{0,s,i}} + \bar{\lambda}_{Z_{0,s,j}}} \right) \qquad (61) $$

one obtains

$$ \langle \hbar\omega_{Z_{0,s}} \rangle = \frac{\alpha^2 m_e c^2}{\sqrt{6 \sum_i \bar{\lambda}_{Z_{0,s,i}}^{-2} \bar{A}_{Z_{0,s,i}}}} \\ \times \left[ (Z_s - N_s) \sum_i \bar{\lambda}_{Z_{0,s,i}} \bar{A}_{Z_{0,s,i}} + \right. \\ \left. + N_s \left( \sum_i \bar{A}_{Z_{0,s,i}}^2 \frac{\bar{\lambda}_{Z_{0,s,i}}}{2} + \sum_i \sum_{j \neq i} \bar{A}_{Z_{0,s,i}} \bar{A}_{Z_{0,s,j}} \frac{\bar{\lambda}_{Z_{0,s,i}}^2}{\bar{\lambda}_{Z_{0,s,i}} + \bar{\lambda}_{Z_{0,s,j}}} \right) \right]^{1/2} \qquad (62) $$

by reporting Eq. (61) in the expression Eq. (60). In Eq. (62), parameters $(\bar{A}_{0,s,i}, \bar{\lambda}_{Z_{0,s,i}})$ are those obtained from the method of moments discussed in Sec. 2.3.

For a single exponential, Eq. (62) simplifies to $\langle \hbar\omega_{Z_{0,s}} \rangle = \alpha^2 m_e c^2 \sqrt{\bar{\lambda}_{Z_{0,s,1}}^3 N_s / 12}$, and for neutral atoms, i.e. when $N_s = Z_s$, $\langle \hbar\omega_{Z_{0,s}} \rangle = 9.44 Z_s \text{ eV}$, considering the Thomas-Fermi model for which $\bar{\lambda}_{Z_{0,s,1}} = 1.13 Z_s^{1/3}$, as defined in Sec. 2.2. With a similar approach based on an approximate description of the Thomas-Fermi model [21], the same value is found, as shown in Ref. [36]. Both relations are very close to the heuristic Bloch relation which can be also well reproduced by the LPA model [53]. Another approximate description of the Thomas-Fermi model given in Ref. [22] gives $\langle \hbar\omega_{Z_{0,s}} \rangle \simeq 12.10 Z_s \text{ eV}$, about 12% larger than the value given by the exact Thomas-Fermi model, but still close to the Bloch relation [53].
As shown for argon in Fig. 18, the quantitative agreement between the non-relativistic variational quantum model and the results of the MCSCF code is good, especially above $Z_{0,s} = 8$. The analytical model has the correct dependency up to the hydrogen-like atom, which is an important assessment of the method. With the multi-Yukawa description of the atomic potential, small departures are observed at low $Z_{0,s}$, even if it never exceeds 20% approximately for the neutral atom. This discrepancy arises from the large sensitivity of $\langle\hbar\omega_{Z_{0,s}}\rangle$ to the atomic model in this limit. Indeed, when $\langle\hbar\omega_{Z_{0,s}}\rangle$ is calculated using the approximate Thomas-Fermi atomic models instead of the multi-Yukawa one, the agreement with the Bloch relation for neutral atoms is much better, within 5%.

When applied to the case of tungsten, as displayed in Fig. 19, $\langle\hbar\omega_{Z_{0,s}}\rangle$ exhibits globally a consistent agreement with the expected limits for a neutral atom and an almost fully stripped one. As for the argon, the agreement is less accurate near $Z_{0,s} = 0$, with a similar relative error. Conversely to the LPA model, $\langle\hbar\omega_{Z_{0,s}}\rangle$ has a correct variation with $Z_{0,s}$ when its value is close to $Z_s$, making the variational quantum description more appropriate, even if relativistic effects are not considered.

From estimates of $\langle\hbar\omega_{Z_{0,s}}\rangle$ using the LPA and the variational quantum descriptions, both using the multi-Yukawa atomic model for the ground-state, a trade-off for an accurate estimate of $\langle\hbar\omega_{Z_{0,s}}\rangle$ whatever $Z_{0,s}$ would be to consider the largest of the values given by both models. Such an approach would allow to describe accurately atomic physics in kinetic calculations, either for a cold plasma like after a major disruption, and a standard hot magnetized plasma expected during regular tokamak operation. This option is considered in the 3-D linearized relativistic bounce-averaged Fokker-Planck code LUKE for studying both the physics of post-disruptive runaway and slide-away RF-driven electrons [10, 45].

## 6. Conclusion

The incorporation of the atomic physics in kinetic calculations is becoming mandatory in order to study the impact of high-Z elements in fusion plasmas. The multi-Yukawa approach for describing the atomic potential, regardless the ionization state, is particularly convenient to obtain consistent analytical solutions for both elastic and inelastic scattering processes that occur in a plasma. It allows fast and accurate kinetic calculations while the full atomic physics can be incorporated, over a very wide range of plasma regimes and electron kinetic energies. With this approximate and accurate atomic model, the dynamics of electrons in a plasma can be studied from the runaway energy range (few ten MeV) in very cold post-disruptive plasmas to slide-away electrons in hot plasmas without changing the atomic model depending upon the studied physics.

The great advantages of the method proposed here are its robustness and flexibility. Indeed, the calibration procedure against advanced numerical atomic codes is rigorous and the parameters defining the multi-Yukawa atomic potential are unique, as their identification do not rely on a non-linear least-square fit procedure which is inappropriate for the non-linear problem here addressed. The method, initially restricted to neutral atoms, has been extended here to any ionization state of any types of elements, making it universal. While it was initially developed for up to three exponentials, it has been extended for an arbitrary number of them. However, for all the elements with an atomic number less than 74 (tungsten), the method does not find more than three exponentials, regardless their ionization states. This method is also flexible, since the impact of most advanced atomic simulations can be investigated without changing the structure of the kinetic codes, but just by modifying the coefficients of the multi-Yukawa potential. Though simple, the
method allows having a more realistic description of the atomic physics as compared to simplified
atomic models like the well-known Thomas-Fermi model and all its approximate representations.

The work presented here has been restricted to microscopic collision processes for kinetic calculations. However, it can be extended to other physical quantities, in particular those who are derived in the first Born approximation, as already shown for the bremsstrahlung of fast electrons on neutral atoms, where partial screening effect may be also important. It is likely that this method can be extended to many more processes, like ionization and knock-on collisions by energetic electrons (beyond Fokker-Planck approximation), opening the possibility of a unified and rigorous description of most of the atomic physics in kinetic descriptions of plasmas.

Acknowledgments

This work has been partially funded by the National Science Centre, Poland (NCN) grant HARMONIA 10 no. 2018/30/M/ST2/00799. We gratefully acknowledge Poland's high-performance computing infrastructure PLGrid (HPC Centers: ACK Cyfronet AGH) for providing computer facilities and support within computational grant no. PLG/2022/015994. This work has been published in the framework of the international project co-financed by the Polish Ministry of Education and Science, as program "PMW". This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 — EUROfusion). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.
<table><caption>Table 1: Table of multi-Yukawa coefficients (3 exponentials) for the different ionization states of tungsten, based on DFT calculations done with GAUSSIAN code for the density of reference [12]. Note that $\bar{A}_{W,3} = 1 - \bar{A}_{W,2} - \bar{A}_{W,1}$ by definition and in the case of two exponentials, $\bar{A}_{W,1} + \bar{A}_{W,2} = 1$.</caption><thead><tr><th>Z<sub>0</sub></th><th>A&#773;<sub>W,1</sub></th><th>A&#773;<sub>W,2</sub></th><th>&lambda;&#773;<sub>W,1</sub></th><th>&lambda;&#773;<sub>W,2</sub></th><th>&lambda;&#773;<sub>W,3</sub></th><th>Z<sub>0</sub></th><th>A&#773;<sub>W,1</sub></th><th>A&#773;<sub>W,2</sub></th><th>&lambda;&#773;<sub>W,1</sub></th><th>&lambda;&#773;<sub>W,2</sub></th><th>&lambda;&#773;<sub>W,3</sub></th></tr></thead><tbody><tr><td>0</td><td>0.0964</td><td>0.7058</td><td>40.7776</td><td>4.7464</td><td>1.2612</td><td>37</td><td>0.0657</td><td>0.9342</td><td>76.5450</td><td>9.0381</td><td>-</td></tr><tr><td>1</td><td>0.1129</td><td>0.6490</td><td>36.2792</td><td>4.8617</td><td>1.5970</td><td>38</td><td>0.0668</td><td>0.9332</td><td>75.9433</td><td>9.3208</td><td>-</td></tr><tr><td>2</td><td>0.1134</td><td>0.5416</td><td>35.9580</td><td>5.4016</td><td>2.1100</td><td>39</td><td>0.06723</td><td>0.9328</td><td>75.6537</td><td>9.6390</td><td>-</td></tr><tr><td>3</td><td>0.0872</td><td>0.4628</td><td>42.7872</td><td>6.4409</td><td>2.4903</td><td>40</td><td>0.09735</td><td>1.0166</td><td>58.6355</td><td>9.3206</td><td>6.8896</td></tr><tr><td>4</td><td>0.0447</td><td>0.3844</td><td>66.8955</td><td>8.5866</td><td>2.8794</td><td>41</td><td>0.07858</td><td>0.9464</td><td>67.6093</td><td>10.1206</td><td>6.1217</td></tr><tr><td>5</td><td>0.2404</td><td>0.7595</td><td>22.7136</td><td>3.4038</td><td>-</td><td>42</td><td>0.05895</td><td>0.9532</td><td>82.4553</td><td>10.8159</td><td>6.0350</td></tr><tr><td>6</td><td>0.2001</td><td>0.7998</td><td>25.9976</td><td>3.6903</td><td>-</td><td>43</td><td>0.04067</td><td>1.0524</td><td>108.3990</td><td>11.2870</td><td>8.6963</td></tr><tr><td>7</td><td>0.1916</td><td>0.8084</td><td>27.0553</td><td>3.8056</td><td>-</td><td>44</td><td>0.07209</td><td>-0.3479</td><td>83.3834</td><td>21.2223</td><td>13.5253</td></tr><tr><td>8</td><td>0.1691</td><td>0.8309</td><td>29.60311</td><td>4.0558</td><td>-</td><td>45</td><td>1</td><td>-</td><td>16.3046</td><td>-</td><td>-</td></tr><tr><td>9</td><td>0.1540</td><td>0.8460</td><td>31.8430</td><td>4.2397</td><td>-</td><td>46</td><td>1</td><td>-</td><td>16.7547</td><td>-</td><td>-</td></tr><tr><td>10</td><td>0.1336</td><td>0.8673</td><td>35.4275</td><td>4.4711</td><td>1.8417</td><td>47</td><td>1</td><td>-</td><td>17.1551</td><td>-</td><td>-</td></tr><tr><td>11</td><td>0.1321</td><td>0.8786</td><td>36.0567</td><td>4.5467</td><td>2.7534</td><td>48</td><td>1</td><td>-</td><td>17.5828</td><td>-</td><td>-</td></tr><tr><td>12</td><td>0.0781</td><td>0.9219</td><td>53.8146</td><td>5.0492</td><td>-</td><td>49</td><td>1</td><td>-</td><td>18.0408</td><td>-</td><td>-</td></tr><tr><td>13</td><td>0.0735</td><td>0.9265</td><td>56.8908</td><td>5.1789</td><td>-</td><td>50</td><td>1</td><td>-</td><td>18.5337</td><td>-</td><td>-</td></tr><tr><td>14</td><td>0.0703</td><td>0.9297</td><td>59.4299</td><td>5.3019</td><td>-</td><td>51</td><td>1</td><td>-</td><td>19.0654</td><td>-</td><td>-</td></tr><tr><td>15</td><td>0.0683</td><td>0.9317</td><td>61.3319</td><td>5.4187</td><td>-</td><td>52</td><td>1</td><td>-</td><td>19.6402</td><td>-</td><td>-</td></tr><tr><td>16</td><td>0.0464</td><td>0.9536</td><td>84.5389</td><td>5.7162</td><td>-</td><td>53</td><td>1</td><td>-</td><td>20.2658</td><td>-</td><td>-</td></tr><tr><td>17</td><td>0.0217</td><td>0.9783</td><td>65.8142</td><td>6.0570</td><td>-</td><td>54</td><td>1</td><td>-</td><td>20.9494</td><td>-</td><td>-</td></tr><tr><td>18</td><td>0.0195</td><td>0.9805</td><td>84.0652</td><td>6.1945</td><td>-</td><td>55</td><td>1</td><td>-</td><td>21.7004</td><td>-</td><td>-</td></tr><tr><td>19</td><td>0.0181</td><td>0.9819</td><td>98.7309</td><td>6.3284</td><td>-</td><td>56</td><td>0.01499</td><td>0.9850</td><td>383.7327</td><td>17.0332</td><td>-</td></tr><tr><td>20</td><td>0.0173</td><td>0.9827</td><td>208.1839</td><td>6.4596</td><td>-</td><td>57</td><td>0.03214</td><td>0.9679</td><td>203.7360</td><td>17.4369</td><td>-</td></tr><tr><td>21</td><td>0.0172</td><td>0.9828</td><td>211.7881</td><td>6.5890</td><td>-</td><td>58</td><td>0.05357</td><td>0.9464</td><td>139.5996</td><td>17.9070</td><td>-</td></tr><tr><td>22</td><td>0.0179</td><td>0.9821</td><td>205.6840</td><td>6.7123</td><td>-</td><td>59</td><td>0.08041</td><td>0.9196</td><td>106.6274</td><td>18.4671</td><td>-</td></tr><tr><td>23</td><td>0.0191</td><td>0.9809</td><td>195.1947</td><td>6.8342</td><td>-</td><td>60</td><td>0.1128</td><td>0.8872</td><td>87.0463</td><td>19.1924</td><td>-</td></tr><tr><td>24</td><td>0.0208</td><td>0.9791</td><td>181.9760</td><td>6.9551</td><td>-</td><td>61</td><td>0.1501</td><td>0.8499</td><td>74.3800</td><td>20.1958</td><td>-</td></tr><tr><td>25</td><td>0.0232</td><td>0.9768</td><td>166.6728</td><td>7.0732</td><td>-</td><td>62</td><td>0.1813</td><td>0.8187</td><td>67.0734</td><td>21.8502</td><td>0.0069</td></tr><tr><td>26</td><td>0.0260</td><td>0.9740</td><td>152.0133</td><td>7.1919</td><td>-</td><td>63</td><td>0.05328</td><td>0.9327</td><td>113.6346</td><td>27.6518</td><td>12.2918</td></tr><tr><td>27</td><td>0.0293</td><td>0.9707</td><td>138.2625</td><td>7.3109</td><td>-</td><td>64</td><td>1</td><td>-</td><td>34.3729</td><td>-</td><td>-</td></tr><tr><td>28</td><td>0.03330</td><td>0.9667</td><td>125.1738</td><td>7.4282</td><td>-</td><td>65</td><td>1</td><td>-</td><td>36.1418</td><td>-</td><td>-</td></tr><tr><td>29</td><td>0.0371</td><td>0.9629</td><td>115.1210</td><td>7.5611</td><td>-</td><td>66</td><td>1</td><td>-</td><td>38.3347</td><td>-</td><td>-</td></tr><tr><td>30</td><td>0.0409</td><td>0.9591</td><td>106.6938</td><td>7.7017</td><td>-</td><td>67</td><td>1</td><td>-</td><td>41.1288</td><td>-</td><td>-</td></tr><tr><td>31</td><td>0.0449</td><td>0.9551</td><td>99.6147</td><td>7.8511</td><td>-</td><td>68</td><td>1</td><td>-</td><td>44.8359</td><td>-</td><td>-</td></tr><tr><td>32</td><td>0.0488</td><td>0.9511</td><td>93.6306</td><td>8.0104</td><td>-</td><td>69</td><td>0.03841</td><td>0.9616</td><td>333.5709</td><td>38.6642</td><td>-</td></tr><tr><td>33</td><td>0.0529</td><td>0.9471</td><td>88.4122</td><td>8.1797</td><td>-</td><td>70</td><td>0.1713</td><td>0.8287</td><td>152.9629</td><td>38.0016</td><td>-</td></tr><tr><td>34</td><td>0.0569</td><td>0.9431</td><td>84.1192</td><td>8.3633</td><td>-</td><td>71</td><td>0.5829</td><td>0.4171</td><td>94.1997</td><td>34.8665</td><td>-</td></tr><tr><td>35</td><td>0.0603</td><td>0.9397</td><td>80.8589</td><td>8.5669</td><td>-</td><td>72</td><td>-0.1893</td><td>1.1893</td><td>311.7230</td><td>127.7285</td><td>-</td></tr><tr><td>36</td><td>0.0634</td><td>0.9366</td><td>78.2547</td><td>8.7894</td><td>-</td><td>73</td><td>-0.2028</td><td>1.2028</td><td>305.1938</td><td>128.8713</td><td>-</td></tr></tbody></table>
<table><caption>Table 2: Table of the coefficients ($\bar{\lambda}_{Z_{0,s,i}}$, $\bar{A}_{Z_{0,s,i}}$) for the neutral tungsten, as determined by the Molière's method of the Thomas-Fermi model [32], by the DHFS method from Ref. [32], and by the DFT using the GAUSSIAN code as the density of reference [12]. All methods use three exponentials.</caption><thead><tr><th>Method</th><th>Molière (Thomas-Fermi)</th><th>DHFS (3 exp.)</th><th>DFT (3 exp.)</th></tr></thead><tbody><tr><td>λ̃<sub>W,1</sub></td><td>28.4633</td><td>28.6330</td><td>40.7776</td></tr><tr><td>λ̃<sub>W,2</sub></td><td>5.6926</td><td>4.2426</td><td>4.7464</td></tr><tr><td>λ̃<sub>W,3</sub></td><td>1.4231</td><td>1.2340</td><td>1.2612</td></tr><tr><td>Ã<sub>W,1</sub></td><td>0.1</td><td>0.15</td><td>0.0964</td></tr><tr><td>Ã<sub>W,2</sub></td><td>0.55</td><td>0.6871</td><td>0.7058</td></tr><tr><td>Ã<sub>W,3</sub></td><td>0.35</td><td>0.1629</td><td>0.1978</td></tr></tbody></table>
<figure><img src="image_2.png" /><figcaption>Figure 1: Relative fraction of different screened ion charges for tungsten at different plasma electron temperatures using the *OPEN-ADAS* database [5].</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2: Deviation of the bound electron density of from spherical symmetry as estimated by the parameter $\Delta\Theta$ for all the screened charges of tungsten, from DFT calculations using GAUSSIAN code [12]. When $\Delta\Theta = 0$, the cloud of bound electrons is spherically symmetric around the nucleus. $\Delta\Theta$ is very small for all nobel gas-like electronic configurations.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: Density of bound electrons for all ionization states of tungsten calculated by the DFT method using the *GAUSSIAN* code [12]. The upper red line corresponds to neutral atom. Details of the simulations are given in the Appendix B.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: Comparison between radial densities of bound electrons for neutral tungsten $W^0$ and ionized states $W^{+10}$, $W^{+42}$ and $W^{+56}$, as calculated by GAUSSIAN (DFT method, red lines) and GRASP (MCDHF method, blue symbols) codes. An excellent agreement is found between the two quantum relativistic codes. Details of the simulations are given in the Appendix B.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Comparison between densities of bound electrons for neutral tungsten W⁰ (upper plot) and the ionized state W⁺⁴² (lower plot) as calculated by GAUSSIAN (DFT method, red line) and simple atomic models: Thomas-Fermi (blue dotted line), Yukawa with the inverse screening length λ₀,ₛ ≈ 0.9Zₛ⁰.⁴²a₀⁻¹[38] (blue dotted-dashed line) and Yukawa with the inverse screening length λ₀,ₛ⁻¹ ≈ 1.13Zₛ¹/³a₀⁻¹[24, 27, 38, 61, 62] (blue dashed line). Details of the simulation for the DFT calculation are given in the Appendix B.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6: Atomic density for neutral tungsten element, $Z_{0,s} = 0$, calculated from DFT (red full line) using GAUSSIAN code [12] and approximated using the method of moments with a single exponential (blue dashed line), two exponentials (green dotted dashed line) and three exponentials (pink dashed line) [29]. The Moliere's solution, as given in the Appendix A is also displayed (black full line) [32].</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Atomic linear density $4\pi\bar{r}^2\bar{\rho}_{Z_{0,s}}(\bar{r})$ for neutral tungsten element, $Z_{0,s} = 0$, calculated from DFT (red full line) using GAUSSIAN code [12] and approximated using the method of moments with a single exponential (blue dashed line), two exponentials (green dotted dashed line) and three exponentials (pink dashed line) [29]. The Moliere's solution, as given in the Appendix A is also displayed (black full line) [32].</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: Atomic form factor for neutral tungsten element, $Z_{0,s} = 0$, calculated from DFT (red full line) using GAUSSIAN code [12] and approximated using the method of moments with a single exponential (blue dashed line), two exponentials (green dotted dashed line) and three exponentials (pink dashed line) [29]. The Moliere's solution, as given in the Appendix A is also displayed (black full line) [32]. Only the three exponentials case is very close to the DFT solution when $\bar{q} \ge 0.08$.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9: Atomic density for charged tungsten element, Z₀,₃ = 42, calculated from DFT (red full line) using GAUSSIAN code [12] and approximated using the method of moments with a single exponential (blue dashed line), two exponentials (green dotted dashed line) and three exponentials (pink dashed line) [29].</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10: Atomic linear density $4\pi\bar{r}^2\bar{\rho}_{Z_{0,s}}(\bar{r})$ for charged tungsten element, $Z_{0,s} = 42$, calculated from DFT (red full line) using GAUSSIAN code [12] and approximated using the method of moments with a single exponential (blue dashed line), two exponentials (green dotted dashed line) and three exponentials (pink dashed line) [29].</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11: Atomic form factor for charged tungsten element, $Z_{0,s} = 42$, calculated from DFT (red full line) using GAUSSIAN code [12] and approximated using the method of moments with a single exponential (blue dashed line), two exponentials (green dotted dashed line) and three exponentials (pink dashed line) [29]. Only the three exponentials case is very close to the DFT solution when $\bar{q} \ge 0.08$.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 12: Coefficient of determination *R*² as a function of all screened ion charges *Z*₀,ₛ for tungsten to illustrate how well the multi-Yukawa reproduces the results of DFT calculated by the *GAUSSIAN* code [12]: red circles (single exponential), blue crosses (two exponentials) and green x-marks (three exponentials). When points exactly overlap in the figure, it means that despite the method is searching three exponentials, only solution with two or one exponentials are found.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 13: Linear density (lower plot) $4\pi\bar{r}^2\bar{\rho}_{Z_{0,s}}(\bar{r})$ for gold ion, $Z_{0,s} = 1$, calculated from DFT (red full line) using GAUSSIAN code [12] and approximated using the method of moments with the use of atomic linear density $4\pi\bar{r}^2\bar{\rho}_{Z_{0,s}}(\bar{r})$ for charged gold element, $Z_{0,s} = 1$, calculated from DFT (red full line) using GAUSSIAN code [12] and approximated using the method of moments with a single exponential (blue dashed line), two exponentials (green dotted dashed line), three exponentials (pink dashed line) and four exponentials (cyan full line).</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 14: Relative evolution of the square of the inverse normalized screening length as a function of the normalized ion charge. Red circles : numerical value from the multi-Yukawa description of the atomic potential with a single exponential; blue dashed line : approximate Thomas-Fermi model from Kirillov et al [22]; green dotted dashed line : $\phi_s(x) = (1 - x^{n_s+1}) / (1-x)$ from a fit of the Hartree-Fock-Slater potential [26, 28]; magenta dotted line : modified formulation $\phi_s(x) = (1 - x^{n_s+1}) / (1-x)^{3/2}$ of the fit of the Hartree-Fock-Slater (HFS) atomic potential.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 15: Normalized Fokker-Planck screening function for the neutral tungsten element, $Z_{0,s} = 0$, as a function of the normalized momentum $\bar{p} = p/(m_e c)$, calculated from DFT results (red circles) using GAUSSIAN code [12] for $\bar{\rho}_{Z_{0,s}}^{num}$ and Eq. (53) and approximated using the method of moments (multi-Yukawa) given by Eq. 50 with a single exponential (blue dashed line), two exponentials (green dotted dashed line) and three exponentials (pink dashed line).</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 16: Normalized Fokker-Planck screening function for the neutral tungsten element, $Z_{0,s} = 0$, as a function of the normalized momentum $\bar{p} = p/(m_e c)$ between $\bar{p} = [0, 1]$ calculated from DFT results (red circles) using GAUSSIAN code [12] for $\bar{\rho}_{Z_{0,s}}^{num}$ and Eq. (53) and approximated using the method of moments (multi-Yukawa) given by Eq. 50 with a single exponential (blue dashed line), two exponentials (green dotted dashed line) and three exponentials (pink dashed line). When $\bar{p} < 0.4$, the numerical value of the normalized Fokker-Planck screening function falls off more rapidly than the multi-Yukawa expression with the use of three exponentials, and does not converge towards zero for $\bar{p} \approx 0$, as expected from theory.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 17: Normalized Fokker-Planck screening function for the charged tungsten element, $Z_{0,s} = 42$, as a function of the normalized momentum $\bar{p} = p/(m_e c)$, calculated from DFT results (red circles) using GAUSSIAN code [12] for $\bar{\rho}_{Z_{0,s}}^{num}$ and Eq. (53) and approximated using the method of moments (multi-Yukawa) given by Eq. 50 with a single exponential (blue dashed line), two exponentials (green dotted dashed line) and three exponentials (pink dashed line).</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 18: Variation of the logarithm of the mean excitation energy in eV units for argon, as a function of the level of ionization. Black dotted dashed line: exponential interpolation determined from neutral atom and hydrogen-like ion according to Eq. (57); green circles: numerical results from MCSCF quantum code [14]; blue dashed line: LPA model from Eq. (59) with $\gamma_{LPA} = \sqrt{2}$ using the multi-Yukawa atomic model with the use of three exponentials calibrated against DFT calculations (GAUSSIAN code) [37]; red full line: variational quantum model from Eq. (62) using the multi-Yukawa atomic model with the use of three exponentials calibrated against DFT calculations (GAUSSIAN code) [36]; black star: mean excitation energy for the neutral atom from NIST database [60].</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure 19: Variation of the logarithm of the mean excitation energy in eV units for tungsten, as a function of the level of ionization. Black dotted dashed line: exponential interpolation determined from neutral atom and hydrogen-like ion according to Eq. (57); blue dashed line: LPA model from Eq. (59) with $\gamma_{LPA} = \sqrt{2}$ using the multi-Yukawa atomic model with the use of three exponentials calibrated against DFT calculations (GAUSSIAN code) [37]; red full line: variational quantum model from Eq. (62) using the multi-Yukawa atomic model with the use of three exponentials calibrated against DFT calculations (GAUSSIAN code) [36]; black star: mean excitation energy for the neutral atom from NIST database [60].</figcaption></figure>
## Appendix A. Molière's model

From the Molière's description of the Thomas-Fermi model [32], the atomic charge density may be expressed as $\rho_{Z_0,s}^M(\hat{r}) = (Z_s - Z_{0,s}) \sum_i \beta_i^2 B_i \exp(-\beta_i \hat{r}) / (4\pi \hat{r} b_s^3)$ with $\beta_i$ and $B_i$ independent of the atomic number $Z_s$, where the distance $r$ to the nucleus is normalized to the Thomas-Fermi atomic characteristic length $b_s$, i.e. $\hat{r} = r/b_s$. The conversion to atomic units may be easily obtained, i.e. $\beta_i \hat{r} \rightarrow \beta_i (r/a_0) (a_0/b_s) = \bar{\lambda}_{Z_{0,s,i}}^M \bar{r}$ such that $\bar{\lambda}_{Z_{0,s,i}}^M = \beta_i a_0 / b_s = \beta_i / \bar{b}_s = \beta_i Z_s^{1/3} / \alpha^{TF}$ or $\beta_i = \bar{b}_s \bar{\lambda}_{Z_{0,s,i}}^M$, where $\bar{b}_s = \alpha^{TF} Z_s^{-1/3}$ and $\alpha^{TF} = [9\pi^2/2]^{1/3} / 4 \approx 0.885$. Much in the same way, $\beta_i^2 B_i / (\hat{r} b_s^3) = (\bar{\lambda}_{Z_{0,s,i}}^M)^2 B_i / (\bar{r} a_0^3)$ and $A_{Z_{0,s,i}}^M = B_i$.

Therefore, for tungsten ($Z = 74$), since $Z_s^{1/3} = 4.1983$ while $B_1 = 0.1$, $B_2 = 0.55$, $B_3 = 0.35$, and $\beta_1 = 6.0$, $\beta_2 = 1.20$, $\beta_3 = 0.30$ from Ref. [32], coefficients which are used in plots for Figs. 6, 13 and 8 are $(\lambda_1^{W-M}, A_1^{W-M}) = (28.4633, 0.1)$, $(\lambda_2^{W-M}, A_2^{W-M}) = (5.6926, 0.55)$ and $(\lambda_3^{W-M}, A_3^{W-M}) = (1.4231, 0.35)$. Comparisons with values for a three exponentials representation of the atomic density obtained with DHFS and DFT models are given in Table 2.

## Appendix B. Quantum simulations

### Appendix B.1. GAUSSIAN code

Density functional theory (DFT) is a computational quantum mechanics modelling method used in physics, chemistry and materials science to investigate the electronic structure, principally the ground state, of many-body systems, in particular atoms. Using this theory, the properties of a many-electron system can be determined by using functionals, i.e. functions of another function, which is the spatially dependent electron density in this case.

The DFT is the method of reference in order to calculate the number density of bound electrons averaged over solid angle as a function of the radius for all ionization states of any atom. There are many tools dedicated to DFT calculations, and among them, the commercial code *GAUSSIAN* is one of reference [12]. The version g09 has been used for the calculations described in the present paper. Regarding the simplicity of the atomic configuration in a tokamak plasma, i.e. field-free atom, results obtained with this tool are likely independent of the code version. The excellent agreement with the results obtained with the latest version of the *GRASP* code described in Appendix B.2 validates this assumption.

The settings of the DFT calculations can be summarized as follows: the atomic model PBE1PBE describing the hybrid-exchange correlation functional is chosen. The basis set on which the solution is determined may be internal to the *GAUSSIAN* code (6-311G, cc-pVDZ or AUG-cc-pVDZ, AUG- standing for augmented), as for most atoms which do not require relativistic quantum calculations (up to Krypton approximately). It may be also external to the *GAUSSIAN* code, like the natural orbital-relativistic correlation consistent basis set ANO-RCC, when quantum relativistic calculations must be performed (accessible from the www.basissetexchange.org website) [63]. In the latter, the calculations are performed by solving the Douglas-Kroll-Hess second order scalar relativistic Hamiltonian for the Dirac equation, instead of solving the usual Schrödinger equation. More advanced details may be obtained from the on-line *GAUSSIAN* code documentation accessible from the gaussian.com website.

In the calculations, the spin multiplicity requires a special care. It may be obtained from Hund's rules for low-Z elements, but is usually obtained from the NIST database (physics.nist.gov/asd)
[60]. It is important to note that for six ionization states of the tungsten ranging from $W^{49+}$ to $W^{50+}$, no spin multiplicity is given likely because of the energy closeness of the different shells due to strong spin-orbit coupling. These ions require therefore a specific treatment to perform GAUSSIAN calculations.

Once GAUSSIAN calculations are carried out, results are post-processed using the Multiwfn program that can be downloaded from the sobereva.com/multiwfn website [64].

## Appendix B. GRASP code

MCDHF calculations have been done with the General Relativistic Atomic Structure Package (GRASP), version 2018 [13, 65]. The FORTRAN 95 code can be downloaded from the website github.com/compas/grasp and easily compiled. Specific scripts have been written from the documentation, using predefined ion configurations [1]: He (1s(2) = 2 electrons; [2]: Ne ([He] + 2s(2)2p(6) = 10 electrons; [3]: Ar ([Ne] + 3s(2)3p(6) = 18 electrons; [4]: Kr ([Ar] + 3d(10)4s(2)4p(6) = 36 electrons; [5]: Xe ([Kr] + 4d(10)5s(2)5p(6) = 54 electrons; [6]: Rn ([Xe] + 4f(14)5d(10)6s(2)6p(6) = 86 electrons, in order to minimize duration and memory requirements for the calculations. This is especially important for weakly ionized high-Z elements like tungsten. Once the radial wavefunctions have been calculated, the atomic density of bound electrons is determined using the dedicated module RDENSITY which can be downloaded from the CPC Library and compiled like all other modules of the GRASP code [66].

## Appendix C. Coulomb logarithm

The Coulomb logarithm $\ln \Lambda_{e,Z_{s,0}} = \ln \left( \frac{\bar{b}_{\max}}{\sqrt{1 + \bar{b}_{\min}^2}} \right)$ from Eq. (24) may be explicitly evaluated, taking into account of the plasma conditions, the type of element and its net charge. Since the Coulomb potential is screened at a distance larger than the Debye length $\lambda_D$, the upper limit is $\bar{b}_{\max} = \lambda_D/b_{90}$. For multispecies plasmas, $\lambda_D \rightarrow \lambda_D^{e-i} \approx \lambda_D^e (1 + Z_{eff}T_e/T_i)^{-1/2}$, where $\lambda_D^e$ is the usual electron Debye length, and $Z_{eff} = \sum_s \sum_{Z_{0,s}} n_{Z_{0,s}} Z_{0,s}^2 / n_e$ is the effective charge, knowing that $\sum_s \sum_{Z_{0,s}} n_{Z_{0,s}} (t, \mathbf{x}) Z_{0,s} = n_e$ from electroneutrality. It is assumed usually in most kinetic calculations that all ion species have the same temperature $T_i$ whatever their net charge.

The value of $\bar{b}_{\min}$ depends of the ratio $b_q/b_{90}$ where $b_q$ is deduced from uncertainty principle. If $b_q/b_{90} \gg 1$, quantum effects predominate and $\ln \Lambda_{e,Z_{0,s}}^q \simeq \ln(\lambda_D/b_q)$, otherwise the classical limit corresponding to $\bar{b}_{\min} = 0$ can be taken, and $\ln \Lambda_{e,Z_{0,s}}^c \simeq \ln(\lambda_D/b_{90})$ [67, 68]. Consequently, $\ln \Lambda_{e,Z_{0,s}}^q$ is less than $\ln \Lambda_{e,Z_{0,s}}^c$. The quantum limit may be determined from the uncertainty principle $\Delta p \Delta x > \hbar/2$ where the momentum increment is approximated by $\Delta p = \mu_s u_s$, $\mu_s = m_e m_s / (m_e + m_s)$ being the reduced mass between colliding particles. Therefore, assuming $\Delta x \simeq b_q$, the impact parameter is $b_q \approx \hbar / (2\mu_s u_s)$ or $b_q \approx (\lambda_C / 4\pi) (m_e / \mu_s) / \bar{u}_s$, where $\lambda_C$ is the Compton length. The ratio $\bar{b}_{\min} = b_q/b_{90} = \bar{u}_s / (2Z_{0,s}\alpha)$ since $b_{90} = r_e (Z_{0,s} / \bar{u}_s^2) (m_e / \mu_s)$ for Coulomb collisions and $\lambda_C / r_e = 2\pi / \alpha$. A rough estimate of the smooth transition between classical and quantum limits may be obtained by averaging $\bar{u}_s$ over the electron and ion distribution functions. In this case, the square root of the mean square velocity is $\langle\langle\bar{u}_s\rangle\rangle \simeq \sqrt{3\bar{T}_e}$ where $\bar{T}_e = T_e / (m_e c^2)$ and the normalized plasma temperature threshold above which quantum effects are significant is $\bar{T}_{e,Z_{0,s}}^q = 4Z_{0,s}\alpha^2 / 3$. For $Z_{0,s} = 1$, $T_{e,Z_{0,s}}^q = 36$ eV, such that the quantum limit must be always taken in standard tokamak plasma conditions with isotopes of hydrogen.
For $Z_{0,s} = 42$, corresponding to the net ionization of tungsten at $T_e = 3$ keV, then $T_{e,Z_{0,s}}^q$ is much larger, $T_{e,Z_{0,s}}^q = 64$ keV, and the classical limit is conversely always valid in tokamak plasmas. Since $m_e \ll m_s$, $\mu_s \simeq m_e$, and $\ln \Lambda_{e,Z_{0,s}}^c \simeq \ln \lambda_D/r_e + 2 \ln \bar{u}_s - \ln Z_{0,s}$, while $\ln \Lambda_{e,Z_{0,s}}^q \simeq \ln \lambda_D/r_e + \ln \bar{u}_s + \ln (2\alpha)$. If $\ln \Lambda_{e,Z_{0,s}}^q$ and $\ln \Lambda_{e,Z_{0,s}}^c$ are both heavily weighted by the ratio $\lambda_D/r_e$, the regime dominated by quantum effects concerns principally fast electrons for partially ionized high-Z elements whose kinetic energy is greater than $\bar{T}_{e,Z_{0,s}}^q$. In the quantum limit, $\Lambda_{e,Z_{0,s}}^a$ is independent of the ion net charge $Z_{0,s}$.

In standard MKSA units, with $\lambda_D = \lambda_D^e$, $\ln \Lambda_{e,Z_{s,0}}^q = 0.5 \ln T_e [\text{keV}] - 0.5 \ln n_e [10^{20} \text{m}^{-3}] + \ln \bar{u}_s + 18.61$ and the thermal value is $\ln \Lambda_{e,Z_{0,s}}^{q-th} \simeq \ln T_e [\text{keV}] - 0.5 \ln n_e [10^{20} \text{m}^{-3}] + 16.04$. Using $\lambda_D^{e-i}$, the Coulomb logarithm must be reduced by the term $-0.5 \ln (1 + Z_{eff}T_e/T_i)$, and for a pure hydrogen plasma with $T_e = T_i$, $\ln \Lambda_{e,Z_{0,s}}^{q-th} \simeq \ln T_e [\text{keV}] - 0.5 \ln n_e [10^{20} \text{m}^{-3}] + 15.7$, a value very close to those found in the literature [67, 69]. Additional small differences may arises from the choice of the averaged velocity. For relativistic electrons, $\bar{u}_s \simeq \bar{v}$, since ions may be considered at rest, $\ln \Lambda_{e,Z_{0,s}}^{q-rel} \simeq \ln \Lambda_{e,Z_{0,s}}^{q-th} + \ln \bar{p} - 0.5 \ln \bar{T}_e - 0.5 \ln 3$, as far as $\bar{p} > \sqrt{3\bar{T}_e}$ or $\bar{E}_c \gg 3\bar{T}_e/2$, where $\bar{E}_c$ is the kinetic energy normalized to the electron rest mass energy. Following Ref. [7], $\ln \Lambda_{e,Z_{0,s}}^q$ may be approximated by $\ln \Lambda_{e,Z_{0,s}}^q \simeq \ln \Lambda_{e,Z_{0,s}}^{q-th} + \ln \left(1 + \left(\bar{p}/(\sqrt{3\bar{T}_e})\right)^k\right)/k$ with $k=5$, in order to have a smooth transition from the thermal limit of the Coulomb logarithm.

Conversely to the quantum limit, $\ln \Lambda_{e,Z_{0,s}}^c$ in the classical limit is weakly dependent of the ion charge $Z_{0,s}$ and is more sensitive to $T_e$, since $\ln (\Lambda_{e,Z_{0,s}}^c/\Lambda_{e,Z_{0,s}}^q) = \ln (\bar{u}_s/(2\alpha Z_{0,s}))$. In MKSA units, $\ln \Lambda_{e,Z_{0,s}}^{c-th} = 1.5 \ln T_e [\text{keV}] - 0.5 \ln n_e [10^{20} \text{m}^{-3}] - \ln Z_{s,0} + 17.69$ and $\ln \Lambda_{e,Z_{0,s}}^{c-rel} = \ln \Lambda_{e,Z_{0,s}}^{c-th} + 2 \ln \bar{p} - \ln \bar{T}_e - \ln 3$. It can be also approximated by the expression $\ln \Lambda_{e,Z_{0,s}}^c \simeq \ln \Lambda_{e,Z_{0,s}}^{c-th} + \ln (1 + (\bar{p}^2/(3\bar{T}_e))^k)/k$. The choice of the electron-ion Coulomb logarithm, $\ln \Lambda_{e,Z_{0,s}}(\bar{p}, \bar{T}_e)$, depends therefore of the type of elements $s$, its local ionization state $Z_{0,s}$ and the temperatures $T_e$ and $T_i$ at the same location in the plasma.
References

[1] Bourdelle C. et al 2015 *Nucl. Fusion* 55 063017

[2] Garofalo A. M. et al 2017 *Nucl. Fusion* 57 076037

[3] Gribov Y. 2006 ftp.jp.iter.org, Available online

[4] Peysson Y et al 2021 IAEA FEC 2020 cea-03249410

[5] Summers H. P. 2004 *The ADAS user manual, version 2.6*

[6] Hesslow L. et al 2017 *Phys. Rev. Letter* 118 255001

[7] Hesslow L. et al 2018 *J. Plasma Phys.* 84(6) 905840605

[8] Hesslow L. et al 2018 *Plasma Phys. Control. Fusion* 60 074010

[9] Stahl A. et al 2016 *Nucl. Fusion* 56 112009

[10] Peysson Y and Decker J. 2014 *Fusion Science and Technology* 65 22

[11] Ostuni V. et al 2022 *Nucl. Fusion* 62 106034

[12] Frisch M. J. et al 2016 *Gaussian 09, revision e.01, Technical report, Gaussian, Inc.*, Wallingford
*CT*

[13] Froese Fischer C. et al 2019 *Computer Physics Communications* 237 184

[14] Sauer S. P. A. 2015 *Advances in Quantum Chemistry: The Mean Excitation Energy of Atomic Ion, Elsevier Inc.* 71(3) 29

[15] Sauer S. P. A. et al 2020 *Molecular Physics* 119(5)

[16] Belkhiri M. et al 2015 *Phys. Rev. A* 92 032501

[17] Gu. M. F. 2008 *Can. J. Phys.* 86 675

[18] Koch H. W. and Motz. J. W. 1959 *Rev. Mod. Phys.* 31(4) 920

[19] Peysson Y and Decker J. 2008 *Phys. Plasmas* 15(9) 092509

[20] Amore P. et al 2014 *Applied Mathematics and Computation* 232 929

[21] Green A. E. S. et al 1969 *Phys. Rev.* 184(1) 1

[22] Kirillov V. D. et al 1975 *Fizika Plazmy* 1 218

[23] Edwards J. P. et al 2017 *Prog. Theor. Exp. Phys.* 108A01

[24] Pratt R. H. and Tseng H. K. 1972 *Phys. rev. A* 5(3) 1063

[25] Pratt R. H. et al 1977 *Atomic Data and Nuclear Data Tables* 2e 175

[26] Avdonina N. B. and Pratt R. H. 1993 *J. Quant. Spectrosc. Radiat. Transfer* 50(4) 349

[27] Pratt R. H. et al 1995 *Nuclear Instruments and Methods in Physics Research* B99 156

[28] Lamoureux M. and Avdonina N. 1997 *Phys. Rev. E* 55(1) 912

[29] Salvat F. et al 1987 *Phys. Rev. A* 36(2) 467

[30] Fronsdal C. and Uberall H. 1958 *Phys. Rev.* 111(2) 580

[31] Haug E. 2008 *Radiation Physics and Chemistry* 70 207

[32] Moliere G. 1947 *Zeitschrift Naturforschung* 2a 133

[33] Nilsson E. et al 2015 *Plasma Phys. Control. Fusion* 57(9) 095006

[34] H. Bethe H. 1930 *Annalen der Physik*, 325
[35] Jackson J. 1998 *Classical Electrodynamic*, John Wiley & Sons, Inc

[36] Garbet X. et al 1987 *Journal of Applied Physics* 61 907

[37] Linhard J. and Scharff M. 1953 *Dan. Mat. Fys. Medd* 27(15)

[38] Botto D. J. et al 1978 *Phys. rev. A* 18(2) 580

[39] Sauer S. P. A. 2011 *Molecular Electromagnetism. A Computational Chemistry Approach*, Oxford University Press

[40] Kohn C. and Ebert U. 2014 *Atmospheric Research* 135-136 432

[41] Walkowiak J. et al 2022 *Phys. Plasmas* (29) 022501

[42] Karney C.F.F. 1986 *Comp. Phys. Rep.* 4 183

[43] Shoucri M. and Shkarofsky I. 1994 *Comp. Phys. Comm.* 82 287

[44] Cannoni M. 2017 *International Journal of Modern Physics A* 32 1730002

[45] Decker J. et al 2016 *Plasma Phys. Control. Fusion* 58 025016

[46] Heitler W. 1957 *Clarendon Press, Oxford, UK*, 3rd edition

[47] Solodov A. A. and Betti R. 2008 *Phys. Plasmas* 15 042707

[48] Decker J. and Peysson Y. 2004 *EUR-CEA-FC-1736, Euratom-CEA*

[49] Beliaev S.T. and Budker G.I. 1956 *Sov. Phys. Dokl.* 1 218

[50] Bethe H. and Heitler W. 1934 *Proceedings of the Royal Society of London A* 146 83

[51] Abramowitz M. and Stegun I. A. 1970 *Handbook of mathematical functions*, Dover Publication, New York, USA, 9th edition

[52] Fano U. and Cooper J. W. 1968 *Rev. Mod. Phys.* 40(3) 441

[53] Chu W.K. and Powers D. 1972 *Phys. Letters* A 40A(1) 23

[54] Inokuti M. and Turner J. E. 1978 *Physics of Elementary Particles and Fields (A3200)* CONF-780534-2

[55] Wyckoff H. O. 1984 *Technical report, International Commission on Radiation Units and Measurements (ICRU) REPORT 37*

[56] Kamakura S. et al 2006 *J. Appl. Phys.* 100 064905

[57] Sauer S. P. A. et al 2018 *J. Chem. Phys.* 148 174307

[58] Rosendorff S. and Schlaile H. G. 1989 *Phys. Rev. A* 40(12) 6892

[59] Dehmer J. L. et al 1975 *Phys. Rev. A* 12(1) 102

[60] Kramida A. et al 2022 *Nist atomic spectra database (ver. 5.10)*, [online], National Institute of Standards and Technology

[61] Hyatt W. J. 1956 *Phys. Rev.* 104(5) 1298

[62] Bunaciu D. et al 1980, *Nuclear Physics A* 339 329

[63] Pritchard B. P. et al 2019 *J. Chem. Inf. Model.* 59(11) 4814

[64] Lu T. and Chen F. 2012 *J. Comput. Chem.* 33 580

[65] Jonsson P. et al 2023 *Atoms* 11(7) atoms11010007

[66] Schiffmann S. et al 2022, *Computer Physics Communications* 278 108403

[67] Wesson J. 2004 *Tokamaks, The International Series of Monographs on Physics*. Clarendon Press
[68] Mulser P. *et al* 2014, *Phys. Plasmas* 21 042103

[69] Sauter O. *et al* 1999 *Phys. Plasmas* 6(7) 2834