A unified description of atomic physics for electron Fokker–Planck calculations

Most realistic kinetic calculations for tokamak plasmas are now required to incorporate the effect of partially ionized high-Z elements arising either from uncontrolled influxes of metallic impurities, such as 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 therefore be modified, because all plasma atoms are not entirely ionized, as is the case for light elements. This represents a challenge, in order to perform fast but also accurate calculations, regardless of the type of element present in the plasma, but also their local levels of ionization while covering a wide range of electron energies in a consistent way, from a few keV to tens of MeV in plasmas whose electron temperature may itself vary from 10 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 one to derive analytical formulations for both elastic and inelastic scattering, which can then be easily incorporated in kinetic calculations. The impact of the number of exponentials in the description of the atomic potential is discussed, and a comparison with simple and advanced atomic models is also performed.


Introduction
The use of tungsten (W) as the plasma-facing material in present-day experimental fusion devices, such as WEST [1], Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
EAST [2] and the International Thermonuclear Experimental Reactor (ITER) that is currently being constructed [3], has raised the question of the impact of partially ionized high-Z impurities on the performance of hot plasmas.For example, the ability to efficiently drive the toroidal plasma current by radio-frequency (RF) electromagnetic waves may be reduced by an enhanced electron pitch-angle scattering and electronion slowing-down due to tungsten, thus limiting the capability of control for improved plasma performance [4].In standard tokamak fusion plasmas, the usual electron temperature T e is supposed to be in the range between 1 and 10 keV, so Relative fraction of different screened ion charges for tungsten at different plasma electron temperatures using the OPEN-ADAS database [5].
that most of the 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 Z 0,s ≃ 42 in a plasma whose electron temperature is T e = 3 keV according to the OPEN-ADAS database, as shown in figure 1, so that N s ≃ 32 electrons are still bound [5].Even at T e = 10 keV, as expected in ITER plasmas, N s is still large for the tungsten element, in the order of twenty 4 .
The role played by the screening of partially ionized high-Z elements has first been considered in order to accurately describe the dynamics of runaway electrons in very cold postdisruptive plasmas, and 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 compared to the traditional picture because of the partial screening, with a significant impact on the critical electric field (Hastie-Connor) beyond which electrons may run away [6][7][8].This original work, implemented in the CODE code dedicated to runaway electron physics in almost zero-temperature postdisruptive plasmas [9], was later extended to standard tokamak regimes in the LUKE solver of the 3D linearized bounceaveraged relativistic electron Fokker-Planck equation [10].This allows us to describe the consequences of uncontrolled impurity influxes of high-Z metallic elements, such as tungsten 4 The mean screened ion charge is defined as Z 0,s = ∑ 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, ∑ i f 0,s,i = 1, Z 0,s,i ∈ {0, Zs}.The number of bound electrons is Ns,i = Zs − Z 0,s,i and the mean value is Ns = ∑ i f 0,s,i Ns,i.
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 that must be described in the 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 nonrelativistic codes describing the ground-state and mean excitation energies for different ionization states [12][13][14][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 a consequence of the relativistic effects and the resulting complex orbital 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 γ s significantly exceeds unity, where γ 2 s = (αZ s ) 2 + 1, and α is the usual fine structure constant.While for argon, relativistic corrections are negligible since γ Ar s ≃ 1.0086, they become more significant for tungsten, as γ W s ≃ 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 5 implemented in the GAUSSIAN and GRASP codes, respectively [12,13], 6 the mean excitation energies for all ionization states that play an important role in the inelastic electron-ion scattering processes are still not available.In much 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 continuously cover the range of kinetic energies from a few keV to several tens of 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 a considerable numerical task.
In this context, the accurate incorporation of 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 bremsstrahlung, a major moment of the non-thermal distribution function for diagnosing fast electron dynamics.The use of analytical formulas based on simple parameterized models is consequently a more suitable approach, with absolute calibration against results obtained by advanced numerical quantum codes.However, simple atomic models usually have limited applicability, which prevents the systematic and consistent use of all quantities that must be modified to take into account 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 usually give less accurate results when the number N s of bound electrons is small compared to Z s .The charge density of the inner shell is usually better described by a Yukawa electrostatic potential [23], which itself is not relevant for neutral or weakly ionized atoms.This is because the charge density fall-off is generally too sharp at large distances from the nucleus 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 usually involves deep electronic shells to calculate the radiation emission [24][25][26][27][28]. Consequently, a unified and accurate description of the atomic electrostatic potential that can be used either for kinetic calculations or for bremsstrahlung without a significant degradation of the numerical performance of the kinetic code, whatever the plasma conditions (cold or hot) and the type of element, is of great interest.
The purpose of this paper is therefore to propose a general and global approach for incorporating atomic physics in Fokker-Planck electron solvers, allowing existing codes to be easily and robustly updated for realistic simulations, whatever the consequences on 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.There is a consequent trade-off between keeping codes fast and accurate over a wide range of electron kinetic energies, regardless of the ionization state of the elements, and with a reasonably correct description of the physics involving both outer and inner electron shells.This approach has already been successfully considered for modeling results of Hartree-Fock-Slater calculations for neutral elements only, from hydrogen to uranium, using up to three exponentials [29].In addition, it was successfully applied for deriving an analytical formula for the bremsstrahlung with screening effects, valid from the classical to the fully relativistic limits [30,31].The modeling here considered, the so-called multi-Yukawa (MY), is basically an extension of Moliere's approach initially used to describe 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 due 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 section 2, the MY atomic model is introduced, and the calibration procedure is explained in detail for an arbitrary number of exponentials.Comparison with DFT and 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 the presence of partially ionized high-Z elements, inelastic electron-ion collisions must also be taken into account, since free electrons in 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 the runaway electron population in postdisruptive plasmas [33].The incorporation of atomic physics in the Fokker-Planck collision operator is first detailed in section 3, for both elastic and inelastic processes.The screening function describing the impact of partially ionized atoms on elastic Coulomb collisions is derived in section 4, using the MY atomic model.The inelastic electron-ion collisions are then considered in section 5, with the approximate approach based on Bethe's formula for electron energy losses per unit length [34,35].In this case, the atomic physics is described by the mean excitation energy of the ground state, which can be calculated by several methods.This quantity can be derived from a non-relativistic variational quantum approach, but also with the classical local plasma approximation (LPA), both using the MY atomic model [36,37].Conclusions are given in section 6.
Although the impact of the screening on 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.

Radial distribution of charge in the ground state and form factor
Kinetic calculations with partially ionized high-Z atoms require an atomic model that accurately describes the spatial distribution of bound electrons ρ Z0,s (r) in the ground state, regardless of the type of atom and its level of ionization, where 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 7 .By definition, ´V ρ Z0,s dr = N s = Z s − Z 0,s where V is the 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 by evaluating the form factor F Z0,s (q) ≡ ´V exp (−iq • r/h) ρ Z0,s dr, where q is the usual recoil momentum in Coulomb collisions.Since the kinetic energy is conserved in the elastic scattering process and assuming that small-angle scattering predominates, |q| ≃ 2 |p| sin (θ/2), with p being the incoming electron momentum and θ the deflection angle.Here, q = p f − p i , where V of the incoming |i⟩ and outgoing electron |f ⟩, respectively, both being considered as plane waves (first Born approximation) 8 .Based on Fermi's golden rule, the relativistic Mott crosssection that describes Coulomb collisions in kinetic calculations must be modified according to the simple rule Z s → Z s − F Z0,s (q) in order to account for the partial atomic screening.For low-energy electrons, since lim ∥q∥→0 F Z0,s (q) = N s , the ion is fully screened, i.e.Z s → Z 0,s = Z s − N s , while conversely, for very energetic electrons, it is fully stripped, since lim ∥q∥→∞ F Z0,s (q) = 0. 9 Within this framework, the form factor F Z0,s (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.

Description of the MY 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 Z0,s (q) = N s / 1 + q a Z0,s /2 m , where āZ0,s ≡ 2a Z0,s /α following the notation used in [7], and q = q/ (m e c) with q = ∥q∥.Here, a Z0,s may be considered as an effective radius of the ion of charge Z 0,s which depends on the chosen atomic model, α is the fine structure constant, c is the speed of light and m e is the electron rest mass.While m = 3/2 with a Z0,s = 3N 2/3 s / (4Z s ) for the approximate Thomas-Fermi model [22], m = 2 and a Z0,s = λ −1 Z0,s for the Yukawa electrostatic potential, where λ −1 Z0,s 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][25][26][27][28]. 10  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 of the rank-two tensor, directly obtained from DFT calculations [39], which essentially measures the deviation of the charge distribution ρ Z0,s (r) from spherical 8 The volume V is chosen so that ⟨i | i⟩ = ⟨f | f⟩ = 1, a condition necessary to have a probabilistic interpretation of the wave functions. 9The procedure is general and may be applied to all cross-sections derived within the first Born approximation as for bremsstrahlung. 10For a neutral atom of atomic number Zs and if the Thomas-Fermi model is used, s a −1 0 , a value frequently found in the literature [38].Here, bs is the atomic radius in the Thomas-Fermi model.Deviation of the bound electron density of from spherical symmetry as estimated by the parameter ∆Θ for all the screened charges of tungsten, from DFT calculations using the GAUSSIAN code [12].When ∆Θ = 0, the cloud of bound electrons is spherically symmetric around the nucleus.∆Θ is very small for all nobel gas-like electronic configurations.symmetry.It is evaluated by a global parameter, ∆Θ = |(max (XX, YY, ZZ) − min (XX, YY, ZZ)) / max (XX, YY, ZZ)|, where the diagonal elements of the tensor are XX, YY and ZZ.According to this simple definition, ∆Θ = 0 corresponds to a perfect spherical symmetry for which all diagonal elements are identical.It turns out that this parameter progressively increases with the ionization level, as shown in figure 2. It is always much lower than 0.15 for Z 0,s ⩽ 40, and very small for all noble gas-like electronic configurations, regardless of the Z 0,s value.Above Z 0,s = 40, some electronic configurations exhibit larger departure from spherical symmetry, but they concern primarily a few values for Z 0,s ⩾ 56, which will hardly ever be found in tokamak plasmas.
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 only be obtained with m = 2, as demonstrated for the case of a neutral atom [30,40].Therefore, with the constraint of performing fast and accurate kinetics, but also bremsstrahlung calculations, based on analytical formulas within 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 of neutral or weakly ionized atoms, the simplest approach is to consider, instead, a generalized Yukawa potential, here named MY, which can be expressed as a series of exponentials: each of which accurately describes the charge distribution around the nucleus, either close to or far from it.From the Poisson's equation △U Z0,s = ∇ 2 U Z0,s = −ρ Z0,s /ϵ 0 , the radial normalized distribution of bound electrons is, where the density ρ Z0,s (r) = ρ Z0,s (r) a 3 0 is in atomic units, λ Z0,s,i ≡ λ Z0,s,i a 0 and A Z0,s,i = A Z0,s,i / (Z s − Z 0,s ), λ Z0,s,i being the inverse of the normalized characteristic length and A Z0,s,i the weight of the i th Yukawa exponential, respectively.Here, r ≡ r/a 0 , where a 0 is the classical Bohr radius.The corresponding form factor is therefore, where a Z0,s,i ≡ 2λ This approach has long been considered appropriate to describe the Thomas-Fermi atomic potential by Moliere using three exponentials [32].The correspondence between coefficients (B i , β i ) found in the literature and A Z0,s,i , λ Z0,s,i is given in appendix A. This method has also been used to fit the density of bound electrons calculated by a DHFS code for neutral atoms only, whose Z s value ranges from 1 (hydrogen) to 92 (uranium) [29].For most elements above argon, approximately three exponentials are necessary to accurately reproduce the radial distribution of charges when ionization is weak.Naturally, the analytical density given by equation (2) can only partially reproduce the oscillations of the DHFS linear density 4πr 2 ρ Z0,s (r) associated with different inner shell contributions.However, the approximate form factor F Z0,s (q) given by equation (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.

Absolute model calibration
In the present study, the method used in [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 of finding a full set of positive λ Z0,s,i values.For some elements with an atomic number larger than tungsten, such as gold (Z s = 79), up to four exponentials can be found using the calibration procedure, but for lower Z s values, the number of exponentials usually never exceeds three, as found for tungsten.The determination of A Z0,s,i and λ Z0,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 under 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 those derived from 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 A Z0,s,i , λ Z0,s,i of the MY description are determined from the condition ⟨r l ⟩ = ⟨r l ⟩ num , with r l num being the moment of order l calculated numerically from the density of bound electrons ρ num Z0,s (r) obtained using advanced atomic quantum codes.Here, from the MY density given by equation (2), where Γ (z) is the Gamma function.Defining R l ≡ ⟨r l ⟩/ (l + 1)!, a set of 2l equations depending upon the number of parameters A Z0,s,i , λ Z0,s,i to be determined is obtained, and A Z0,s,i , λ Z0,s,i are calculated by solving the equation R l = R num l , where the number l being an integer greater than −1.
For a fit with a single Yukawa potential, only two terms remain, and since i A Z0,s,i = 1, it can be deduced that A Z0,s,1 = 1 and λ Z0,s,1 = R N −1 .For two exponentials, four equations with four unknowns must be considered.By grouping the equations, and the values λ are therefore solutions of the quadratic equation, If both roots are real and positive, they correspond to λ Z0,s,1 , λ Z0,s,2 , respectively, so that and A Z0,s,2 = 1 − A Z0,s,1 , otherwise, a single exponential must be considered for the modelization.In this case, there is some loss of accuracy in the modeling of the bound electron density, but it has a moderate impact on the form factor, regarding its definition given by equation ( 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.This arises principally for highly ionized atoms because the highest moments R num l are too small compared to the lowest ones.The weight of ρ N Z0,s (r) at large r is therefore unable 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 λ values must be real and positive solutions of the polynomial equation to ensure that the atomic potential may be well described by a MY potential with the use of three exponentials.Otherwise, two exponentials must be considered in the modeling procedure, thus removing useless moments associated with three exponentials.Defining the vector the coefficients A Z0,s,i are determined from the matrix relation A = N −1 MX, where parameters A Z0,s,i are components of the vector, with and Formally, it is possible to easily extend this method by recurrence to any number of exponentials.The dimensions of X, N and M must be adjusted according to the number of exponentials, as well as the degree of the polynomial equation in λ to be solved, its coefficients being determined by expanding the product i λ − λ Z0,s,i .In the numerical implementation of the method of moments, the possibility of finding up to four exponentials has been considered.However, a larger number of exponentials is rarely found, only for a few low ionization states of elements heavier than tungsten, such as gold.For tungsten, the maximum number of exponentials never exceeds three, regardless of its ionization state.As discussed previously, if no solution is found for a given set of multiple exponentials, a solution is sought for a number of exponentials decremented by one unity, and the procedure is repeated until a set of positive and real λ values is found.The case with a single exponential corresponding to the standard Yukawa potential is the ultimate solution if a MY potential cannot be found numerically.

Comparisons between MY 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 of 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 of 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 figure 3, the radial profiles of the density of bound electrons for all ionization states of tungsten calculated using the GAUSSIAN code are displayed.For low ionization states, the density clearly exhibits 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 tungsten atoms and the ions W 10+ , W 42+ and W 56+ in figure 4. Consequently, numerical densities of reference ρ N Z0,s (r) given by one or the other of the two codes can be used indifferently for determining the coefficients A Z0,s,i , λ Z0,s,i of the MY description.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 the GAUSSIAN (DFT method, red lines) and GRASP (MCDHF method, blue symbols) codes.Excellent agreement in found between the two quantum relativistic codes.Details of the simulations are given in appendix B.
In the literature, a comparison between DFT and the simple atomic models frequently used in publications has been carried out.Here, the radial dependencies of the density of bound electrons are evaluated for neutral tungsten and the ion W 42+ using the Thomas-Fermi model as well as the standard Yukawa one (single exponential).For the latter, two inverse screening lengths have been considered: λ B 0,s ≃ 0.9Z 0.42 s a −1 0 from a fit of the Herman-Skillman potentials determined by solving the DHFS equations [38] and , where b s is the reference length in the Thomas-Fermi model.As shown in figure 5, a good quantitative agreement is observed between DFT and the Thomas-Fermi model for neutral tungsten, as expected from the theory, while the agreement is poor with the standard Yukawa model, regardless of the inverse screening length.Conversely, the agreement for W 42+ between DFT and the 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 applications for accurately describing the atomic physics in kinetic and radiation calculations since the types of elements in the plasma may change with operating conditions, while their ionization states can also vary considerably with the temperature of the plasma.
Using results obtained with the DFT model, the set of coefficients A Z0,s,i , λ Z0,s,i has been determined for up to four exponentials, but for elements lighter than gold, the maximum number of exponentials never exceeds three, as for tungsten.The full list of all ionization states of tungsten is given as a reference in table 1.For very weakly ionized states, Z 0,s ⩽ 5, three exponentials are found by the numerical procedure, because of the different slopes in the radial density, as shown in figure 3, while the number of exponentials is usually lower for larger Z 0,s values since the decrease in the radial density from the nucleus becomes 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 figure 6.
With the use of three exponentials, an excellent quantitative agreement is found between ρ N Z0,s (r) determined by DFT and the approximate MY description ρ Z0,s (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 Molière's description of the Thomas-Fermi model Table 1.MY coefficients (3 exponentials) for the different ionization states of tungsten, based on DFT calculations done with the GAUSSIAN code for the density of [12].Note that A W,3 = 1 − A W,2 − A W,1 by definition and in the case of two exponentials, Atomic density for neutral tungsten element, Z 0,s = 0, calculated from DFT (red full line) using the 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].
Moliere's solution, as given in appendix A is also displayed (black full line) [32].
is also in very good agreement with the DFT results, which is consistent with the results shown in figure 5.Even if the radial dependence of the density of the nucleus is well reproduced, in the core of the atom, oscillations of the linear density 4πr 2 ρ num Z0,s (r) cannot be well reproduced by a series of exponentials, as shown in figure 7.This intrinsic limitation has already been observed in [29], where the same approach is considered but with another reference atomic model (DHFS).However, since the discrepancy occurs for r < 0.5, its impact on the Atomic linear density 4πr 2 ρ Z0,s (r) for neutral tungsten element, Z 0,s = 0, calculated from DFT (red full line) using the 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].Moliere's solution, as given in appendix A is also displayed (black full line) [32].
form factor remains small, as displayed in figure 8, when the normalized recoil momentum q is less than 0.1.Knowing that most of the Coulomb collisions occur principally for very low q values corresponding to the first Born approximation, the MY description is therefore remarkably robust for very weakly ionized atoms.
For higher ionization states, the advantage of the MY description is that it remains accurate for describing both the density and the form factor.An example is given in figure 9 for Atomic form factor for neutral tungsten element, Z 0,s = 0, calculated from DFT (red full line) using the 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].Moliere's solution, as given in appendix A is also displayed (black full line) [32].Only the three exponentials case is very close to the DFT solution when q ⩾ 0.08.Atomic density for charged tungsten element, Z 0,s = 42, calculated from DFT (red full line) using the 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].
W 42+ .In this case, the agreement is very good for both two and three exponentials, while the standard Yukawa description with a single exponential has rather poor agreement with the radial dependence of the density determined by DFT.However, as expected, the discrepancy is less pronounced compared to the case of neutral atoms due to the small remaining bumps in the radial density.For the neutral atom case, the small oscillations of the numerical linear density 4πr 2 ρ num Z0,s (r) are not well reproduced by the series of exponentials, as shown Figure 10.Atomic linear density 4πr 2 ρ Z0,s (r) for charged tungsten element, Z 0,s = 42, calculated from DFT (red full line) using the 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].Atomic form factor for charged tungsten element, Z 0,s = 42, calculated from DFT (red full line) using the 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 q ⩾ 0.08. in figure 10, but the departure from ρ num Z0,s (r) has again a very small impact on the form factor (see figure 11).
In table 2, the set of values A Z0,s,i , λ Z0,s,i obtained for the neutral tungsten from the DFT and DHFS methods are given for comparison [29].The coefficients of Moliére's method are also reported [32].Even if the methodology is similar to the one detailed in [29], the coefficients for the three exponentials case exhibit a quite significant difference.Nevertheless, Table 2.The coefficients λ Z0,s,i , A Z0,s,i for the neutral tungsten, as determined by Molière's method of the Thomas-Fermi model [32], by the DHFS method from [32] and by DFT using the GAUSSIAN code as the density of [12].All methods use three exponentials.
Method Molière (Thomas-Fermi) DHFS ( Coefficient of determination R 2 as a function of all screened ion charges Z 0,s for tungsten to illustrate how well MY 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, this means that despite the method seeking three exponentials, only solutions with two or one exponentials are found.
as shown in figure 13, the approximate MY linear densities remain fairly close to the values obtained by DFT, even if not all the oscillations can be well reproduced.The differences between the coefficients result from their large sensitivity due to small changes, illustrating the ill-conditioned nature of the problem here addressed.This justifies a posteriori the chosen method, compared to a standard least-squares fit procedure that cannot converge when multiple close solutions exist [43].
In figure 12, the usual coefficient of determination R 2 is displayed for all the ionization states of tungsten to illustrate how well the MY reproduces the results of DFT for the atomic density 11 .The fact that R 2 > 0.44 regardless of the number of exponentials indicates that the MY is an appropriate simplified 11 The coefficient of determination is calculated according to the standard formula R 2 = 1 − SSres/SStot where SSres = ∑ j (ρ num Z0,s (rj) − ρ Z0,s (rj)) 2 is the residual sum of squares and SStot = ∑ j (ρ num Z0,s (rj) − ⟨ρ num Z0,s ⟩) 2 is the total sum of squares.The sum is performed on all radial locations rj.Here, ⟨ρ num Z0,s ⟩ = Figure 13.Linear density (lower plot) 4πr 2 ρ Z0,s (r) for gold ion, Z 0,s = 1, calculated from DFT (red full line) using the GAUSSIAN code [12] and approximated using the method of moments with the use of atomic linear density 4πr 2 ρ Z0,s (r) for charged gold element, Z 0,s = 1, calculated from DFT (red full line) using the 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).
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.This also illustrates the limits of this model for very weakly ionized atoms, which justifies the need for the MY description.
Conversely, for up to three exponentials, R 2 always remains higher than 0.8, and often close to unity even for neutral or weakly ionized atoms.Between Z 0,s = 45 and Z 0,s = 55, as well as 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 ρ Z0,s (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 figure 12, and the density of bound electrons is therefore less accurately described compared to other ionization states.Nevertheless, detailed calculations have shown that the impact on the atomic form factor still remains moderate.The fact that R 2 ⩾ 0.7 indicates that the MY model remains well consistent with the results of DFT and the MCDHF codes, even for these more difficult cases.
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 figure 13.Some differences can be seen in the inner part of the bound electron density between three and four exponentials, especially when is the mean of the calculated bound electron density.With this definition, 0 ⩽ R 2 ⩽ 1 if the model is consistent with numerical data.Zero indicates a very poor agreement, and if R 2 = 1, the agreement is perfect.

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 MY 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: φs(x) = (1 − x ns+1 )/(1 − x) from a fit of the Hartree-Fock-Slater potential [26,28]; magenta dotted line: modified formulation φs(x) = (1 − x ns+1 )/(1 − x) 3/2 of the fit of the Hartree-Fock-Slater (HFS) atomic potential.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 MY solution with the highest possible number of exponentials gives a better agreement with DFT and MCDHF calculations.
Although 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 λ 2 Z0,s /λ 2 0,s = φ 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], φ s (x) = (1 − x) −4/3 , while from a fit of DHFS calculations, a heuristic dependency of the form φ s (x) = 1 − x ns+1 / (1 − x) was guessed with n s = Z s (1/3 − 0.0020 × Z s ) [26,28].As shown in figure 14, λ 2 Z0,s /λ 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 [26,28] is far from the numerical calculations based on DFT, it is found that a modified formulation φ s gives a much better agreement.This improvement is valid for all elements, whatever Z s is, 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 λ 2 Z0,s /λ 2 0,s with Z 0,s /Z s can be reasonably well reproduced.

Elastic scattering
In kinetic calculations, the Coulomb collision operator may be expressed formally as df e /dt| coll ≡ s Z0,s C e,Z0,s f e , f Z0,s + C e,e (f e , f e ) where C e,Z0,s f e , f Z0,s describes the interactions between the momentum distribution function f e (t, x, p) of test electrons and the momentum distribution function f Z0,s (t, x, p s ) of atoms of species s with an ionization state Z 0,s , while C e,e (f e , f e ) is the linearized electronelectron collision operator 12  [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 Z0,s (t, x) = ´fZ0,s (t, x, p s ) d 3 p s of ions with a net charge Z 0,s at the spatial location x results from the local balance between ionization and recombination processes, assuming in general that f Z0,s is a Maxwellian distribution.The relative fraction of partially ionized atoms is given by the ratio n Z0,s (t, x) /n s (t, x), where n s (t, x) = Z0,s n Z0,s (t, x) may be obtained by considering a local collisional-radiative equilibrium, as for the OPEN-ADAS database in the LUKE code [5,10].
The incorporation of the partial screening effects for elastic scattering in kinetic calculations requires that we re-express the friction vector and diffusion tensor A Z0,s and D Z0,s of the Fokker-Planck formulation of the collision operator, assuming that small angle scattering still predominates for Coulomb collisions, which remains a good assumption even in the presence of high-Z impurities [7], C e,Z0,s f e , f Z0,s ≃ −∇p A Z0,s fe (t, x, p) + ∇p∇p D Z0,s fe (t, x, p) , (15) with and 12 The Fokker-Planck collision operator used here describes electron dynamics in plasmas whose temperatures range 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 for selfcollisions between fast electrons and the thermal bulk is particularly important for an accurate quantitative estimate of the rf-driven or Ohmic current source [44,45].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 selfconsistency between the bulk electron temperature and the fast electron energy losses may be obtained, but always requires 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, due to the linearization.If the 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.
where P Z0,s ∆t is the transition probability describing the fact that an electron is at phase space point (x, p) and time t, given that it was at point (x − △x, p − △p) at time t − ∆t, due to a collision with a partially ionized atom.By definition, ´d△pP Z0,s ∆t (x, △p, p) = 1, which states that all electrons are taken into account, irrespective of the initial phase space location (x − △x, p − △p).Here, △p T is the transposed vector of △p, with △p = p − p s , for all p s values and all scattering directions with respect to p directions, where p s is the momentum of the ion of net charge Z 0,s .Since the transition probability is proportional to the product of the elementary crosssection dσ e,Z0,s (p) with the density of targets per unit surface u s ∆tf Z0,s (t, x, p s ), where u s is the relative velocity before the scattering process between the test electron and the atoms of species s with a net ionization state Z 0,s , the friction vector is, while the diffusion tensor is, Here, the Møller relative velocity u s normalized to the speed of light c is given by the relation, where v = p/γ e is the electron (or test particle) velocity 13and v s = p s m e /m Z0,s is the ion velocity of mass m s,Z0,s ≃ m s since m e is much less than the ion nucleus mass m s .Here, p s = p s / (m e c) and p = p/ (m e c).However, since |v s | ≪ |v|, because of the large difference in mass between m e and m s , the relative velocity may be simplified and u s ≃ |v − v s |, even for run-away electrons.Indeed, in tokamak plasmas, for an ion temperature of 5.11 keV, v th s ≃ 2 × 10 −3 for hydrogen and ten times less for tungsten.Møller corrections to u s are therefore always negligible since the energy of these electrons cannot exceed 30 MeV, because of synchrotron radiation losses [47].
The fully screened relativistic Mott cross-section of collision between an electron and an ion of charge Z 2 s,0 is, where x = sin (θ/2), with θ being the usual deflection angle of the electron and r e the classical electron radius.Spin and relativistic corrections are negligible in the non-relativistic limit, and when p 2 ≪ 1, the Mott cross-section merges with the usual Rutherford expression [48,49].Since p = γβ, where β = v = 1 − 1/γ 2 and γ is the Lorentz factor, for 200 keV slide-away electrons, p 2 ≃ 0.93, while for 20 MeV run-away electrons, p 2 ≃ 1610.Knowing that |△p| = 2 |p| sin (θ/2) and that the angular integral is taken over ´dΩ = ´xmax x min sin θdθ ´2π 0 dϕ, where ϕ is the azimuthal angle in the center of the mass frame, the fully screened friction vector A Z0,s and diffusion tensor D Z0,s are, where Υ = 4π r 2 e c, while x Since both definite integrals ´xmax x min . . .dx/x in equations ( 22) and ( 23) give the same value, one obtains, where ln Λ e,Zs,0 = ln(b max / 1 + b 2 min ) is the Coulomb logarithm.Therefore, taking x min = 1/Λ e,Z0,s and x max = 1, equations ( 22) and ( 23) may be approximated by, (25) and The calculation of the Coulomb logarithm ln Λ e,Z0,s is discussed in appendix C. The partial screening is taken into account by replacing Z 2 0,s → Z s − F Z0,s (q) 2 in equations ( 25) and (26), where q = 2p sin (θ/2) = 2px, and following the definition in [7], the Fokker-Planck screening function g Z0,s (p) is defined as, since The formulation (28) guarantees that ´1 1/Λe,Z 0,s |Z s − F Z0,s (q)| 2 dx/x = Z 2 0,s ln Λ e,Z0,s for weakly energetic electrons, while conversely, for very energetic ones, it is Z 2 s ln Λ e,Zs .However, as pointed out in [7], partial screening cannot be described in a strict Fokker-Planck sense other than in the complete and no screening limits.In order to maintain dominant screening terms and avoid unphysical behavior for partial screening, only terms to the lowest order in x must be considered, which allows q to be significant for large electron energies, and consequently take the full form of F Z0,s (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 Λ e,Z0,s in the ultra-relativistic limit.
In the limit of an almost zero ion temperature, as in the post-disruptive regime in tokamaks, s where δ (v s ) is the Dirac function, thus assuming that ions are at rest.In this case, the integration over v s may be performed analytically, and expressions in [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 that A Z0,s and D Z0,s be expressed in terms of Rosenbluth potentials, allowing a convenient conservative formulation of the collision operator.Knowing that s , ∇ v u s = ûs and ∇ v ∇ v u = (I − ûs ûs ) /u s , the integral ´d3 v s may be permuted with the derivatives ∂/∂v, and the term where while with In the coordinate system (p, ξ, φ) used by the LUKE Fokker-Planck solver in momentum space [50], where ξ is the cosine of the pitch-angle, the expression of the collision operator in terms of the divergence of the electron flux in momentum space ∇ p • S coll p (f e ) is, assuming a local axisymmetric plasma.By symmetry, D coll pξ = D coll ξ p = F coll ξ = 0.The contribution of elastic electron-ion collisions to the non-zero diffusion and friction terms D coll pp and F coll p , which describes momentum exchange between particles, is always very small compared to the one of the electron-electron collisions, because of the very large difference in mass between electrons and ions, and may therefore be neglected.The single large non-zero term arising from electron-ion collision is D coll ξξ , which is proportional to Z 2 0,s ln Λ e,Z0,s .Consequently, introducing the partial screening in kinetic calculations requires simply to make the transformation Z 2 0,s ln Λ e,Z0,s → Z 2 0,s ln Λ e,Z0,s + g Z0,s (p) for the pitch-angle diffusion D coll ξξ , with a careful account for the Coulomb logarithm ln Λ e,Zs,0 with Z 0,s , as discussed in appendix C. Here, D coll ξξ incorporates the contribution of all ion species present in the plasma and their respective ionization states.

Inelastic scattering
In the 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 whichconsequently jumps into 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,51], the Abraham-Lorentz-Dirac reaction force for very energetic electrons arising from synchrotron radiation losses [47] and the new one from e − i excitation.The latter can be deduced from Bethe's stopping-power formula describing the losses of energy dE per unit length dx [35,52]: with where hω Z0,s is the mean excitation energy for the ion of net charge Z 0,s .Since the energy loss over a distance △x is equivalent to the work of an effective drag force F excitation p (p) over that distance, its expression in Fokker-Planck calculations in simply, 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 hω Z0,s .In this case, the logarithm term always predominates over the small spin corrections given by the −β 2 term or possible very small additional terms.For electrons whose kinetic energy becomes low compared to the mean excitation energy, the Bethe formula indicates that the stopping power tends to decrease.Within this limit, the dominant inelastic term comes from e − e collisions, as the Bethe formula goes to zero.Unfortunately, the Bethe formula may reverse sign, a non-physical effect that is an intrinsic limitation of the Bethe approach.This problem was identified in [7] and bypassed by performing an interpolation.In this case, ln B Z0,s in equation ( 36) is replaced by ln(1 + B nB Z0,s )/n B , with n B an integer that is chosen heuristically to be 5.When B Z0,s ≫ 1, ln(1 + B nB Z0,s )/n B ≃ ln B Z0,s , and the Bethe formula is well retrieved.The Bethe-like expression guarantees that in the limit p = γβ → 0, ln(1 + B nB Z0,s )/n B − β 2 is always positive and smoothly becomes very small.The exact value of ln(1 + B nB Z0,s )/n B − β 2 is not critical since the dominant inelastic term is from e − e collisions.Another approach, is also to enforce inelastic collisions from the excitation of high-Z elements to zero, when ln B Z0,s − β 2 becomes negative 14 .Both methods are equivalent numerically.
Regarding the formulation of the Fokker-Planck solver in the LUKE code, as shown in equations ( 34) and (35) in [10,47,50], the drag force F excitation p (p), as given by equation (38), may be readily incorporated in F coll p .

Fokker-Planck screening function
From the definition of the Fokker-Planck screening function given by equation ( 28), and making the change of variable y = q/x = 2p/α, where p = p/ (m e c), g Z0,s (p) may be expressed as the sum of two terms g Z0,s,1 (p) and g Z0,s,2 (p) where, and Here, Λ ≡ Λ Z0,s , in order to simplify notations. 14For low-energy electrons, γ The term ln BZ 0,s − β 2 becomes negative if ln BZ 0,s < β 2 , which leads to a transcendental equation in E for determining this threshold.
The form factor given by equation ( 3) may be recast in the simple form, since a Z0,s,i = αā Z0,s,i /2 and αq = q.Therefore, assuming that the condition 2pa Z0,s,i / (αΛ) = pa Z0,s,i /Λ ≪ 1 holds.This is always valid in tokamak plasmas since the Debye sphere has numerous particles.
In much the same way, where and which can be integrated analytically so that, with Gathering all terms, with and, as expected, lim p→0 g Z0,s (p) = 0 is verified regardless of the element and its ionization state Z 0,s, .
In the case of a single exponential corresponding to the standard Yukawa atomic potential, equation (50) simplifies to the usual form: where n = 2, since A 0,s,1 = 1, while g Z0,s,2,i,j (p) = 0 by definition.Using the Thomas-Fermi-Kirillov model [22], n = 3/2, and equation ( 6) in [6] is well retrieved.The difference is generally small for tungsten, a few percent, between the MY and Thomas-Fermi-Kirillov models.
The analytical expression of g Z0,s (p) may be easily implemented in Fokker-Planck solvers, allowing fast and accurate kinetic calculations, whatever p and the type of elements and their level of ionization.As shown in figure 15, g Z0,s (p) from equation ( 50) for neutral tungsten with the use of three exponentials is very close to the numerical estimate g num Z0,s (p) directly obtained in the limit p ≫ 1, from DFT calculations using results of the GAUSSIAN code.The formula is, Figure 15.Normalized Fokker-Planck screening function for the neutral tungsten element, Z 0,s = 0, as a function of the normalized momentum p = p/ (mec), calculated from DFT results (red circles) using the GAUSSIAN code [12] for ρnum Z0,s and equation ( 53) and approximated using the method of moments (multi-Yukawa) given by equation ( 50) with a single exponential (blue dashed line), two exponentials (green dotted dashed line) and three exponentials (pink dashed line). where with and from [7], γ EM being the Euler-Mascheroni constant [53].Since equation ( 53) is derived in the limit p ≫ 1, it is consequently not valid at low p, and g num Z0,s (p) does not converge towards zero when p ⩽ 0.5, as shown in figure 16 for the neutral atom of tungsten.When the number of exponentials is reduced, g Z0,s (p) is always lower than g N Z0,s (p).The relative error is about 13% at p = 1 for a single exponential.This tendency is similar for an ionized atom, as shown for W 42+ in figure 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 Figure 16.Normalized Fokker-Planck screening function for the neutral tungsten element, Z 0,s = 0, as a function of the normalized momentum p = p/ (mec) between p = [0, 1] calculated from DFT results (red circles) using the GAUSSIAN code [12] for ρnum Z0,s and equation ( 53) and approximated using the method of moments (multi-Yukawa) given by equation ( 50) with a single exponential (blue dashed line), two exponentials (green dotted dashed line) and three exponentials (pink dashed line).When p < 0.4, the numerical value of the normalized Fokker-Planck screening function falls off more rapidly than the MY expression with the use of three exponentials, and does not converge towards zero for p ≈ 0, as expected from theory. in section 2, the relative deviation of g Z0,s (p) from g num Z0,s (p) remains small whatever the p value, in the order of a few percent, even if the coefficient of determination R 2 , shown in figure 12, is lower.This results from the fact that g Z0,s (p) is itself an integral, which smooths out possible errors.

Mean excitation energy
The mean excitation energy hω Z0,s 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 hω Z0,s = (1/Z 0,s ) ik f ik ln (hω ik ), where f ik is the dipole oscillator strength of the transition ω ik for the atomic system between quantum states |i⟩ and |k⟩, according to Bethe's theory [34,54].Its determination from first principles calculations is a considerable challenge, so except for elements that do not require relativistic corrections, hω Z0,s is generally obtained from empirical laws constrained by measurements for neutral atoms only [55][56][57][58].Recent advanced calculations carried out by a non-relativistic multi-configurational self-consistent field (MCSCF) code have allowed us to estimate hω Z0,s for all ionization states of the elements lighter than argon Z s ⩽ 18 [14,15,59].Although this result represents considerable progress, the accurate determination of hω Z0,s for many higher-Z elements, such as tungsten, is still missing, which makes it difficult to study the impact of inelastic processes by electronion interaction in a hot plasma.In this context, several simple Normalized Fokker-Planck screening function for the charged tungsten element, Z 0,s = 42, as a function of the normalized momentum p = p/ (mec), calculated from DFT results (red circles) using the GAUSSIAN code [12] for ρnum Z0,s and equation ( 53) and approximated using the method of moments (multi-Yukawa) given by equation ( 50) with a single exponential (blue dashed line), two exponentials (green dotted dashed line) and three exponentials (pink dashed line).models have been introduced to compare their impact on kinetic calculations.
In general, MCSCF calculations show that hω Z0,s has an exponential-like dependence with Z 0,s , which can be easily determined within 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, hω Z0,s has a rather complex structure, while it becomes almost proportional to Z s for Z s > 18, i.e. h ω Z0,s ≃ 10Z s eV, this being known as the Bloch relation [55].For Z s < 18, hω Z0,s oscillates with Z s and tends to increase up to approximately 50% as its value decreases.The Z s dependence of hω Z0,s for neutral atoms can be well described by a statistical approach to the energy loss process, known as the LPA [37,55].In the other limit corresponding to a single bound electron, hω Z0,s = Z 2 s I H eV, where I H = 14.9916 eV is obtained from non-relativistic quantum calculations [14,60].Consequently, hω Z0,s may be approximated by the simple heuristic relation, where Z 0,s is the charge of the fully screened ion, and hω Z0,s is expressed in eV units.This relation can be considered as an upper bound of hω Z0,s , since electron-electron correlations tend to reduce the mean excitation energy [14].This method has been used to quantify the impact of tungsten on the toroidal plasma current driven by RF waves at the lower hybrid frequency in tokamaks [4].It is shown that for argon, in figure 18, the exponential interpolation given by equation ( 57) is in good agreement for both weakly and highly ionized states compared  57); green circles: numerical results from the MCSCF quantum code [14]; blue dashed line: LPA model from equation (59) with γ LPA = √ 2 using the MY atomic model with the use of three exponentials calibrated against DFT calculations (GAUSSIAN code) [37]; red full line: variational quantum model from equation (62) using the MY 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 the NIST database [62].
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 the density of bound charges calculated in the groundstate.Indeed, since excited states are transient with a characteristic time that is 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 hω Z0,s .Two models are interesting for this purpose, the LPA approach [37], dedicated principally to very weakly ionized atoms and a variational quantum description [36].Although 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 MY 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.

LPA
The LPA has been widely used to calculate the mean excitation energies of neutral atoms [55].The LPA formula can be extended to any ionization state, according to the relation, ln hω Z0,s = 4π Ns ˆ∞ 0 r 2 ρ Z0,s (r) ln γ LPA √ 4πα 2 mec 2 ρ Z0,s (r) dr.(58) The LPA formula gives generally poor results when the ionization state is high and a fewer 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 [61].In particular, the choice of γ LPA is rather arbitrary, and its value, from heuristic arguments, is generally set at √ 2. Incorporating equation ( 2) into equation (58) with As shown for argon in figure 18, which is the highest-Z element for which advanced numerical quantum calculations are available [14], a good quantitative agreement is found between LPA calculations using equation (59) and MCSCF ones for the neutral atom.The value from the NIST database is also consistent with the LPA level [62].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 prevents 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 figure 19.

Variational quantum description
A non-relativistic variational quantum model to calculate the mean excitation energy, initially derived for inertial fusion experiments, is an interesting approach to obtain a more accurate estimate hω Z0,s as a function of Z 0,s .According to [36], the mean excitation energy of the ground state is given by the relation ln hω Z0,s = 1 2 ln S (1) /S (−1) where S (−1) = 2m e a 2 0 r 2 / 3h 2 and S (1) = 4K 0 /3.The functions S are moments of the strength distribution of oscillators [54], which may be expressed as a function of K 0 , the averaged kinetic energy of the cloud of bound electrons, and r 2 = (4π/N s ) ´∞ 0 r 4 ρ Z0,s (r) dr or r 2 = 6R 2 from equation (6).Therefore,  59) with γ LPA = √ 2 using the MY atomic model with the use of three exponentials calibrated against DFT calculations (GAUSSIAN code) [37]; red full line: variational quantum model from equation ( 62) using the MY 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 the NIST database [62].
where α 2 m e c 2 ≃ 27.21 eV, which is about twice the Rydberg unit of energy.The calculation of K 0 is performed using the virial theorem, 2K 0 ≃ − U Z0,s , where U Z0,s is the atomic potential related to the density of bound charge ρ Z0,s by the Poisson's equation, as given by equation (1), using the MY description.Therefore, since, one obtains, by reporting equation ( 61) in the expression equation (60).In equation (62), parameters (A 0,s,i , λ Z0,s,i ) are those obtained from the method of moments discussed in section 2.3.For a single exponential, equation ( 62) simplifies to hω Z0,s = α 2 mec 2 λ 3 Z0,s,1 Ns/12, and for neutral atoms, i.e. when Ns = Zs, hω Z0,s = 9.44Zs eV, considering the Thomas-Fermi model for which λZ0,s,1 = 1.13Z 1/3 s , as defined in section 2.2.Using a similar approach based on an approximate description of the Thomas-Fermi model [21], the same value is found, as shown in [36].Both relations are very close to the heuristic Bloch relation, which can also be well reproduced by the LPA model [55].Another approximate description of the Thomas-Fermi model given in [22] gives hω Z0,s ≃ 12.10Zs eV, about 12% larger than the value given by the exact Thomas-Fermi model, but still close to the Bloch relation [55].
As shown for argon in figure 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 MY 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 hω Z0,s to the atomic model in this limit.Indeed, when hω Z0,s is calculated using the approximate Thomas-Fermi atomic models instead of the MY 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 figure 19, ⟨hω Z0,s ⟩ exhibits globally a consistent agreement with the expected limits for a neutral atom and an almost fully stripped one.For argon, the agreement is less accurate near Z 0,s = 0, with a similar relative error.Conversely, the LPA model, hω Z0,s has a correct variation with Z 0,s when its value is close to Zs, making the variational quantum description more appropriate, even if relativistic effects are not considered.
From estimates of hω Z0,s using the LPA and the variational quantum descriptions, both using the MY atomic model for the ground-state, there is a trade-off for an accurate estimate of hω Z0,s whatever Z 0,s would be to consider the largest of the values given by both models.This approach would allow us to accurately describe atomic physics in kinetic calculations, either for a cold plasma such as after a major disruption or for a standard hot magnetized plasma expected during regular tokamak operation.This option is considered in the 3D linearized relativistic bounce-averaged Fokker-Planck code LUKE to study both the physics of post-disruptive runaway and slideaway RF-driven electrons [10,47].

Conclusion
The incorporation of atomic physics in kinetic calculations is becoming mandatory in order to study the impact of high-Z elements in fusion plasmas.The MY approach for describing the atomic potential, regardless of the ionization state, is particularly convenient to obtain consistent analytical solutions for both elastic and inelastic scattering processes that occur in a plasma.This 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 (a few tens of 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 MY atomic potential are unique, as their identification does not rely on a non-linear least-squares 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 type of element, making it universal.While it was initially developed for up to three exponentials, it has been extended to 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 of 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 code, but just by modifying the coefficients of the MY potential.Although simple, the method allows for having a more realistic description of atomic physics compared to simplified atomic models, such as the well-known Thomas-Fermi model and all its approximate representations.
The study presented here has been restricted to microscopic collision processes for kinetic calculations.However, it can be extended to other physical quantities, in particular, those that are derived in the first Born approximation, as already shown for the bremsstrahlung of fast electrons on neutral atoms, where partial screening effects may also be important.It is likely that this method can be extended to many more processes, such as 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.
the GAUSSIAN code (6-311G, cc-pVDZ or AUG-cc-pVDZ, AUG-standing for augmented), as for most atoms that do not require relativistic quantum calculations (up to Krypton approximately).It may also be 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.orgwebsite) [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 online GAUSSIAN code documentation accessible from the https://gaussian.comwebsite.
In the calculations, the spin multiplicity requires special care.It may be obtained from Hund's rules for low-Z elements, but is usually obtained from the NIST database (https:// physics.nist.gov/asd)[62].It is important to note that for six ionization states of 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.Therefore, these ions require 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 https://sobereva.com/multiwfnwebsite [64].

Appendix C. Coulomb logarithm
The Coulomb logarithm ln Λ e,Zs,0 = ln(bmax/ 1 + b 2 min ) from equation ( 24) may be explicitly evaluated, taking into account 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 λ D , the upper limit is bmax = λ D /b 90 .For multispecies plasmas, λ D → λ e−i D ≃ λ e D (1 + Z eff Te/T i ) −1/2 , where λ e D is the usual electron Debye length, and Z eff = s Z0,s n Z0,s Z 2 0,s /ne is the effective charge, knowing that s Z0,s n Z0,s (t, x) Z 0,s = ne from electroneutrality.In most kinetic calculations it is usually assumed that all ion species have the same temperature T i whatever their net charge.
The value of b min depends of the ratio bq/b 90 where bq is deduced from the uncertainty principle.If bq/b 90 ≫ 1, quantum effects predominate and ln Λ q e,Z0,s ≃ ln (λ D /bq), otherwise the classical limit corresponding to b min = 0 can be taken, and ln Λ c e,Zs,0 ≃ ln (λ D /b 90 ) [67,68].Consequently, ln Λ q e,Z0,s is less than ln Λ c e,Z0,s .The quantum limit may be determined from the uncertainty principle △p△x > h/2 where the momentum increment is approximated by △p = µsus, µs = mems/ (me + ms) being the reduced mass between colliding particles.Therefore, assuming △x ≃ bq, the impact parameter is bq ≈ h/ (2µsus) or bq ≈ (λ C /4π) (me/µs) /us, where λ C is the Compton length.The ratio b min = bq/b 90 = us/ (2Z 0,s α) since b 90 = re Z 0,s /u 2 s (me/µs) for Coulomb collisions and λ C /re = 2π/α.A rough estimate of the smooth transition between classical and quantum limits may be obtained by averaging ūs over the electron and ion distribution functions.In this case, the square root of the mean square velocity is ⟨⟨ūs⟩⟩ ≃ 3Te where Te = Te/ mec 2 and the normalized plasma temperature threshold above which quantum effects are significant is T q e,Z0,s = 4Z s,0 α 2 /3.For Z 0,s = 1, T q e,Z0,s = 36 eV, so that the quantum limit must always be taken under standard tokamak plasma conditions with isotopes of hydrogen.For Z 0,s = 42, corresponding to the net ionization of tungsten at Te = 3 keV, then T q e,Z0,s is much larger, T q e,Z0,s = 64 keV, and the classical limit is conversely always valid in tokamak plasmas.Since me ≪ ms, µs ≃ me, and ln Λ c e,Z0,s ≃ ln λ D /re + 2 ln ūs − ln Z 0,s , while ln Λ q e,Z0,s ≃ ln λ D /re + ln ūs + ln (2α).If ln Λ q e,Z0,s and ln Λ c e,Z0,s are both heavily weighted by the ratio λ D /re, the regime dominated by quantum effects concerns principally fast electrons for partially ionized high-Z elements whose kinetic energy is greater than T q e,Z0,s .In the quantum limit, Λ a e,Z0,s is independent of the ion net charge Z 0,s .
In standard MKSA units, with λ D = λ e D , ln Λ q e,Zs,0 = 0.5 ln Te [keV] − 0.5 ln ne 10 20 m −3 + ln ūs + 18.61 and the thermal value is ln Λ q−th e,Z0,s ≃ ln Te [keV] − 0.5 ln ne 10 20 m −3 + 16.04.Using λ e−i D , the Coulomb logarithm must be reduced by the term −0.5 ln (1 + Z eff Te/T i ), and for a pure hydrogen plasma with Te = T i , ln Λ q−th e,Z0,s ≃ ln Te [keV] − 0.5 ln ne 10 20 m −3 + 15.7, a value very close to those found in the literature [67,69].Additional small differences may arise from the choice of the averaged velocity.For relativistic electrons, ūs ≃ v, since ions may be considered at rest, ln Λ q−rel e,Z0,s ≃ ln Λ q−th e,Z0,s + ln p − 0.5 ln Te − 0.5 ln 3, as far as p > 3Te or Ec ≫ 3Te/2, where Ec is the kinetic energy normalized to the electron rest mass energy.Following [7], ln Λ q e,Z0,s may be approximated by ln Λ q e,Z0,s ≃ ln Λ q−th e,Z0,s + ln(1 + (p/( 3Te)) k )/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 Λ c e,Z0,s in the classical limit is weakly dependent on the ion charge Z 0,s and is more sensitive to Te, since ln(Λ c e,Z0,s /Λ q e,Z0,s ) = ln (ūs/ (2αZ 0,s )).

Figure 1 .
Figure 1.Relative fraction of different screened ion charges for tungsten at different plasma electron temperatures using the OPEN-ADAS database [5].

Figure 2 .
Figure 2. Deviation of the bound electron density of from spherical symmetry as estimated by the parameter ∆Θ for all the screened charges of tungsten, from DFT calculations using the GAUSSIAN code[12].When ∆Θ = 0, the cloud of bound electrons is spherically symmetric around the nucleus.∆Θ is very small for all nobel gas-like electronic configurations.

Figure 3 .
Figure 3. Density of bound electrons for all ionization states of tungsten calculated by the DFT method using the GAUSSIAN code [12].Upper red line corresponds to neutral atoms.Details of the simulations are given in appendix B.

Figure 4 .
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 the GAUSSIAN (DFT method, red lines) and GRASP (MCDHF method, blue symbols) codes.Excellent agreement in found between the two quantum relativistic codes.Details of the simulations are given in appendix B.

Figure 5 .
Figure 5.Comparison between densities of bound electrons for neutral tungsten W 0 (upper plot) and the ionized state W +42 (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 λ B 0,s ≃ 0.9Z 0.42 s a −1 0[38] (blue dotted-dashed line) and Yukawa with the inverse screening length λ TF 0,s ≃ 1.13Z 1/3 s a −1 0 [24, 27, 38, 41, 42] (blue dashed line).Details of the simulation for the DFT calculation are given in appendix B.

Figure 7 .
Figure 7. Atomic linear density 4πr 2 ρ Z0,s (r) for neutral tungsten element, Z 0,s = 0, calculated from DFT (red full line) using the 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].Moliere's solution, as given in appendix A is also displayed (black full line)[32].

Figure 8 .
Figure 8. Atomic form factor for neutral tungsten element, Z 0,s = 0, calculated from DFT (red full line) using the 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].Moliere's solution, as given in appendix A is also displayed (black full line)[32].Only the three exponentials case is very close to the DFT solution when q ⩾ 0.08.

Figure 9 .
Figure 9. Atomic density for charged tungsten element, Z 0,s = 42, calculated from DFT (red full line) using the 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].

Figure 11 .
Figure11.Atomic form factor for charged tungsten element, Z 0,s = 42, calculated from DFT (red full line) using the 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 q ⩾ 0.08.

2 [
[min] max = (1 + b max] min ) −1/2 , with 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 2 0,s β −2 .Here, the electron velocity normalized to the speed of light is linked to the Lorentz factor by the relation γ = (1 − β 2 ) −1/2 .The values of b [max] min are discussed in appendix C. In equation (23), the term (I − ûs ûs ) is the usual perpendicular collision operator.

Figure 17 .
Figure17.Normalized Fokker-Planck screening function for the charged tungsten element, Z 0,s = 42, as a function of the normalized momentum p = p/ (mec), calculated from DFT results (red circles) using the GAUSSIAN code[12] for ρnum Z0,s and equation (53) and approximated using the method of moments (multi-Yukawa) given by equation (50) with a single exponential (blue dashed line), two exponentials (green dotted dashed line) and three exponentials (pink dashed line).

Figure 18 .
Figure18.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 equation (57); green circles: numerical results from the MCSCF quantum code[14]; blue dashed line: LPA model from equation(59) with γ LPA = √ 2 using the MY atomic model with the use of three exponentials calibrated against DFT calculations (GAUSSIAN code)[37]; red full line: variational quantum model from equation(62) using the MY 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 the NIST database[62].

Figure 19 .
Figure19.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 equation (57); blue dashed line: LPA model from equation(59) with γ LPA = √ 2 using the MY atomic model with the use of three exponentials calibrated against DFT calculations (GAUSSIAN code)[37]; red full line: variational quantum model from equation (62) using the MY 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 the NIST database[62].
c−th e,Z0,s = 1.5 ln Te [keV] − 0.5 ln ne 10 20 m −3 − ln Z s,0 + 17.69 and ln Λ c−rel e,Z0,s = ln Λ c−th e,Z0,s + 2 ln p − ln Te − ln 3. It can be also approximated by the expression ln Λ c e,Z0,s ≃ ln Λ c−th e,Z0,s + ln(1 + p 2 / 3Tek )/k.The choice of the electronion Coulomb logarithm, ln Λ e,Z0,s p, Te , therefore depends on the type of element s, its local ionization state Z 0,s and the temperatures Te and T i at the same location in the plasma.