Numerical study of limits of neoclassical theory in the plateau regime in the presence of impurities

This paper presents a numerical simulation to test the corrections to the radial electric field in the presence of impurities for tokamak devices whose main ions are in the plateau regime. The effects due to the presence of impurity density variation in a pedestal tokamak are known to generate conditions where conventional neoclassical theory is no longer valid. It is found that for a specific axisymmetric concentric magnetic field, the presence of steep radial gradients for both temperature and density translates in a correction to the radial electric field. We anticipate that comparisons to the Landreman–Fülöp–Guszejnov model display good agreement for both the analytic approximation and numerically solved value of the normalized density when low impurity concentration and atomic charge state cases are considered, coupled with values of the ordering parameter Δ≤2$$ \Delta \le 2 $$ . Higher values of the latter reveal better agreement to the neoclassical Hazeltine–Hinton formula, suggesting a discrepancy with the validity threshold of the model.


INTRODUCTION
Since the early days of research on toroidal fusion devices, much effort has been devoted to the study of collisional effects that undergo the dynamics of charged particles in fusion devices.When dealing with the description of a plasma through the various analytical models at our disposal, particular care has to be regarded to the approximation used, as their domain of applicability is heavily reliant on the ordering of physical parameters determining the phenomena at hand.This applies particularly to the transport theory in magnetically confined plasmas.The neoclassical theory is able to account for the magnetic field non-uniformities, and the resulting rich orbit topologies introduces the orderings  =   ∕L ⊥ ≪ 1, , where   is the poloidal ion gyroradius, L ⊥ the radial scale length associated with the density and temperature profiles,  = v th ∕qR is the transit frequency, and E∕(Bv th ) = O().Here v th is the ion thermal velocity, p the pressure, E and B are the electromagnetic fields with corresponding safety factor q. The first-ordering ensures that the length of the poloidal radius is much smaller than the scale length for macroscopic parameters, while the following rules out nonequilibrium processes and rapid fluid phenomena.In tokamak plasmas, it is often the case that a very steep gradient of temperature and pressure may not be accommodated in the conventional neoclassical theory domain, and as such one is forced to relax the associated orderings.Yet the task to construct a theory under these assumptions incurs in nontrivial complexities, similarly as one considers, for example, the -ordering to not be satisfied, and nonlocal transport fluxes emerge. [1]As in impure tokamak plasmas, the development of poloidal variation is firstly manifested by the impurity density n z , the neoclassical requirement of  to be the smallest parameter in the transport problem cannot be always met, as in the plasma edge.To reconcile the theory, one allows  to increase while satisfying the neoclassical ordering, allowing in turn poloidal asymmetries to arise.Such choice translates in assuming the parameter Δ ≡ ν ii Z 2 of order unity.We designate Z as the impurity charge state number, νii = L || ∕ ii the ion collisionality,  ii the mean-free path for the bulk ions, and L || the connection length.[4] From a physical standpoint, Δ expresses the ratio between the ion-impurity friction force and the parallel impurity pressure gradient. [2]As will be shown later, the presence of impurity ions is such that the flux surface averaged radial electric field roughly converges to the Hazeltine-Hinton formula [5,6] in the limit of vanishing impurities, while in the opposite limit, it does not as the steep radial gradients of the density and temperature profiles generate a poloidal variation in the impurity density.This variation can lead to a substantial influence of the impurity ions on the plasma flow when the bulk ions are in the plateau collisionality regime and the impurities in the Pfirsch-Schlüter collisionality regime. [4]We will restrict our attention to a transport problem for an axisymmetric concentric magnetic field with orderings  ≪ 1, Δ ∼ 1, and Z ≫ 1 while assuming a large aspect ratio  = r∕R ≪ 1, in line with the analytical models offered in literature. [1,3,4]Several numerical tools have also been developed for investigating neoclassical physics either with Eulerian approach [7] or with PIC approach. [8,9]In the present manuscript, we are comparing the analytic theory to the numerical simulations with the results of full-f code ELMFIRE. [9]he remainder of the paper is organized in the following manner.Section 2 summarizes the analytical neoclassical transport theory results of Ref. [4] applicable to regimes of high impurity concentration for background ions in the plateau regime.Section 3 describes the setup of the simulation to test the validity of the analytical model and presents the numerical results obtained.The results show that for low impurity concentration, the approximated solution of the normalized density yields a better agreement of the radial electric field when compared to its full analytical form.Finally, the paper is summarized and the results discussed in detail in Section 4.

ANALYTICAL MODEL
[3] When the gradients are chosen to be sufficiently steep, the impurities are pushed toward the inside of each flux surface while the particle flux will develop nonlinear terms proportional to the gradients.As the dynamics of the impurities becomes nonlinear, the value of the radial electric field changes by means of the impurity ions density variation on each flux surface and consequently deviates from the Hazeltine-Hinton formula [5,6] : In the remainder of the paper, we will assume Z i = 1 as we consider the gamma factor in the latter expression, to approximate the dependence over all collisionality regimes on the ion collisionality parameter  * i = √ 2rB 0 ∕B p v th,i  i  3∕2 as derived by [10] Here, v th,i = √ 2eT i ∕m i is the ion thermal velocity, 4 ln Λ the ion collision time, Z i the impurity ions atomic charge state, and ln Λ the Coulomb logarithm.
As the radial gradients of the density and temperature profiles become more steep, the impurity ions are induced to develop a poloidal variation which substantially affects the current plasma flow.We will restrict ourselves to the case where the bulk hydrogen ions are in the plateau regime while the impurities are in the Pfirsch-Schlüter regime.In order to perform a comparison to (1) and to our numerical simulations data, we here express the radial electric field of Ref. [4]  in quasi-toroidal coordinates (see Appendix A).With this in mind, as shown by Landreman et al., [4] the impurities radial electric field satisfies the equation where I = R 0 B 0 for our field choice (A1); see Appendix A, B(r) = b √ < B 2 > is the field magnitude, with b 2 = 1 − 2 cos  while K z (r) is proportional to the poloidal impurity rotation.
Here n = n z ∕ < n z > is the normalized impurity density, while represents the ratio of the temperature and pressure scale lengths while prime variables are to be considered radial derivatives.Here, the factor B p R 0 in Equation ( 3) arises from expressing derivative in quasi-toroidal coordinates, d∕d = ( 1∕B p R 0 ) d∕dr; see Appendix B for further details.
The coefficient y is computed by requiring the ambipolar condition for the particle fluxes to be satisfied.Due to the electron ions mass ratio being small, such condition reduces to a more tractable form, which mathematically reads Γ i + zΓ z = 0.A thorough derivation of the ion and impurity fluxes may be found in Ref. [4].In line with our notation choice, its theoretical expression will be where we have chosen  = ⟨n z ⟩Z 2 T e ∕ ( n i,0 T e + n e,0 T i ) , following the same reasoning as for (4), while the coefficients  * and g read In the latter formulas, ln Λ is the ion-impurity collision time.Note that the equivalence between (6) and the expression for  * defined in Ref. [4] follows directly from the identity 2n 0 ∕T 0 = n e0 ∕T e + n i0 ∕T i .If we were to apply the limit n z → 0, the expected pure plasma limit, y → 1 would be recovered.In such limit, the poloidal impurity rotation (4) becomes while the electric field (3) reduces to The absence of the electric charge in the latter equation is due to the temperature being expressed in eV.Note that given the small difference in magnitude between the magnetic fields B and B 0 , order of 1%, convergence of ( 9) to its neoclassical counterpart ( 1) is guaranteed for plateau regime ions.Given the general analytical form of the parallel momentum equation for impurities, see eq. ( 20) in Ref. [4], it is clear that a solution to the former may be achieved only for specific limits.Relevant to our case of study is the limit of a plasma with small inverse aspect ratio.The impurity density may be expanded in powers of the inverse aspect ratio  by means of the expression , where both coefficients n c and n s are of O().As shown in Ref. [4], both coefficients may be computed by solving the appropriate equation and read It will be shown later in the manuscript that the -expansion of n, has the effect of smoothing the contribution of the impurity density to the radial electric field and consequently increases the convergence to the conventional neoclassical value.
Given the Equations ( 10) and (11), depending on the steepness of the ion density or temperature gradient, an up-down and in-out asymmetry of the impurity density is generated while being in the Pfirsch-Schlüter regime.In contrast when in the plateau regime, the sign of the asymmetry is heavily reliant on the magnitude of the parameters ,  * and y, a detailed account on the matter can be found in Ref. [1,2,4].The asymmetry in the impurity density is closely interlinked to the poloidal impurity rotation exhibited in both regimes, which as one might expect will differ from the conventional neoclassical theory.As a matter of fact, if one assumes the impurity density gradient to be nonconstant on each flux surface, considering the main ion to be in the plateau regime one finds a modified expression for the impurity poloidal rotation where A straightforward limit procedure shows that in the absence of impurities, the value of V pl z converges to the conventional neoclassical prediction, as X → 1.

NUMERICAL RESULTS
In the following numerical analysis, we refer to FT-2 tokamak plasmas as our benchmark. [11]In the first subsection, we compare the theoretical predictions of Section 2 for a plasma composed of hydrogen main ions and carbon impurity species in various distinct charge states, while in the second subsection, we compare the analytical results with the numerical results of the particle-in-cell code ELMFIRE.

Comparison of analytic equations
In our numerical simulations, we set a = 0.0785 m, R = 0.553 m, B t = 2.2 T, and I p = 55 kA, where with a, we refer to the plasma minor radius, R the major radius, B t the toroidal magnetic field, and I p the plasma current.The initial radial profiles for the temperature and density are defined as hyperbolic functions n, T = n 0 , T 0 ( where r 0 = (r L + r R )∕2 is the mid-radius with r L = 0.02 m and r R = 0.07 m being the radial limits of the simulation domain, T 0 = 120 eV and n 0 = 5 × 10 19 m −3 are the temperature and density at r = r 0 and gradient scale length X 0 , while L X = X∕(dX∕dr) = 0.05 m is identical for both temperature and density.In the present results, we assume the temperature ratio between ions and impurities to be ∼ 1.
Figure 1 exhibits the radial density and temperature profiles at the beginning of simulation.In both plots is rather apparent the presence of steep gradients for both quantities, in line with the theoretical framework introduced in Section 2. In particular, the impurity concentration has been selected to be much lower with respect to the ion and electrons counterpart.Such choice derives from the domain bounds of the model, where we require the parameter Δ = ν ii Z 2 to be of order of unity, while  =   ∕L ⊥ , expressing the ratio between the poloidal Larmor radius of the bulk ions and the radial scale length, to be much lesser than 1.Here,

F I G U R E 2
Radial profiles of the Δ = ν ii Z 2 and X parameters for different atomic charge states.
In Figure 2, the dependency of the parameters X and Δ upon the charge state number Z depicted.The former, as introduced previously through the formula (13) serves as a correcting factor to the poloidal impurity rotation, as compared with the conventional neoclassical prediction.The latter parameter is responsible for the change in ordering needed to apply the analytic theory of Ref. [4] aimed to account for large gradients when impurities are considered.The requirement Δ ∼ 1, as shown by Figure 2a, is what allows to relax the conventional neoclassical ordering Δ ≪ 1, given the parameters chosen to run the numerical simulations.As expected, both parameter values agree with the plots presented in Ref. [4].In Figure 3, the parameter  is depicted for different charge states as well as the main ions collisionality.The plasma main ions lie in the plateau collisionality regime according to the initial profiles as shown by Figure 3. Figure 4 depicts the radial electric field (3) and the poloidal velocity in the plateau regime (12) for different charge states Z of the carbon impurity ions.In the first plot, the increase of the impurity poloidal velocity for higher values of the atomic charge state is traced to the mass conservation requirement, which in turns enforces a rearrangement of the impurity ions on each flux surface.In a similar fashion, the second plot makes apparent the dependency of the radial electric on the temperature and density gradients, which as introduced in Ref. [1-3] generate an increase in the impurity density asymmetry.

F I G U R E 3
Radial profile of the product of the collisionality of the plasma main ions  * i with  3∕2 and of the parameter  for different atomic charge states.As shown by the plot, the main ions are in the plateau regime being  * i  3∕2 < 1.

(km/s) (kV/s)
F I G U R E 4 Radial profiles of the poloidal velocity and electric field in the plateau regime for variable atomic charge states.

Comparison to numerical simulation results
As stated at the beginning of Section 3, a numerical comparison of the radial electric field between the Hinton-Hazeltine form, [12] the Landreman-Fülöp-Guszejnov form, [4] and the one simulated through ELMFIRE has been performed to test the analytical equations of Section 2. In the following, we tested the Landreman-Fülöp-Guszejnov form in the two cases.The first case, referred to as the "weak density approximation"-case ("Weak Approx" in the figures), involves the -expansion with poloidal density dependencies ( 10)- (11).In order to obtain the variable y, a numerical computation is performed by substituting Equations ( 10)- (11) to Equation ( 5) at each radial position.Since y(n s , n c ), n s (y), and n c (y), the value of y in framework of this expansion is solved by Newton's iteration.This expansion requires [4] n − 1 ≪  which is not generally true as the density variation (see Figure 5) can be close to local value of .Thus, it is worthwhile to test the validity of expansion numerically.In contrast, the second case, referred to as the "No approximation" case, utilizes the numerical density data obtained from the ELMFIRE simulation.This choice tests both the analytical and numerical tool at the same time.
The code used to perform the numerical simulations, ELMFIRE, [9] has been designed at Aalto University.ELMFIRE is a full-f gyrokinetic particle-in-cell (PIC) code that is able to simulate both neoclassical and turbulent physics.The static magnetic equilibrium with concentric circular flux-surfaces is used according to a user-defined current profile.Collisions

F I G U R E 5
Poloidal density variations for different values of the parameter  and g corresponding to the cases 2B (red) and 1D (blue) in the left plot and 1C (red) and 1B (blue) in the right plot; see parameters in Table 1.
are evaluated using a momentum-conserving binary collision model [13] while the momentum conservation interpolation for the electric field is based on [14] algorithms, which enable the inclusion of neoclassical effects.An electric field solver is used in the code as it allows to obtain ion polarization from implicit treatment of the polarization drift velocity, without requiring explicit use of the quasi-neutrality equation as polarization terms in the Poisson equation cancel to order O ( (k ⊥ ) 2 ) ; see Ref. [15].Further numerical details can be found in recent paper, [16] while a revised set of equations of motion for the new version of the code can be found in Ref. [17].In this section, we employ the same physical parameters as in 3.1.simulate quasi-stationary state of the neoclassical electric field which is compared to the analytic equations of Section 2. The code and its predecessor have been already earlier been used for studying neoclassical poloidal rotation, [18,19] the effects of poloidal asymmetries in neoclassical physics [20] as well as poloidal localized particle sources. [21,17]For our scope, we run ELMFIRE in "neoclassical mode" by means of a coarse resolution (N  = 100) in the poloidal direction.Such resolution enables the study of density poloidal variations effects affecting the conventional neoclassical physics while prohibiting turbulence growth as the proper scales are not resolved.At the same time, loop voltage and LH heating which were used in heated cases [22] are here turned off to investigate the quasi-stationary state development of the neoclassical radial electric field.In the present manuscript, the impurity concentration was selected to have similar -parameter values as in Ref. [4].A total of 25.7 million electrons (1070 per cell on average) and an almost equal amount of main ions were simulated with a number of impurities ions totaling 1.5 × 10 6 (64 per cell on average).As a first test of ELMFIRE simulation, we look at the poloidal variation of density from the code for different -values as, according to Ref. [4], this has clear effect of inboard-outboard symmetry.As shown in Figure 5, we get similar peaking of impurity density toward outboard (inboard) equator for  ≈ 0.7 ( ≈ 0.2) as observed in fig. 1 of Ref. [4].After that, we compared the radial electric fields obtained from the four different approaches.Figure 6 shows a comparison of the radial electric field for a plasma with fully ionized (Z = 6) carbon impurity species.The ordering parameters  = 0.361 and  = 0.321 were selected for case a and case b.Good agreement between all radial electric fields was found for the cases L n = 5 cm and L t = 20 cm with parameters Δ = 1.11 and  = 0.235.Conversely, for the cases L n = 10 cm and L t = 5 cm, with Δ = 1.38 and  = 0.209, more discrepancies were found in the results for values of r∕a < 0.60.Such results are in contrast to the ones depicted in Figure 8, where a comparison between a plasma in the plateau regime with partly (Z = 2) and fully (Z = 6) ionized carbon is shown, for three times higher density n 0 = 1.5 × 10 20 m −3 corresponding to a normalized collisionality of  * i = 15.In this case, the ordering parameter  = 0.49 was selected for both cases.Relatively good agreement between simulation and all three theory versions was found for both cases Z = 2 (Δ = 0.63 and  = 0.02) and Z = 6 case (Δ = 5.3 and  = 0.2).This leads to infer that the numerical results obtained replicate fairly well the radial electric field behaviour of the Landreman-Fülöp-Guszejnov model, for low impurity concentration, atomic state and ordering parameter Δ close to unity.Moreover, even with the higher n 0 , the plasma is still in the plateau regime.The gradient lengths were increased,  1.
which in turn yielded a lower value of the ordering parameter , thereby providing additional support to the analytical theory presented in Ref. [4].In Figures 7, 9, and 10 we modified the temperature and density gradients scale lengths around the base case to modify the Δ value.Here, we take the average of the E r profiles over 2 cm in the mid-radius where the gradients are steepest and plot these values as a function of  =   ∕L ⊥ for different L ⊥ values, where The ordering parameters for each case are collected to Table 1.As expected, the results start to deviate when  increases due to the steepening of the density and temperature profiles.What is quite surprising is that the numerical results agree relatively well with the conventional Hinton-Hazeltine estimate for all cases with steep gradients ( ≥ 0.4)  which was also the case in the numerical testing of Hinton-Hazeltine formula in Ref. [18].At that time, the predecessor of ELMFIRE which had similar neoclassical physics was used and showed good agreement of results in the plateau regime.For smaller  s -values, the numerical results agree better with the Landreman-Fülöp-Guszejnov estimates.In the Z = 6 case, see Figure 7, the difference is large for all cases where relative difference does not change much as a function of gradient scale length but for lower -values is in the favor of Landreman-Fülöp-Guszejnov with 10%-20% difference to ELMFIRE value.For Z = 2-case in Figure 9, agreement is good between all four curves even up to  ≈ 0.5.In the cases with very gently sloping gradient (3A-3E) leading to small -values in Figure 10, absolute difference between different

CONCLUSION AND DISCUSSION
In this paper, we have presented our simulations results related to the findings of Ref. [4] for a specific axisymmetric concentric magnetic field where the presence of large gradients affects the functional form of the radial electric field for a non pure plasma.We have discussed how in the pure plasma limit, the conventional neoclassical Hazeltine-Hinton predictions are recovered.We obtained good agreement with the Landreman-Fülöp-Guszejnov form for small impurity concentration and atomic charge state.The comparisons of the radial electric field expression in Ref. [4] yielded better agreement for values of the Δ parameter in the interval [0.5,2] for both the analytic approximation and numerically solved value of the normalized density, where the use of terms averaged over the numerical data from the simulations displayed good agreement with the electric field value from that same simulation.This may indicate that the validity threshold of the model orderings is more restrictive.As far as we know, there has not yet been any numerical validation of the analytic expressions presented in Ref. [4]; therefore, before making further conclusions, we encourage to verify the results presented in the manuscript by means of other independent numerical simulation tools to assess the limits of the analytic theory and to confirm the present numerical findings.

6 F I G U R E 7
Radial profiles of the electric field for charge state number Z = 6, base case (a) L n = 10 cm, L t = 5 cm and (b) L n = 5 cm, L t = 20 cm.(kV/m) (kV/m) Radial averaged E r over mid-radius for Z = 6 by varying (a) density and (b) temperature gradient length over values L n,T = 5, 10, and 20 cm corresponding to the cases (a) 1C-1B-1A and (b) 1E-1D-1A in Table

8
Radial profiles of the electric field for base case L T = L n = 5 cm for charge state numbers (a) Z = 2 and (b) Z = 6.Here, three times higher density, n 0 = 1.5 × 10 20 m −3 is used.(kV/m) (kV/m) F I G U R E 9 Radial averaged E r over mid-radius for Z = 2 by varying (a) density and (b) temperature gradient length over values L n,T = 2.5, 5, and 10 cm corresponding to the cases (a) 2C-2A-2B and (b) 2E-2A-2D in Table 1.Here, density is n 0 = 1.5 × 10 20 m −3 .

(
kV/m) (kV/m) F I G U R E 10 Radial averaged E r over mid-radius for Z = 6 by varying (a) density and (b) temperature gradient length over values L n,T = 10, 20, and 40 cm and keeping the other gradient length fixed L n,T = 20 cm corresponding to the cases (a) 3C-3A-3B and (b) 3E-3A-3D in Table1.Here, density is n 0 = 1.5 × 10 20 m −3 .
15213986, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/ctpp.202300066 by Aalto University, Wiley Online Library on [04/04/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Radial profiles of the plasma species density and temperature.
F I G U R E 1 15213986, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/ctpp.202300066 by Aalto University, Wiley Online Library on [04/04/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Parameter values for the test cases in Figures6-10.the same level as in other cases, leading to larger relative differences.Here, "No approximation" version of Landreman-Fülöp-Guszejnov equation gives best agreement with the ELMFIRE result.