Asymmetry between deuterium and tritium turbulent particle flows

Understanding the scaling of energy and particle confinement for mixed deuterium–tritium plasmas is important for fusion research as experiments ramp up from pure deuterium to 50–50 deuterium–tritium (D–T) plasmas in reactor operation. Theoretical studies of the scaling of turbulent transport with hydrogenic isotope have focused primarily on energy transport and for a single bulk ion species.^1–7 While the turbulent radial energy fluxes are generally positive in a tokamak (outward-directed from the core to the wall consistent with monotonically decreasing temperature profiles), the particle fluxes can be negative or positive. Furthermore, while the ion and electron turbulent energy fluxes can be quite different even for a pure plasma, the turbulent particle fluxes are constrained by ambipolarity such that $∑iziΓi=Γe$⁠, i.e., the charge-weighted sum of the ion particle fluxes must equal the electron particle flux. Thus, in a pure plasma with a single hydrogenic ion species $Γi=Γe$⁠. The sign of the ambipolar particle flux depends on turbulence drive^8,9 [e.g., ion temperature gradient (ITG)-driven vs trapped electron mode (TEM)-driven turbulence] and is influenced by collisions, which reduce the inward particle pinch,¹⁰ and E × B shear, which may enhance the inward particle pinch.¹¹ In a multi-ion plasma, the ion particle fluxes can differ while still maintaining ambipolarity. In early gyrokinetic simulations of ITG-driven turbulence in 50–50 D–T plasmas, it was observed that the deuterium and tritium particle fluxes were unequal.¹ This observation implied that a plasma with initially equal D and T concentrations would experience so-called flow separation, whereby the concentration would evolve to an unequal state. This buildup of one species relative to the other could potentially result in a loss in fusion performance, although no assessments of this phenomenon have been reported.

The present work generalizes the original analysis of Estrada-Mila¹ to include a more comprehensive analysis of the particle flux asymmetry in two-ion plasmas with arbitrary density fraction. In what follows, for a mixed D–T plasma, nonlinear gyrokinetic simulations are carried out in ITG and TEM-dominated turbulence regimes. Variations in both the relative ion density fractions and density gradients are explored, and conditions for restoring equal D and T fluxes are identified. To further provide an intuitive basis for the numerical results, an analytic quasilinear theory of the particle flux symmetry breaking is developed. The quasilinear theory correctly predicts the numerically computed deviation in the ion density fraction (from 50–50), or ion density gradient (from the electron gradient), required to restore equal deuterium and tritium particle fluxes.

In this section, we carry out gyrokinetic simulations of turbulent transport to explore the asymmetry between the deuterium and tritium particle fluxes in mixed D–T ion-species plasmas with varying ion density fraction and density gradient.

Here we summarize the parameters of the numerical simulations done with the gyrokinetic code CGYRO.¹² Kinetic electrons and one or two hydrogenic (z_i = 1) ion species (deuterium with $mD/me=3670$ and tritium with $mT/me=5506$⁠) are used throughout. The GA standard case parameters¹³ are used as a baseline, given by: $R0/a=3, r/a=0.5$⁠, q = 2, s = 1, $Ti=Te=T0, a/LTi=a/LTe=3$⁠, $a/Lni=a/Lne=1$⁠. Here r is the midplane minor radius, a is the midplane minor radius of the last closed flux surface, q is the safety factor, $s=(r/q)(dq/dr)$ is the magnetic shear, and $1/LTi≐−d ln Ti/dr$ and $1/Lni≐−d ln ni/dr$ are the inverse temperature and density gradient scale lengths, respectively. An unshifted Miller local circular equilibrium^14,15 is used for the geometry. Weak collisions are included with $ν¯e=(a/vtD)τee−1=0.02$⁠, where $τee−1=2πe4ne ln Λ/(me1/2Te3/2)$ is the electron collision rate and $vtD=TD/mD$ is the deuterium thermal speed. Rotation effects are not included (i.e., $γE=0, γp=0$⁠). Transverse electromagnetic fluctuations with small but finite β (⁠$βe,unit=0.05%$⁠) are included to increase the maximum stable linear time step but have a negligible effect on the transport. Particle and energy fluxes are normalized by the deuterium gyroBohm reference levels, given, respectively, by

where $ρ*D=ρsD/a, ρsD=vtD/ΩcD,unit$⁠, and $ΩcD,unit=eBunit/(mDc)$ with $Bunit(r)=(q/r)(dψ/dr)$ defined with respect to the safety factor, q, and the poloidal flux divided by $2π$⁠, ψ.¹⁵ It is also useful to define the related species gyroradius, $ρa=vta/Ωca$⁠, where $vta=Ta/ma$ and $Ωca=zaeB(θ)/(mac)$⁠. Positive (negative) particle flux implies a radially outward (inward) flow of particles from the core, which is generally unfavorable (favorable) for the confinement of fuel ions and electrons.

For the numerical resolution, the nonlinear CGYRO simulations use a radial box of length $Lx=90ρ¯$ with N_r = 128 radial modes and a binormal box of length $Ly=94ρ¯$ with $Nα=16$ complex toroidal modes. With these choices, the simulations retain complex poloidal wavenumbers,

where $kθ≐nq/r$ and n is the toroidal mode number. The radial and binormal resolutions are fixed with respect to the ion gyroradius defined by the effective mass of the ion species in the simulation, $m¯≐(nDmD+nTmT)/ne$⁠, where $ρ¯=v¯t/Ω¯c$⁠, $v¯t=TD/m¯$⁠, and $Ω¯c=eBunit/(m¯c)$⁠. Other resolution parameters are $Nθ=28$ (field line resolution), $Nξ=16$ (pitch-angle resolution), and N_u = 8 (energy resolution) with maximum energy $umax2=8$⁠. Definitions of the CGYRO numerical resolution parameters can be found in Ref. 12.

In what follows, we will vary the deuterium and tritium ion density fractions and density gradients from the GA standard case parameters, keeping the electron values fixed in a manner that preserves ambipolarity: $∑izini=ne$ and $∑i(zini/Lni)=(ne/Lne)$⁠. Three nominal values of the electron density gradient are considered, $a/Lne={0.5,1,3}$⁠. By varying the electron density gradient, a transition from ion temperature gradient (ITG)-driven turbulence to trapped electron mode (TEM)-driven turbulence occurs. This transition is illustrated in Fig. 1. In reality, the transition between the two regimes is not abrupt. However, as shown in a similar scan performed in Ref. 16, the real frequency of the various modes shifts from the ion to the electron direction as $a/Ln$ increases. We refer to the range of $0≤a/Ln≲1.5$ where the ion energy transport dominates as the ITG-dominated regime, and $a/Ln≳1.5$ where the electron energy transport is comparable to or larger than the ion energy transport as the TEM-dominated regime. In the ITG-dominated regime, the particle flux is low and the inward convective pinch driven by the temperature gradient dominates at low density gradient (⁠$a/Ln≤0.5$⁠). In the TEM-dominated regime, the particle flux is large and positive, indicating an unfavorable outward flow of particles. For this case, the neoclassical transport is negligible due to the weak collisions (⁠$|Γaneo|/ΓGBD<1×10−2$ as computed by the NEO code¹⁷) In steady-state, the transport equations require that the particle flux from turbulent and neoclassical transport balances the particle source according to

where $V′=dV/dr$ is the derivative of the flux-surface volume.¹⁸ Thus, regimes of low particle flux may be expected in the core with low particle fueling, while regimes of larger particle fluxes may be expected in the core with pellet fueling or in the edge due to a wall source.

Figure 2 shows the variation of the turbulent particle flux for deuterium, tritium, and electrons with density fraction of deuterium and tritium ranging from pure deuterium to 50–50 D–T to pure tritium at fixed (equal) density gradient scale lengths for the three cases of $a/Ln={0.5,1,3}$⁠. For each case, at any given density ratio, the deuterium and tritium fluxes are relatively centered about $Γe/2$⁠. Further, Γ_e for the two-ion-species plasma (black curve) is very nearly equal to the single-ion value with effective mass $m¯$ (stars). It has similarly been verified that the nonlinear $kθ$-spectra for the electron energy and particle fluxes are approximately invariant functions of the effective wavenumber $kθρ¯$⁠. From these results, it appears that the nonadiabatic electron dynamics (from the electron parallel motion) responds to the ion-drift timescale defined with respect to the effective ion mass $m¯$⁠.^7,16 For the ITG turbulence cases, the pure deuterium flux is less than the pure tritium flux, indicating that pure deuterium is better confined than pure tritium, while the reverse is true for the TEM case. For all 3 cases, at 50–50 D–T, there is an asymmetry between the D and T particle flows, with the deuterium particle flux more positive (outward or less inward) than the tritium particle flux, indicating that the tritium is better confined than the deuterium. For the ITG turbulence case at $a/Ln=0.5$ [Fig. 2(a)], where the electron and ion fluxes are inward, equal D and T particle fluxes are restored with higher concentration of deuterium than tritium: $nD/ne=0.72, nT/ne=0.28$⁠. But for the ITG case at $a/Ln=1$ [Fig. 2(b)], where the electron and ion fluxes are outward, equal D and T particle fluxes are restored with higher concentration of tritium than deuterium: $nD/ne=0.31, nT/ne=0.69$⁠. For the TEM case at $a/Ln=3$ [Fig. 2(c)], the fluxes are nearly linear with density fraction and the asymmetry at 50–50 D–T is very small, such that the equal flow condition is $nD/ne=0.48, nT/ne=0.52$⁠.

Figure 3 shows the turbulent particle flux as function of ion density gradient scale length. Deuterium, tritium, and electron fluxes are plotted for a 50–50 D–T plasma, varying the ion density gradient scale length around the fixed electron density gradient scale length. As in Sec. II C, both the ITG and TEM turbulence regimes are considered. In all regimes, the ion particle flux Γ_a is generally an increasing function of $a/Lna$⁠. For the ITG turbulence cases [Figs. 3(a) and 3(b)], the variation is nearly linear, consistent with diffusive-type transport, whereas the behavior is nonlinear for the TEM turbulence case [Fig. 3(c)]. In the previous section, it was shown that the ion density fraction can be varied from 50–50 D–T to restore equal D and T particle fluxes (Fig. 2). Similarly, equal flows can be restored for equal 50–50 D–T (in all turbulence regimes) by reducing the deuterium density gradient while increasing the tritium density gradient. For each case, the change in the density gradients required is quite small: 2.5–4.5%. This trend is qualitatively consistent with previous gyrokinetic simulations of 50–50 hydrogen-deuterium plasmas for Alcator C-Mod H-mode parameters, which showed that the null-flux condition could be met by increasing the deuterium (heavier species) gradient and reducing the hydrogen (lighter species) gradient.¹⁹ 

An analytical scaling of the D and T particle fluxes is developed to further understand the nonlinear simulation results of Sec. II. We write the gyrokinetic equation for the distribution function in nonorthogonal coordinates $(ψ,θ,α)$ using a spectral representation in the radial (⁠$x≐2π(r−r0)/L∈[0,2π)$⁠) and binormal (α) spatial dimensions with discrete wavenumber $k⊥$⁠.^12,15 We further employ the ballooning transformation^20,21 with extended angle $θp≐θ+2πp∈(−∞,∞)$⁠, such that

Here $f0a$ is the zeroth-order equilibrium Maxwellian distribution function, θ₀ is the ballooning angle (a linear eigenmode label), n is the toroidal mode number (also an eigenmode label), and p is the radial wavenumber. We consider the linear gyrokinetic equation in the collisionless, electrostatic limit, and use $(λ,ε)$ velocity-space coordinates, where $λ=v⊥2/(v2B)$ and $ε≐v2/(2vta2)$⁠. Without loss of generality, we consider only $θ0=0$⁠. For simplicity, we assume equal temperatures for all species, $Ti=Te=T0$⁠. We write the equation in dimensionless form by normalizing to the deuterium thermal speed, $vtD$⁠, and the midplane minor radius of the last closed flux surface, a, thus defining the normalized time as $τ≐(vtD/a)t$⁠. The gyrokinetic and Poisson equations²² become

where $H̃a=h̃a+zaJ0(γa)δϕ̃$ is the nonadiabatic distribution and $γa≐k⊥v⊥/Ωca$ is the argument of the Bessel functions. The dimensionless operators representing parallel streaming, perpendicular drift, and gradient drive are given, respectively, by

Here we define the normalized scale lengths,

The magnetic equilibrium geometry parameters Θ, G_q, $Gθ, gsin$⁠, $gcos1$⁠, and $gcos2$ are defined by Eqs. (109)–(111), (101), and (103)–(104), respectively, in Ref. 15. Using the ansatz $h̃a=ĥae−iωτ$ in Eq. (5), we can write that

where the operator $La$ is defined as

Substituting Eq. (14) into Eq. (6), the Poisson equation becomes

where the density response operators are given by

Using the artificial ordering parameter δ, we consider the following orderings for the ion species:

Physically, the ordering in Eq. (18) is expected to hold at a large aspect ratio, the finite Larmor radius ordering in Eq. (19) is valid for ion-scale (long-wavelength) turbulence regimes, and the drift ordering in Eq. (19) and the gradient drive ordering in Eq. (20) are consistent with strongly driven (more slab-like) systems where the frequency is close to diamagnetic and is larger than the magnetic drift frequency. Expanding the ion density response operators to $O(δ)$⁠, we find that

where

and $ω̂≐ω/(kθρsD)$⁠. Substituting Eq. (21) into Eq. (16), we get the following dispersion relation for the linear frequency:

where $ĉ1≐∑i(zini/ne)c1i$ and $ĉ2≐∑i(zini/ne)c2i$⁠. Solution of this matrix equation for ω is complicated as the electron response function is not diagonal in θ_p and furthermore depends on ω.

The turbulent radial particle flux is given by

where the angle brackets denote a flux-surface average. Substituting in Eq. (14), this becomes

Using Eq. (21), we can write the ion particle flux as

Simplifying by using the dispersion relation in Eq. (25) to write $1/ω2$ in terms of $1/ω$ and R_e, we find that

where $γ̂≐Im(ω̂)$ is related to the normalized linear growth rate. Using quasineutrality and assuming hydrogenic ion species (z_i = 1) with $a/LTi=a/LTe$⁠, the ion particle flux reduces to

where

Note that Δ_t and Δ_k are independent of ion species and are positive (⁠$Δt>0$ and $Δt>0$⁠) for positive temperature and density gradient scale lengths. Δ_k, which is smaller than Δ_t by order $O(kθρsD)2$⁠, arises from finite Larmor radius (FLR) corrections.

We now consider Eq. (30) for a two-ion-species plasma composed of deuterium and tritium nuclei. Analytic formulas for the D and T particle fluxes are derived in two simple limits: equal density gradients with arbitrary ion density fraction and equal ion density fraction (50–50) with arbitrary ion density gradient.

First, we consider the simplest limit of equal logarithmic density gradients, $a/Lni=a/Lne$⁠, and equal 50–50 D–T composition. In this case, Eq. (30) reduces to

This result is consistent with Eqs. (29) and (30) of Ref. 1. Because Δ_k is positive, this important result not only shows that the deuterium and tritium fluxes are asymmetric, but also that for arbitrary turbulence regime the deuterium particle flux is more positive (outward) than the tritium particle flux. Thus, this implies that for 50–50 D–T, the tritium is always better confined than the deuterium. This is consistent with the simulation results in Fig. 2 for the different turbulence regimes.

Generalizing to allow for arbitrary ion density composition, the Δ_t term in Eq. (30) vanishes due to the assumption of equal ion and electron logarithmic density gradients and we find that the deuterium and tritium turbulent fluxes are given by

In Fig. 2, for the TEM turbulence case, a nearly linear scaling of the ion fluxes with their respective density ratio was observed, while a more complex scaling was observed for the ITG turbulence cases. In the context of the analytic results in Eqs. (35)–(36), this implies that for regimes with strongly nonadiabatic electrons like the TEM case, Γ_e dominates over Δ_k, while for regimes with weakly nonadiabatic electrons like the ITG cases the two terms are comparable and a quadratic scaling is expected.

The formulas furthermore predict that the asymmetry between the fluxes is $O(kθρsD)2$ relative to the electron flux. This indicates that equal D and T fluxes (⁠$ΓD=ΓT$⁠) can be restored with a slight deviation from equal (50–50) density fraction. To lowest order, we can calculate that this occurs for

Thus, for outward electron particle flux (⁠$Γe>0$⁠), like the ITG case in Fig. 2(b) and the TEM case in Fig. 2(c), equal D and T fluxes are obtained with lower density of deuterium than tritium, while for inward electron particle flux (⁠$Γe<0$⁠), like the ITG case in Fig. 2(a), higher density of deuterium than tritium is needed. The change in the density required decreases as the electron particle flux increases, such that it becomes insignificant in the strongly nonadiabatic electron regime, such as for Fig. 2(c).

In the limit of equal local deuterium and tritium density (50–50 D–T), i.e., $nD/ne=nT/ne=1/2$ with arbitrary local density gradient, Eq. (30) reduces to

where

Here we have used the quasineutrality relation for the density gradients: $nDκDn+nTκTn=neκen$⁠. While we have written the fluxes in terms of diffusive-like and convective-like terms, it is important to note that in general Γ_e depends on the ion gradients, as do $δϕ̂$ and thus Δ_t and Δ_k. This means that D₀, V₀, and V_k are nonlinear functions of the ion density gradients and each can be positive or negative. In Fig. 3, for the ITG turbulence cases, nearly linear scaling of the ion fluxes with their respective ion density gradient was observed with equal slopes, while a more complex scaling was observed for the TEM turbulence case. In the context of the analytic results in Eqs. (39)–(40), this implies that for regimes with weakly nonadiabatic electrons like the ITG cases, the nonlinearity of these terms is weak, unlike for cases with strongly nonadiabatic electrons like the TEM turbulence case.

As noted above, the asymmetry between the fluxes is $O(kθρsD)2$ relative to the electron flux at equal density gradients. This indicates that equal deuterium and tritium fluxes (⁠$ΓD=ΓT$⁠) can be restored with a small deviation from equal ion density gradients. To lowest order, we can calculate that this occurs for

Thus, for outward electron particle flux (⁠$Γe>0$⁠), like the ITG case in Fig. 3(b) and the TEM case in Fig. 3(c), D₀ is positive, so equal D and T fluxes are obtained with flatter density gradient for deuterium and steeper density gradient for tritium, relative to the electron density gradient. For inward electron particle flux (⁠$Γe<0$⁠), this is also the case if the electron particle flux is relatively small (⁠$|Γe|<Δt$⁠). The results of Fig. 3(a) suggest that this may be the case for ITG turbulence regimes with weakly nonadiabatic electrons.

The asymmetry (or flow separation) between deuterium and tritium turbulent-driven particle fluxes in mixed D–T plasmas has been studied numerically with nonlinear gyrokinetic simulations of two-ion-species plasmas in both ITG and TEM-dominated turbulence regimes. The trends were verified with an analytic quasilinear theory of the flux symmetry breaking. The main results are as follows:

1. The electron particle flux for a mixed D–T plasma follows the flux for a single-ion species with the equivalent density-weighted average mass.

2. For equal density gradients at 50–50 D–T ion density composition, the deuterium particle flux is more positive than the tritium flux, implying that the tritium is better confined than deuterium.

3. The asymmetry between the D and T particle fluxes is $O(kθρsD)2$ relative to the electron particle flux. Thus, equal fluxes can be restored with a small deviation in ion density fraction or density gradient, such as

4. • lower density of deuterium than tritium when the electron flux is outward;

5. • higher density of deuterium than tritium when the electron flux is inward;

6. • flatter density gradient for deuterium and steeper density gradient for tritium relative to the electron density gradient when the electron flux is outward or if the electron flux is inward yet relatively small (e.g., weakly nonadiabatic electron regimes).

While these results provide insights for understanding drivers of local D–T turbulent flow separation, development of practical solutions for D–T isotope ratio is beyond the scope of this work. Modeling isotope ratio control scenarios in a reactor requires solving the steady-state transport equations with the turbulent and neoclassical fluxes and the particle sources in an integrated source-driven, core-to-edge modeling approach.^23–25 Fueling from NBI and other heating sources as well as by impurities (α particles, helium ash, and wall impurities) and neutrals complicate the transport modeling process. Control methods may ultimately involve core pellet fueling and/or edge pellet and gas fueling used for mitigation of ELMs.

This work was supported by the U.S. Department of Energy under Award Nos. DE-FG02–95ER54309 and DE-FC02–06ER54873 (Edge Simulation Laboratory), and DE-SC0017992 (AToM SciDAC-4 project). An award of computer time was provided by the INCITE program and the ALCC program. This research used resources of the Oak Ridge Leadership Computing Facility, which is an Office of Science User Facility supported under Contract No. DE-AC05–00OR22725. Computing resources were also provided by the National Energy Research Scientific Computing Center, which is an Office of Science User Facility supported under Contract No. DE-AC02–05CH11231. This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise, does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

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