Approximate atomic models for fast computation of the Fokker–Planck equation in fusion plasmas with high-Z impurities and suprathermal electrons

The progress in theoretical and experimental development of fusion devices brings constant changes in their design, where the combination of different improvements can raise new issues that were previously neglected. An example is the transition of the material of plasma facing components from carbon to tungsten (Pitts et al., 2013). Tungsten is a high-Z element that is not fully ionized in the plasma temperature range of fusion devices. The presence of such impurities, due to their high electric charge, may have a significant impact on the electron dynamics (Hesslow et al., 2017). Their proper description is crucial for the design of current drive systems or for the mitigation of the runaway electrons by massive gas injection in post-disruptive plasmas (Reux et al., 2015). If the plasma contains a large fraction of fast electrons, such as runaway electrons or suprathermal electrons created by Lower Hybrid Current Drive (LHCD), it might be necessary to include effects of partial screening in collision description, in order to obtain accurate results (Hesslow et al., 2018).

The numerical tools created to model the fast electrons dynamics in modern tokamaks (Peysson et al., 2016; Pokol et al., 2019; Król, 2020) should incorporate the corrections necessary to include partial screening to be able to predict plasma behavior, while keeping the additional computational effort at a reasonable level.

The aim of this study is to compare different models which can be used to evaluate the influence of partial screening on electron collisions. An ab initio model, based on quantum mechanics, is used as a reference in terms of accuracy when comparing different classical and empirical approximated models. This reference model, despite being accurate, is not preferred in practice for numerical implementation to solve the Fokker–Plank equation, because of the associated computation cost. Since the calculations can require multiple integrations of the values acquired from the atomic model (integrations over a few parameters, where some of them vary depending on plasma conditions), it is strongly preferable to use a model which can provide an analytic solution of these integrals. Thus, such approximated models can grant a great improvement in computation speed at the cost of a loss of accuracy of the solution. The goal is to find the optimal trade-off between accuracy and computation time and to provide the best solution for future numerical tools designed for the nuclear fusion community.

When calculating the velocity distribution function with the Fokker–Planck equation, the impact of collisions between $a$ and $b$ species on the distribution $fa$ is described by the collision operator $Cab$⁠. In the special case of the two species, a and b, where the distribution of b is Maxwellian, the collision operator takes the following form (Hesslow et al., 2017):

$Cab$—collision operator for collisions between particle species a and b,$Lfa$—Lorentz scattering operator, $νDab$—deflection frequency, $νsab$—slowing-down frequency, $ν||ab$—parallel-diffusion frequency, $p$—normalized momentum, and $fa$—velocity distribution function of species a.

This description can be applied to many cases of tokamak plasmas, where the ion velocity distribution can be assumed as Maxwellian and the electron distribution is analyzed in detail.

If fast electrons with energies in the range of tens of keV and higher are the object of interest, it is necessary to include the effects of the so-called partial screening in the description of the collisions of these electrons with non-fully ionized atoms. Fast electrons are not simply deflected by the Coulomb interaction with the net charge of the ion, but probe its internal electron structure, such that the nuclear charge is not completely screened. This effect modifies the deflection and slowing-down frequencies in the Eq. (1) and must be taken into account for a more accurate description (Hesslow et al., 2018).

A description, within the Born approximation, of the electron-ion elastic collisions can be obtained by including the atomic form factor for the given ion into the differential collision cross section. Assuming that the ions are infinitely heavy stationary targets initially at rest, the cross section can be expressed as follows (Mott and Massey, 1965):

where $r0$—the classical electron radius, θ—deflection angle, $p$ = γ$v/c$—normalized electron momentum, $v$—impacting electron velocity, $c$—the speed of light in vacuum, γ—Lorentz factor, and $Z$— atomic number.

$Fq$ denotes the atomic form factor, corresponding to the Fourier transform of the electron density distribution $ρr$ in the momentum space, with 0 $≤Fq≤N$⁠. The lower limit $Fq=0$ corresponds to the case of no screening, when the interaction between a free electron and the atomic nucleus is not influenced by bound electrons. The upper limit $Fq=N$⁠, on the other hand, represents complete screening, when the nucleus is screened by all bound electrons. The $Fq$ is expressed as

where: $q=2p sin(θ/2)/α$ is the momentum transfer, $α ≈ 1/137$—the fine structure constant, $r$—atomic radius as a spatial coordinate measured from the center of the atom, $a0$—in the Bohr radius, and N—the number of bound electrons.

The electron density is normalized such that the number of bound electrons is

Electron densities used in this work are always spherically averaged, so only the radial coordinate $r$ is used in model descriptions, instead of a full 3D representation in $r$⁠. As for $q$ and $p$⁠, only the magnitudes are the values of interest. Hereafter they will be denoted as $q$ and $p$⁠, respectively, to simplify the notation. In addition, all quantities used in calculations are expressed in atomic units by default unless otherwise noted. In the case of the spherically symmetric electron density (r), the atomic form factor defined by Eq. (3) can be expressed as

where: $q=q$ and $r=r$⁠.

The atomic form factor is related to the electron density of bound electrons $ρr$ around the atom nucleus. Therefore, the calculation of the atomic form factor strongly depends on the model used to describe the electron density distribution $ρr$ in the atom. The most accurate description requires full quantum relativistic calculations, including the electron exchange and correlation effects. Such calculations are computationally demanding and the obtained solution accuracy is greater than required for applications like solving the Fokker–Planck equation (Peysson et al., 2021; Król et al., 2021). However, there are some semi-classical models, which can provide approximated solutions with far less computational effort. They can provide a significant acceleration of the numerical codes used for solving the Fokker–Planck equation without losing the demanded accuracy.

In the present work, in order to better understand and describe the effects of the used atomic model on the calculation of the atomic form factor, such simple models are compared against ab-inito calculations. In the subsequent part of the article, the Density Functional Theory (DFT) is used as the reference model. Simpler models are based on two approaches and their modifications, namely, Thomas–Fermi (TF) and Pratt–Tseng (PT).

DFT is a method to investigate the electronic structure of matter in a wide range of scales (from atoms to solids and large organic molecules). Its basic principle is that the ground state energy of the system is a unique functional of the electron density $ρr$⁠. DFT is based on two theorems by Hohenberg and Kohn (1964). With Kohn–Sham (KS) equations (Kohn and Sham, 1965), it is possible to calculate fictitious orbitals of a non-interacting system, which generates the same electron density as the interacting system. Thus DFT is an ab initio quantum mechanical model providing results as close to reality as possible and is used here as a reference for models comparison.

In the present work, calculations using DFT were performed with the Gaussian software package (Frisch et al., 2019). The hybrid exchange-correlation functional PBE1PBE (known also as PBE0) (Perdew et al., 1996; Adamo and Barone, 1999), along with the natural orbital-relativistic correlation consistent basic set ANO-RCC, was used. The results from Gaussian were then post-processed in the Multiwfn program (Lu and Chen, 2011), which was used to obtain the spherically averaged electron density $ρr$⁠. An example of the resulting electron density for tungsten ions is presented in Fig. 1. The shell structure of the atom can be clearly seen in the results for different ionization stages.

The TF model is a well-known semi-empirical approximation for many electron systems in the external potential $V(r)$⁠. It cannot reproduce quantum effects such as the shell structure, but is relatively simple and uses the electron density $ρr$ as a primary quantity—in this way it was a precursor to DFT. In general, the TF model treats bound electrons as an electron gas. This gas, despite being inhomogeneous in atoms or molecules, can be treated locally with relations developed for the homogeneous electron gas. Such approach makes it possible to eliminate the necessity to calculate electron wave functions and allows calculating the electron density directly. The major drawback of this approach is that the used approximations do not hold in the vicinity of the nucleus, at which the electron density becomes singular.

According to the TF model, the electron density $ρr$ can be calculated as (Lundqvist and March, 1983)

where

The electron density vanishes outside of the ion radius $ρTFr>rmax=0$⁠.

The value of $φ(x)$ can be obtained from the Thomas–Fermi equation

with the following boundary conditions:

In practice, a solution of the Thomas–Fermi equation was found with an iterative approach, because $xmax$ is not known for ions. The bisection method was used to find the $xmax$ value which satisfied the boundary condition: $Z−NZ+φ′xmaxxmax=0$ within 0.000 05 margin. The search was made for $xmax$ between 0.001 and 150. This method is very simple and robust, the only drawback was the computational efficiency. However, in this case it was not a major issue since values of the electron density were calculated only once. Moreover, this case is simple enough to be solved on the order of seconds. The boundary condition for neutral atoms must be changed and instead of $φ∞=0$⁠, $φ$ must vanish for some finite value of $x,$ i.e., $φxmax=0$⁠. In this work, calculation with $xmax=400$ was performed. Influence of the chosen $xmax$ on the solution of the Thomas-Fermi equation was analyzed in other works, e.g., Jardin et al. (2020) and based on them, the selected value is considered as a sufficient approximation. An example of the resulting electron density calculated in the framework of the TF model is presented in Fig. 2. One of the most visible differences between DFT and the TF model is the lack of a shell structure in the TF model.

Over time, many modifications to the TF model were developed to improve its accuracy, by including additional electron interactions and relativistic effects. One of the most basic improvements is the Dirac correction, which introduces a term for exchange energy of electrons.

According to the Dirac modification, the modified Thomas–Fermi–Dirac (TFD) equation takes the following form (Csavinszky, 1969):

with the following boundary conditions:

with the appropriate modification, the final expression of the electron density $ρr$⁠, within the TFD model, takes the following form (Abrahamson, 1961):

where $b=r/x$—factor defined in the TF model.

A solution of the TFD equation is found with an iterative approach, in a similar way as in the previous case of the TF equation. The difference is that, due to the change in the boundary conditions, the TFD model no longer requires a solution at $xmax=∞$ in the case of the neutral atom; thus no approximation for the $xmax$ is needed. An example of the resulting electron density is presented in Fig. 3. The TFD model is very similar to the TF model, but in the TFD the electron density is discontinuous at the ion radius.

The PT model is an alternative method of deriving an approximated analytical expression for the atomic form factor. It is based on the Yukawa potential, approximating the electrostatic potential around the ion nucleus as the sum of the classical Coulomb term and an exponential term accounting for the screening of bound electrons at short distances from the nucleus. The electron density calculated in this way is a rough approximation of the electron density calculated with DFT, nevertheless its analytical expression allows straightforward calculation of the atomic form factor. The electron density $ρr$ is expressed within the PT model as

Four versions of this model are investigated with different approaches for calculating the $a$ parameter,

where $x=Z−NZ$⁠.

Other parameters are defined differently in every version of this model, as presented in Table I. The PT₁ model is the modified version of PT₂, where the parameter $k$ was modified to better fit DFT results. PT₂ is the original model proposed by Botto, for which coefficients are based on the fit to the Hermann–Skillman potential. PT₃ is based on the transfer of the analytic results from the Kirillov approximation (Kirillov et al., 1975) of the TF model into the framework of the PT model. The PT4 model was created for the purpose of bremsstrahlung calculation and features modified parameters to better fit the interior of the atom.

The PT models provide closed sets of equations with no need for any iterative process in order to calculate the electron density. In addition, they allow creating an analytical solution for the atomic form factor, as it will be shown later. Thanks to that, it is definitely the most appealing model when considering computation time optimization, granting a few orders of magnitude faster computation than other models. Its computational simplicity takes a price of less accurate results, but still on a satisfactory level. An example of the resulting electron density calculated using PT model 1 (Botto I) is presented in Fig. 4. A significant difference can be seen in comparison not only to DFT results, but also TF and TFD.

The differences between the previously defined electron density models for the tungsten ion W¹⁰⁺ are depicted in Fig. 5. It shows that, close to the nucleus, the electron density calculated based on both PT and TF models is in a quite good agreement with the results obtained using DFT. However, when considering values of the electron density in the region further from the nucleus, PT models give rather poor approximation of the DFT results. The results from DFT show the effect of the shell structure of the atom by steep changes of the electron density. Unfortunately, none of the considered simplified models can reproduce this effect, which is related to the underlying classical assumptions.

For the most accurate electron density distribution obtained from DFT, there is no analytical solution of the form factor. Therefore, the form factor is numerically integrated, which is a time-consuming operation. In this work, the form factor was calculated for 60 001 values of $q$ for each ion, which guaranteed negligible numerical errors. For practical reasons, it would be possible to use a much coarser grid of $q$ with just a few thousand of points, but for this analysis it was important to avoid any influence of numerical error. Even with the detailed representation of the form factor, the most computationally expensive part was the calculation of DFT electron density, which for tungsten ions can take up to a few minutes on 20 CPUs node of a HPC cluster. Integrations, to obtain form factors and later deflection frequency, were made on the order of seconds. For some models, it is possible to derive an analytical solution of the above integral or at least its approximation, granting a significant improvement in computation time. The solution can be later used to create analytical solutions of other integrals needed to calculate collision frequencies. Thank to this, the final result for elastic collisions for a given ion can be obtained on the order of thousandths of a second on a single CPU.

The TF model does not have an analytical solution, so to calculate the form factor based on this model, it is necessary to use a similar procedure as in the case of DFT. It means that there is no significant benefit from using this model for calculation of the form factor. It is however possible to obtain an approximated solution of the form factor calculated with the TF model, by using some approximations during the integration, as described by Kirillov et al. (1975). This allows to calculate the form factor using the following equation:

Unfortunately, using this approximation provides results which are not so close to DFT calculation as when using the exact TF solution with numerical integration presented in Sec. IV. A comparison of the form factors calculated based on the Kirillov approximation with those calculated by numerical integration of the electron density obtained within TF, TFD, and DFT models is presented in Fig. 6. The results show that a simple classical model can give reasonable approximation of the atomic form factors. The TFK approximation is less accurate than the exact TF model, but it is much faster to compute.

In the case of the PT models, it is possible to give an analytical expression for the form factor, using the electron density in the form defined by Eq. (5),

where $a$ is an individual parameter for every ion as defined in Table I. Its value is common for both the electron density calculation and the form factor calculation. Unfortunately, again the solution obtained with this analytical expression is not close to the form factor obtained with numerical integration of the TF or DFT (but it gives results with similar accuracy to Kirillov approximation). The results of the form factor calculations based on four PT models compared to the form factors calculated within numerical integration scheme of the TF and DFT-calculated electron densities are presented in Fig. 7.

For the presented calculations, the estimates of the atomic form factors $F(q)$ were stored in the computer memory as a table of numerical values with $q$ logarithmically distributed in the range between 10⁻⁶ and 10⁺⁶ with 60 001 samples. To create a figure of merit allowing a quantitative comparison of different approaches to the atomic form factor calculations, the root mean square (RMS) of the difference between the normalized form factors calculated within the DFT formalism (reference) and the TF/TP models was calculated for every ion. The difference between form factors was normalized between 0 and 1 by dividing $F(q)$ by the number of electrons bound to the nuclei. This approach allowed to obtain a similar range of values for every ion, while still using all advantages of the RMS error calculation. Sensitivity to outliers is an advantage in this case, because it promotes solutions which have similar accuracy for all values. The equation for the RMS defined above takes the following form:

where: $n$—number of records used for storing the form factor in computer memory, $FDFT$—value of the first form factor used as a reference for a given $qi$ value from ith record, $F2$—value of the second form factor used in comparison, and $N$—number of bound electrons in the considered ion.

The calculated values of the RMS mentioned above are presented in Figs. 8 and 9 for tungsten and argon ions, respectively. Tungsten was chosen because of its use as a material for plasma facing components (Brezinsek et al., 2019), and the argon ions are presented due to their importance in experiments for mitigation of runaway electrons (Reux et al., 2015). In the case of tungsten, the results obtained by integration of the electron density obtained from TF models are the closest to the DFT results for most ions. This proves that the classical model can approximate quite well the quantum-mechanical model, but gives no further benefits, as calculating procedure is of the similar level of complexity. The Kirillov approximation of the TF model provides less accurate results than the exact TF solution, comparable to the PT models. From the models which have analytical solutions of the form factor, the first version of the PT models (Botto I) seems to be the best choice for the ionization stages which can be expected in the tokamak plasma core. The Kirillov approximation, on the other hand, seems to be a better choice for lower ionization stages, which can be found in plasma edge or post-disruptive plasmas. Taking this into consideration, the choice of the most accurate model should be made only after definition of the most important cases, e.g., elastic collisions are most sensitive for partial screening in the case of low ionization stages.

In the case of argon, the results are similar—the direct TF model proves to be the most accurate in most cases. Only for a low number of bound electrons the Kirillov approximation is more accurate. The PT₁ model shows the best performance for intermediate ionization states. The conclusion is therefore similar that the choice of the most accurate model depends on the ions of interest.

None of the models investigated in Sec. IV is clearly better than others in all regions of interest. Because it was expected that room for improvement exists, an attempt was made to extend some of these models in order to achieve an increased accuracy while maintaining short computation time. Two modifications of the form factors calculation are presented in Secs. V A and V B.

In the case of the PT models, it is possible to give an analytical solution of the form factor using the electron density calculated with more than one exponential. In this way, different parameters for every shell can be used (Peysson et al., 2021). The electron density can then be calculated with the following equation:

with the corresponding analytical solutions of the form factors given as

where $N=∑i=15Ni$ and $Ni$ is the number of electrons included in every summation part, $ai$ are parameters corresponding to each group of electrons.

This means that it is possible to obtain more accurate results of $ρr$ and $Fq$ by calculating $ai$ separately for different groups of electrons and then summing up the contributions with different $ai$ as in Eqs. (18) and (19). In this work, it was decided to use up to five electron groups, each one representing approximation of an electron shell in the atomic structure. The exact number of electrons which can be included in every group is 2, 8, 18, and 26, and the rest are in the 5^th group, as presented in Table II. The first two groups represent exactly the electrons on the first and second shell, respectively. The third group includes the third electron shell plus the 4s subshell, because the latter has similar energy to the 3d subshell and the order of filling of these two subshells differs between elements. The choice of the following groups follows similar logic. The fifth shell appears in elements which have as few as 37 electrons, because electrons on the 4f subshell have similar energy as electrons on the 5s and 5p subshells. Because of that, filling of the mentioned subshells is quite irregular and the configuration of xenon, which is the last element without the electrons on the sixth shell, is chosen as a limit for the fourth group. The heaviest element considered in the scope of this work is tungsten and only ground states of the ions are taken into account. It means that the sixth shell can be occupied at most by two electrons. Including these two electrons into a separate group would complicate the model while bringing little improvement in terms of accuracy, so it was decided to include them in the group with electrons on the fifth shell. Ions with few electrons have less electron groups, as some electron shells remain unoccupied.

Based on the equation used for calculating $a$ in the different versions of PT models, the following approximation of $ai$ coefficients was proposed:

where $c1,i, c2,i, c3,i, c4,i$ are constants which need to be defined for each electron group. In order to find the optimal values of these parameters, a code was written in Matlab to perform optimization with a gradient descent with momentum algorithm (Qian, 1999). The results were fitted to the atomic form factor calculated by integration of the DFT-calculated electron density. As the cost function of the optimization algorithm, a RMS function as presented before was used.

The fitting procedure was as follows: first, the parameters $c1,1, c2,1, c3,1, c4,1$ were optimized for ions of selected elements, with bound electrons only in the first group. With these parameters defined, the parameters $c1,2, c2,2, c3,2, c4,2$ were optimized for ions with a minimum of 3 and a maximum of 10 electrons, i.e., electrons from the second group. With the parameters from the first two groups defined, the next set of parameters for the next group was optimized using ions with the corresponding number of electrons, and so on. This procedure required five optimization runs to define all parameters. Each time, the optimization of the parameters was performed until the RMS function change between iterations was smaller than 10⁻⁷ and the change of the optimized parameters between iterations was smaller than 10⁻⁵. The momentum parameter of the gradient descent algorithm was set to 0.8. The learning rate was in each case manually adjusted by increasing it if the convergence was too slow or decreasing it in the case of divergence or strong oscillations.

As a result, the optimal parameters found for Eqs. (20b) and (20c) are presented in Table III. The values of $c3,i$ and $c4,i$ equal to 1 were rounded to full number, because in these cases the values of $xns, i+1$ from Eq. (20a) is reaching zero and increasing the values of $c3,i$ and $c4,i$ does not change the result.

Such approach is in fact a search for minimum of the very complicated function of 20 parameters. It is almost impossible to determine if the minimum found is a global minimum, or just a local one. During investigations performed in the scope of this work, it was observed that obtaining better parameters for some specific ions can reduce the accuracy when calculating the form factor for other ions. The presented set of parameters grants reasonable accuracy for a wide range of elements, though it was created with focus on high-Z elements. In almost all cases, it gives better results than the models presented in Sec. IV. Nevertheless, it would be possible to obtain different sets of parameters if someone would like to focus on the different range of elements or ions.

To check if the $ai$ values found in this way are close to the best results which can be obtained with the general multi-exponential form of $F(q)$ as in Eq. (19), an additional fitting procedure was performed. The $ai$ parameters were fitted separately for every tungsten ion, without using Eq. (20a). It means that an optimization process was run 74 times with the $a1, a2, a3, a4, a5$ parameters optimized each time to find the best fit of the resulting form factor to the results obtained with DFT electron density.

The results of the atomic form factor approximation created as described in Eq. (20a) (called PT_opt in this work) are presented in Fig. 10. The line described as “Fit” presents the result obtained with fitting the $ai$ directly, without defining them with universal equations Eqs. (20b) and (20c). The vertical lines mark the electron groups used in multi-exponential Yukawa potential definition.

The obtained results show a significant improvement in terms of RMS. Furthermore, time calculation effort is on the similar level as for other PT models presented before. When comparing the PT_opt model with the direct fit of $ai$⁠, the equations for $ai$ used in the PT_opt model obviously generate results less accurate than when using directly fitted $ai$ parameters; however, they provide the ability to calculate the results without storing the fitted parameters for every ion. The resulting form factor can be considered as a relatively good approximation of the DFT model.

As shown previously, the form factor formulation in the case of PT and Kirillov approximation has a common form. It is therefore possible to develop a more general definition with two free parameters that will give the same or better accuracy of the form factor. Based on form factor definitions for PT and TFK in Eqs. (15) and (16), a generalized expression of the form factor is proposed as

where $a$ and E were used as free parameters in a fitting procedure to match the resulting form factor with the form factor based on the DFT results.

In this way, a set of two parameters for every ion is enough to calculate the form factor. It is however impossible to derive analytically the electron density from Eq. (21). The other drawback of this method is that it can only be applied to the ions for which the electron density was previously calculated with DFT and the parameters must be found using the fitting procedure. No extrapolation of the model is possible. Despite this, it nevertheless provides a fast way of computing the atomic form factor based on a limited number of parameters. The results of the proposed approaches are compared in Fig. 10 to the other models.

Values of parameters obtained for the FE model for tungsten ions are presented in Figs. 11 and 12. In the case of E, the value of 1.5 is used in the TFK model, and the value of 2 is used in PT models. Values of a obtained from PT₁ and PT₂ are given for comparison.

The effect of partial screening on the deflection frequency between electrons and ions in Eq. (1) was calculated as in (Hesslow et al., 2017). The goal of these calculations was to show the influence of partial screening on elastic collisions, therefore instead of calculating the absolute values, the ratio between deflection frequency calculated with complete screening and partial screening is used. It is defined with the following equation:

where $νDei$ is a deflection frequency calculated with the partial screening, $νD, csei$ is a deflection frequency calculated with complete screening, $Z0=Z−N$⁠, $lnΛ$ is the Coulomb logarithm and $gp$ is defined as

where $x=qα/2p$⁠.

As the Coulomb logarithm was not the value of interest in this case, it was assumed constant for all tungsten ions $lnΛ=16.3$⁠. It corresponds to conditions in typical tokamak plasma core, with plasma electron density n_e = 5 × $1019 m−3$ and electron temperature T_e = 3 keV.

After integration, Eq. (23) takes the following general form for models which can be generalized with Eq. (21):

where $y=2a p/α$⁠.

In the case of the multi-exponential model described by Eq. (19), the integration gives the following result:

The greatest impact of the partial screening can be observed for weakly ionized atoms and high momentum of the colliding free electron. A comparison of the deflection frequency ratio calculated with different models for tungsten W¹⁰⁺ ion is presented in Fig. 13. The DFT model is calculated with numerical integration of $gp$ with Eq. (23). The other models are calculated with $gp$ as in Eq. (24) or Eq. (25), depending on the model. The thermal and suprathermal ranges of the electron energies are marked in the plot to make the analysis easier. The boundary between these two values was set arbitrarily at 20 keV. The lower limit for the thermal region is set at 10 eV, because at this temperature majority of the electrons will be bound in hydrogen atoms. A vertical line marks the electron momentum for which its kinetic energy is equal to the ionization energy of the tungsten W¹⁰⁺ ion which is about 0.2 keV.

To quantify the average performance of different models, the RMS of the difference between DFT based results and different models was calculated for tungsten. The obtained values are presented in Table IV in two versions: the RMS of the absolute and relative difference, respectively, between a given model and the DFT. The resulting quantities for PT_opt are about four times smaller than in the case of other models.

Most of the approximations of the electron density (r) presented in Secs. III and V A can be used in the calculation of the slowing-down frequency. It can be done with the Local Plasma Approximation (LPA) proposed by Lindhard and Scharff (1953), which uses the electron density as the main quantity to estimate the mean excitation energy $I$⁠. This approximation was derived mainly based on the TF atomic model, treating bound electrons as a homogeneous electron gas. It does not have a strong theoretical basis from the point of view of modern atomic physics, but it can still produce some useful results. The mean excitation energy $I$ can be used for the calculation of the inelastic collisions with bound electrons and, in this case, is a primary parameter describing plasma. According to the Bethe stopping-power formula (Hesslow et al., 2017), the slowing-down frequency $νsee$ in Eq. (1) for the Fokker–Plank collision operator is modified as follows:

where $ne$—density of free electrons in plasma, $ni$—density of ions in plasma, $β=v/c$⁠, $h=γ−1mc2/I$ and $lnΛee=lnΛ+1/k ln1+2γ−1/pTe2k/2$⁠, where $pTe$—thermal momentum, $k$-model parameter.

The mean excitation energy $I$ can be calculated with the LPA using the following form of $lnI$ for ions:

where the corresponding local plasma frequency:

and $e$—electron charge and $m$—mass of the electron.

The values of $lnI$ for neutral atoms were calculated using LPA with different models of electron density. To define their accuracy, the results were compared against experimental values from the NIST database (National Institute of Standards and Technology, 1995). The plots of $lnI$ as a function of the atomic number Z for neutral atoms calculated using different models are presented in Fig. 14. To provide more quantitative comparison, the RMS of the difference between LPA models and the NIST database were calculated. Results are presented in Table V in two versions corresponding to the RMS of the absolute and relative difference, respectively, between a given model and the NIST database. For most elements, LPA can be used to obtain approximation of the mean excitation energy. Result obtained based on the PT_opt model are similar to results based on the DFT calculation.

The values of $lnI$ for ions of some elements were calculated using different models of electron density and compared to the results obtained with Multi-Configurational Self-Consistent Field (MCSCF) ab initio quantum mechanical calculation by Sauer et al. (2015). Results for argon are presented in Fig. 15. Despite the fact that LPA is not well suited for ions, it can be noted that the PT_opt model gives very similar results to DFT electron density and is even able to reproduce some features of the atomic shell structure. However, all the models based on LPA underestimate $lnI$ for high ionization stages—which is not surprising given the assumptions made in the model.

Different classical approaches to calculation of the electron density distribution and atomic form factors were compared with the ab initio results obtained from the DFT method in order to find a computationally efficient way of calculating the atomic form factors that can be implemented into Fokker–Planck equation solvers. The existence of an analytical solution of the integrals needed to calculate form factor and deflection frequency was an important requirement, because it gives the possibility to achieve a much shorter computation time. As none of the investigated models was considered as satisfactory enough for this purpose, an attempt to extend the existing model was made, in order to achieve better approximation of the DFT-based results. A definition of the multi-exponential Yukawa potential model was presented, which can achieve better accuracy than most of the investigated models, while still providing advantageous computational efficiency. In addition, comparison to some fitting models was presented, as an alternative approach to computing the form factors.

The application of the presented models in calculation of the plasma collision parameters was presented. The calculation of the elastic collisions between free electrons and ions shows that use of the proposed approximation should provide quite accurate results, while significantly reducing computation time. The impact of including partial screening can result in a few times higher collisions frequency than while using the complete screening model, it is thus non-negligible. The effect is more important in the case of ions with high number of bound electrons and high energy of colliding free electrons. These two conditions are a priori in contradiction, because higher electron temperatures would result in stronger ionization of atoms. Nevertheless, in the case of runaway electrons in relatively cold post-disruptive plasmas, these conditions are common and the impact of partial screening should be significant. Regarding microwave heating or current drive that can result in the local production of suprathermal electrons, the description of the interaction with non-fully ionized impurities, without the partial screening, could result in an underestimation of the elastic collision frequency and even more significantly of the fast electron Bremsstrahlung (Peysson et al., 2021; Król et al., 2021). The presented models can be implemented into calculation of the fast electron bremsstrahlung in plasma, which also requires atomic form factor integration. In this case, the possibility to find an analytic solution would also be very important, as bremsstrahlung depends on many parameters and numerical integration would be a more time consuming solution.

It is planned to include created models into suite of codes used for calculations of plasma dynamics, to verify the influence on plasma properties. Impact of the partial screening on the runaway electron generation, current drive efficiency, and entropy generation in plasma should be investigated.

In the case of inelastic collisions of free electrons with electrons bound to the ion, the proposed models show similar results to the DFT based calculation. Unfortunately, the LPA approach used in this work is a significant source of uncertainty by itself. More reliable methods of calculating mean excitation energy for high-Z element ions will be a subject of further investigation. This is a region so far mostly unexplored, because experimental measurements of the stopping power and related mean excitation energy are usually restricted to neutral atoms. Values for different ionization stages are acquired from ab initio models, which meet difficulties in the case of high-Z elements.

The further work is planned to investigate different approaches to obtain mean excitation energy. Application of some ab-initio methods, such as Multi-Configuration Dirac–Hartree–Fock (MCDHF), Coupled Cluster (CC) or Time-Dependent Density Functional Theory (TDDFT) should provide better insight into this topic as well as more accurate results.

The presented calculations should help to better understand and simulate the fast electrons population dynamics in plasma with non-fully ionized impurities. This will allow a precise control of the energy deposition in the case of the LHCD and electron cyclotron resonance heating, as well as more accurate strategies of runaway electron mitigation. The precise description of the electron distribution can also improve understanding of other processes in plasma, like tungsten transport. It should help to enhance the plasma confinement in tokamaks and bring us closer to efficient nuclear fusion.

This work has been partially funded by the National Science Centre, Poland (NCN) Grant HARMONIA 10 No. 2018/30/M/ST2/00799. We thank the PLGrid project for computational resources on the Prometheus cluster. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 and 2019–2020 under Grant Agreement No. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.