Prevention of core particle depletion in stellarators by turbulence

In reactor-relevant plasmas, neoclassical transport drives an outward particle flux in the core of large stellarators and predicts strongly hollow density profiles. However, this theoretical prediction is contradicted by experiments. In particular, in Wendelstein 7-X, the first large optimized stellarator, flat or weakly peaked density profiles are generally measured, indicating that neoclassical theory is not sufficient and that an inward contribution to the particle flux is missing in the core. In this Research Letter, it is shown that the turbulent contribution to the particle flux can explain the difference between experimental measurements and neoclassical predictions. The results of this Research Letter also prove that theoretical and numerical tools are approaching the level of maturity needed for the prediction of equilibrium density profiles in stellarator plasmas, which is a fundamental requirement for the design of operation scenarios of present devices and future reactors.


INTRODUCTION
The stellarator is the main alternative to the tokamak in the quest for achieving net energy from controlled nuclear fusion, due to its intrinsic steady-state operation and absence of current-driven instabilities.Nonetheless, neoclassical transport-associated with the inhomogeneity of the magnetic field and with particle collisionshas traditionally handicapped the performance of stellarators.For this reason, its minimization is one of the main targets in stellarator optimization.Wendelstein 7-X (W7-X) has demonstrated that large stellarators with optimized neoclassical transport can be designed [1] and built with high accuracy [2] and that this has been fundamental for achieving triple product record values [3].
However, even in optimized stellarators, neoclassical transport may play a major role in the plasma core [4].In this region, high temperatures are achieved, causing strong neoclassical particle and energy fluxes due to their unfavorable scaling with the temperature in lowcollisionality plasma regimes.For example, in the 1/ν regime, the neoclassical particle and energy fluxes scale as T 7/2 s and T 9/2 s , respectively, with T s being the temperature of species s (see, e.g., Ref. [3]).In what follows, we only deal with plasmas consisting of singly charged ions and electrons, so that s can take the values s = i and e, respectively.
Stellarator plasmas with core particle transport dominated by neoclassical processes have been predicted to present problems of core particle depletion [6,7].If the core particle source is negligible, as is the case in large stellarators with peripheral fueling, the steady-state particle balance equation for each of the species reads Γ s = Γ neo = 0.This relation, which additionally implies that no net radial current exists in the plasma, i.e., that the particle fluxes are ambipolar, imposes a constraint on the plasma profiles.The neoclassical particle flux can be written as = 0, arrive at the constraint on the plasma profiles [7].Namely, Reactor-relevant plasmas are expected to display peaked temperature profiles, dT s /dr < 0, with T e ≈ T i .Since δ s 12 > 0, this automatically leads to a hollow density profile dn e /dr > 0. Because δ s 12 > 1, (i.e., because the so-called thermodiffusion coefficient L s 11 δ s 12 is relatively large), the hollowness can become large enough to make the pressure gradient positive.
The above analytical argument illustrates the generality of the mechanism of core depletion by neoclassical transport.Of course, more accurate relations between the stationary density and temperature gradients can be obtained by means of numerical simulations.Figure 1 shows the neoclassical radial particle flux Γ neo s , obtained with the neoclassical code KNOSOS [8], over a wide range of values for the density and electron temperature gradients.In Fig. 1 and in the remainder of this Research Letter, we define a/L Ts := −a d ln T s /dr and a/L ns := −a d ln n s /dr.Peaked (hollow) density profiles have a/L ns > 0 (a/L ns < 0).For the rest of the parameters, typical values from W7-X hydrogen plasmas sustained by Electron Cyclotron Resonance Heating (ECRH) are employed (see the caption of Fig. 1 and Ref. [9]), and E r is set by ambipolarity of the neoclassical fluxes.Outward fluxes (Γ neo s > 0) are obtained in most of the represented area, driven by a/L Te , indicating the importance of thermodiffusion.Indeed, the contour Γ neo s = 0 (white region) predicts a slightly hollow stationary density profile even in the presence of a practically flat electron temperature profile.For the larger values of a/L Te obtained in W7-X, a distinctive hollow density profile should be measured if particle transport were described by neoclassical theory only.
Despite this robust neoclassical prediction, flat or weakly peaked density profiles have generally been measured in ECRH plasmas of W7-X (see, e.g., Ref. [10]).This disagreement between neoclassical theory and experiments may indicate that a significant inward contribution (Γ s < 0) to the particle flux has not yet been identified.Under this hypothesis, this Research Letter addresses the particle transport problem in W7-X from the perspective of microturbulence and by means of a large number of nonlinear gyrokinetic simulations in stellarator geometry, carried out with the code stella [11].First, a parameter scan is performed, demonstrating that the turbulence resulting from the onset of the ion and the electron temperature gradients produces an inward pinch in broad regions of parameter space.Then, a specific W7-X discharge is studied, performing turbulence simulations throughout the plasma radius.The resulting profile of the turbulent particle flux is compared with the shortfall in the flux inferred from a careful particle balance analysis that includes estimates of the particle source and computations of the neoclassical particle flux.The sign of the missing contribution to the particle flux needed to sustain the experimental density profile agrees with that of the calculated turbulent particle flux, not only at the core but also at the edge of the plasma.In particular, the simulations consistently predict a sign change of the turbulent flux at intermediate radial positions, reflecting the fact that the neoclassical flux is too large at the core and too small at the edge of the plasma to explain the experimentally determined particle flux.

PARAMETRIC DEPENDENCE OF THE TURBULENT PARTICLE FLUX
Let us write the radial turbulent particle flux as the sum of a diffusive and a convective contribution, with D being the diffusion coefficient (which, in general, depends on the density gradient) and V being the convection velocity.As D is always positive, the diffusive term adds a contribution to the flux with opposite direction to that of the density gradient.In other words, the diffusive flux is inwardly or outwardly directed depending on whether the density gradient is hollow or peaked, respectively.On the other hand, V can be either positive or negative, supporting the formation of hollow or peaked density profiles, respectively.More specifically, the sign of V for a flat density profile (a/L ns = 0) is of particular interest, as it predicts whether a hollow (if V > 0) or peaked (if V < 0) density profile will develop.In order to gain insight into this matter, and prior to the comparison between numerical and W7-X experimental results, we provide a comprehensive numerical characterization of the turbulent particle flux, with emphasis on the case with a flat density profile, i.e., a/L ns = 0.
The turbulent transport is modeled with the flux-tube δf gyrokinetic code stella [11], which has been extensively benchmarked for W7-X geometry [12] and applied to the study of turbulent impurity transport in this device [13,14].The gyrokinetic simulations are nonlinear, collisionless, electrostatic and account for kinetic ions and electrons.The calculations are performed at r/a = 0.25 for the standard W7-X configuration.Due to the low magnetic shear, we employed generalized twist-and-shift boundary conditions [15].Further information on the resolution, boundary conditions, flux tube employed and other details of the gyrokinetic simulations presented in this Research Letter, can be found in Appendix A.
We carried out a parameter scan around the reference point in parameter space {T e /T i , a/L Ti , a/L Te , a/L ns } = {1, 3, 3, 0} in order to study the dependence of the turbulent particle flux on the following quantities: a/L ns in Fig. 2(a), the electron-temperature-to-ion-temperature ratio T e /T i in Fig. 2(b), the normalized ion temperature gradient a/L Ti in Fig. 2(c), and the normalized electron temperature gradient a/L Te in Fig. 2(d).Such exhaustive scans cover both reactor-relevant conditions, where ions and electrons are expected to be thermalized and have comparable values of density and temperature gradients, and conditions corresponding to the first W7-X campaigns, where the temperature and its gradient are noticeably larger for electrons than for ions in the core.For the sake of clarity, the reference point in parameter space is represented with a golden square.
Starting with Fig. 2(a), note that turbulence driven solely by either the electron or the ion temperature gradient (being a/L Ti = 3 or a/L Te = 3) cannot lead to peaked density profiles, as convection is outwards (V > 0) for the flat density profile cases (a/L ns = 0).In contrast, when both temperature gradients coexist, turbulence gives rise to a significant inward convection (V < 0), being Γ turb s /Γ gB i ≈ −0.6 at a vanishing density gradient.Moreover, the normalized density gradient at which the particle flux is zero is a/L ns ≈ 0.7.In other words, if all transport were due to turbulence (and for the fixed values of other parameters used in this plot), the equilibrium density profile would be clearly peaked.From Fig. 2(b), it can be observed that a larger electrontemperature-to-ion-temperature ratio enhances the particle pinch.Furthermore, Fig. 2(c) shows that a/L Ti need not be large or comparable to a/L Te in order to turn the −0.02 0.00 0.02 0.04 particle flux inwards.The negative convection starts to manifest from values as modest as a/L Ti ∼ 0.25 when a/L Te = 3, and its magnitude increases roughly linearly with the size of a/L Ti .Similarly, Fig. 2(d) shows that, although the case with only a/L Ti = 3 drives positive flux, adding an electron temperature gradient as moderate as a/L Te ≳ 1.0, turns the particle flux negative.
In summary, the existence of an inward turbulent particle flux is very robust and is found over broad regions of the scanned parameter space, whenever the temperature gradients of ions and electrons are both finite and not too small.In Ref. [16] it is found that this turbulent pinch holds for different stellarator devices, which aligns with other references that report inward turbulent particle fluxes in LHD [17,18].
Finally, we would like to shed some light on why the electron temperature gradient is key for obtaining inward turbulent particle fluxes.Particle fluxes are driven by density and electrostatic potential fluctuations that are out of phase.For example, if the electrons were treated adiabatically, the density and potential would be in phase, and no particle fluxes would be driven.The sign of the particle flux is thus correlated with the sign of the phase shift between density and potential fluctuations.Figure 3 shows the spectra (along the normalized binormal wave number k y ρ i , where ρ i is the ion Larmor radius) of the turbulent particle flux (top row) and of the phase difference between density fluctuations, δn k , and electrostatic potential fluctuations, φ k (bottom row).The spectra are represented for the case where a/L Te = 0 (left column) and a/L Te = 3 (right column), both for a/L Ti = 3 and for a/L ns = 0. Inward particle flux contributions occur when the density fluctuations lag behind on the potential fluctuations (negative values of the phase difference), which is what predominantly happens at almost all represented scales for the case with a/L Te = 3.In contrast, the phase difference is essentially positive at the scales with largest contributions to the total particle flux when a/L Te = 0.

PARTICLE BALANCE ANALYSIS FOR A W7-X PLASMA DISCHARGE
To determine the radial particle flux in an experiment, one needs to compute the neutral ionization source in the confined plasma.In current fusion experiments, the recycling neutrals-i.e., the plasma that is neutralized at the divertor plates and returns to the main plasma-are generally the main source of particles.In Ref. [19] it was concluded that both the main chamber and divertor targets contribute to the fueling of the plasma core in W7-X ECRH plasmas.Once the neutral hydrogen atoms reach the confined plasma they are transported towards the core by charge exchange reactions with the plasma ions, while they ionize by electron impact.
The radial profile of the neutral density, n 0 (r), and the associated ionization source are estimated using a short-mean-free-path, one-dimensional neutral transport model given by (see, e.g., Ref. [20]) where D CX = T 0 /m 0 ν CX is the charge exchange diffusivity coefficient, T 0 and m 0 are the neutral temperature and mass, and ν ion and ν CX are the temperature-and density-dependent ionization and charge exchange frequencies.Neutrals and plasma ions are assumed to be isothermal, T 0 (r) = T i (r), due to frequent charge exchange reactions between them.Equation ( 4) determines the neutral profile up to a multiplicative constant that can be determined from the global particle confinement time τ p .In Ref. [19], a confinement time of τ p = 0.258 s was estimated for a similar discharge using a single-reservoir plasma particle balance.The neutral transport model has recently been compared with neutral density measurements in W7-X (see Fig. 14 in Ref. [9]), showing varying degrees of agreement for τ p between 0.15 and 1.0 s.In Fig. 4(b), τ p is varied from 0.1 to 1.0 s to account for uncertainties in this number.
Let us denote by Γ exp s the flux profiles calculated, for this range of values of τ p , by integrating the neutral ionization source ν ion n 0 over the plasma radius.The neutral density n 0 is obtained from (4) using the experimental profiles shown in Fig. 4(a).The results for Γ exp s are displayed in Fig. 5(a).There is a region in the core (the width of which depends on the assumed τ p ) where the estimated central particle sources are too small to sustain the outward neoclassical radial flux Γ neo s calculated with DKES [21].One is led to conclude that an additional particle flux must exist which moves particles towards the magnetic axis in that region.Further outside, on the contrary, the additional flux must be outward directed and account for essentially all the particle flux, since the neoclassical contribution falls one to three orders of mag-nitude short with respect to Γ exp s .A similar shortfall is observed in other stellarators, such as HSX [22,23].reasonably well, taking into account, among other factors, the uncertainties in the plasma profiles and the simplicity of the neutral transport model employed.The effect of the uncertainties of the plasma profiles on the turbulent particle flux is quantified in Appendix B. In addition, in the future, it would be desirable to carry out this kind of analysis using fullflux-surface or radially global gyrokinetic codes.

CONCLUSIONS
Whereas neoclassical theory predicts strongly hollow density profiles in the plasma core of large stellarators, measured density profiles in Wendelstein 7-X are flat or weakly peaked, indicating that theoretical particle flux calculations are missing an inward contribution.In this Research Letter, we have shown that the discrepancy between theory and experiment disappears when the turbulent component of the particle flux is included in the calculations: Turbulence driven by both the ion and electron temperature gradients provides the missing inward particle flux needed to sustain the experimentally measured density profiles.Moreover, we have performed a one-dimensional particle balance analysis in a Wendelstein 7-X discharge where the sign of the turbulent particle flux agrees, over the entire plasma radius, with the sign of the difference between the experimental and the neoclassical particle fluxes.In particular, the simulations predict a sign change of the turbulent flux, accounting for the fact that the neoclassical flux is too large at the core and too small at the edge of the plasma to explain the experimental particle flux.The phase-space coordinates employed by stella are denoted by {x, y, z, v ∥ , µ}.The coordinates perpendicular to the magnetic field lines (Taylor expanded to first order in the flux tube formalism) are defined as x = r−r 0 and y = r 0 (α − α 0 ), where r 0 and α 0 define the center of the flux tube, r = a ψ t /ψ LCFS is the effective minor radius, a is the minor radius, ψ t is the enclosed toroidal flux divided by 2π and ψ LCFS is the toroidal flux at the last closed flux surface.The field line label is defined as α = θ p − ιζ, with ι being the rotational transform and (ζ, θ p ) being straight-field line coordinates, which are chosen to be the geometrical toroidal angle ζ and the poloidal PEST coordinate θ p [24].The parallel coordinate z = ζ is taken to be the geometrical toroidal coordinate.The velocity coordinates are the component of the velocity parallel to the magnetic field, v ∥ , and the magnetic moment, µ s = m s v 2 ⊥ /2B, where m s is the mass of species s, v ⊥ is the perpendicular velocity and B is the magnitude of the magnetic field.
A Fourier decomposition of the electrostatic potential is performed in the directions perpendicular to the magnetic field, for which radial and binormal wave numbers are defined with minimum values k x,min = 2π/L x and k y,min = 2π/L y , determined by the extent (L x , L y ) of the perpendicular simulation domain.The simulations are performed with a resolution of (N x , N y , N z , N µ , N v ∥ ) = (91, 271, 49, 12, 48).The perpendicular box size in real space is (L x , L y ) = (94ρ i , 94ρ i ), which corresponds to a box in Fourier space of (k x,max ρ i , k y,max ρ i ) = (2,6).Here ρ i = v th,i /Ω i is the ion Larmor radius with v th,i = 2T i /m i being the ion thermal velocity, Ω i = eB r /m i being the ion gyrofrequency, T i and m i being the temperature and mass of the ions, and B r = 2 |ψ LCFS |/a 2 being the reference magnetic field.In stella, lengths in the perpendicular direction are normalized with respect to the ion Larmor radius ρ i .The particle flux is normalized with respect to the gyro-Bohm particle flux given by Γ gB i = n i v th,i (ρ i /a) 2 , where n i is the density of the ions.
The generalized twist-and-shift boundary conditions [15] are employed to deal with the low shear of Wendelstein 7-X.Instead of relying on the global magnetic shear at a given radial location, we will take advantage of the local magnetic shear, which varies along the field lines.Thanks to this flexibility, the length of the flux tube can be chosen such that the local magnetic shear at the ends of the flux tube imposes k x,min = k y,min .For the simulations in this Research Letter, the flux tube extends approximately one poloidal turn to ensure an equal spacing in the perpendicular wave number grid.The flux tube for which the simulations of this Research Letter are run, is commonly known as the bean flux tube, defined by α 0 = 0.It is centered with respect to the equatorial plane, θ p = 0, and the bean-shaped toroidal plane, ζ = 0.

APPENDIX B: UNCERTAINTIES IN THE PLASMA PROFILES AND GRADIENTS
To quantify the effect of the uncertainty in the plasma profiles [Fig.4(a)], gyrokinetic simulations are performed with a 10% (or 20%) change in the plasma parameters for r/a = 0.4, the results of which are shown in Fig. 6.Since stella performs simulations with dimensionless parameters, the uncertainties in the plasma density n i = n e and the ion temperature T i are visualized by denormalizing the particle fluxes accordingly.Figure 6 suggests that at r/a = 0.4, higher density gradients and lower temperature gradients would improve the agreement between Γ exp s − Γ neo s and Γ turb s .Similarly, lower densities and temperatures would improve agreement, although their influence is smaller compared with the effect of changing the gradients.

6 FIG. 1 :
FIG. 1: Neoclassical particle flux as a function of a/L Te and a/L ns for the standard W7-X configuration [5].Here, r/a = 0.25, n e = n i = 6.3 × 10 19 m −3 , T e = 2.2 keV, T i = 1.1 keV and a/L Ti = 0.81.The circle corresponds to the experimental discharge analyzed in the final part of this Letter.The dashed line shows the combination of a/L Te and a/L ns that gives Γ neo s = 0 in Eq. (2).

s
+ Γ turbs = 0. Here, Γ s is the flux-surface-arXiv:2209.04194v3 [physics.plasm-ph]8 Nov 2023 averaged radial particle flux, and it is assumed to be the sum of a neoclassical contribution Γ neo s and a turbulent contribution Γ turb s .Assuming Γ turb s = 0, the above steady-state conditions reduce to Γ neo i = Γ neo e

FIG. 2 :
FIG. 2: Turbulent radial particle flux in gyro-Bohm units as a function of the (a) normalized density gradient, (b) electron-temperature-to-ion-temperature ratio, (c) normalized ion temperature gradient, and (d) normalized electron temperature gradient.The point with T e /T i = 1, a/L Ti = a/L Te = 3, and a/L ns = 0 (golden square) is shared by the four represented scans.

FIG. 3 :
FIG. 3: Spectra of the particle flux (top) and phase difference between the density (δn k ) and electrostatic potential (φ k ) fluctuations (bottom), with a/L Te = 0 (left) and a/L Te = 3 (right).The dominant phase difference for each k y ρ i is denoted by a dot.Positive values are highlighted in red and negative values in blue.
FIG. 5: Γ exp s , Γ neo s and Γ turb s for discharge No. 20180920.017.The green error bars (a) and the gray squares (b) correspond to ±20% variations of the plasma gradients.The turbulent particle flux Γ turb s is modeled with stella and shown in Fig. 5(b) in black.Its sign and size agree with Γ exp s − Γ neo s

FIG. 6 :
FIG. 6: Effect on the turbulent particle flux of 10 and 20% changes in the plasma parameters at r/a = 0.4.The values of Γ exp s − Γ neo s and Γ turb s shown in Fig. 5 are displayed in red and black, respectively.The error bar on Γ exp s − Γ neo s represents the uncertainty in τ p , while the error bar on Γ turb s represents the standard mean deviation of the saturated state of the time trace of the simulation.
The authors acknowledge fruitful discussions with Félix I. Parra.This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 101052200 -EUROfusion).Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission.Neither the European Union nor the European Commission can be held responsible for them.This research was supported in part by Grant No. PGC2018-095307-B-I00, Ministerio de Ciencia, Innovación y Universidades, Spain, and by Grant No. PID2021-123175NB-I00, Ministerio de Ciencia e Innovación, Spain.