Multispecies density peaking in gyrokinetic turbulence simulations of low collisionality Alcator C-Mod plasmas

Reactor designers have long assumed that flat density profiles would be an unavoidable consequence of edge fueling in large dense toroidal plasmas. Although peaked density profiles would increase fusion power¹ at constant stored energy and line average density, creating peaked densities by deep pellet injection is not feasible with present pellet injector technology, and neutral beam fueling of the core would consume too large a fraction of the electrical output to be practical in a power plant.² The neoclassical Ware pinch is too small to be effective in burning plasmas and completely absent with fully non-inductive current drive. Nature has, however, provided a turbulently driven pinch that occurs at low collisionality that is well suited to burning plasmas. Its occurrence is well documented in several tokamaks^3,4 and theoretically understood (see the comprehensive review, Ref. 4, and later work^5,6). Although the initial quasilinear study⁷ found density peaking only for extremely low collisionality, the first nonlinear simulations^8,9 based on measured parameters (for a DIII-D L-mode plasma) found that density peaking was predicted at the collisionality of the experiment.

The unexpected occurrence of low-collisionality H-mode plasmas with peaked density profiles in Alcator C-Mod¹⁰ (which had previously produced high-collisionality flat density profiles in enhanced D_α (EDA) H-modes with ion cyclotron radiofrequency heating (ICRH)) confirmed the empirical association of density peaking and low collisionality observed in AUG and JET, and it removed a competing association with “Greenwald density fraction” that had no theoretical link to density peaking (AUG and JET data had high covariance of collisionality and n/n_G). It also extended the phenomenon to plasmas with approximately equal ion and electron temperatures, and in some cases with T_e > T_i. The density profiles evolve more quickly than the time scale of the Ware pinch so a turbulent pinch is required.¹⁰ 

The first nonlinear turbulence simulations of these C-Mod plasmas¹¹ revealed that the particle flux has a bidirectional k_yρ_s spectrum for conditions in or near the low-collisionality density peaking regime: outward flux for low k_yρ_s and inward for higher k_yρ_s modes, with a transition in the vicinity of 0.2 < k_yρ_s < 0.4. The crossover point moves to higher k_yρ_s as either the collisionality or a/L_n is increased. When these controlling parameters are sufficiently large, all the modes contribute to outward flux, and when they are sufficiently small, all modes produce an inward flux. The first quasilinear study of density peaking⁷ had no similar nonlinear guide, adopted a low k_yρ_s and found that density peaking occurred only for extremely low collisionalities. Later nonlinear simulations for AUG conditions^12,13 reproduced the bidirectional features of the k_yρ_s spectrum found in the C-Mod simulations, and detailed theoretical analysis^12,14 led to better choices for k_yρ_s and successful quasilinear predictions of density peaking for JET as well.¹³ 

Strictly empirical predictions of density peaking $ne(0)/〈ne〉∼1.5$ in ITER were confirmed by nonlinear calculations for parameters similar to those expected in ITER,¹² and this would produce a 30% enhancement in fusion power (at the design β_N and $n¯e$⁠) and the same enhancement of Q at fixed heating power (see Ref. 4 and citations therein). Other benefits of density peaking would be higher bootstrap current, transport of edge-fueled D and T to the plasma center, lower edge density should increase the pedestal temperature at the edge localized mode (ELM) stability limit and reduce the chance of high-density disruptions at high Greenwald density fraction. However, density peaking could produce unwelcome changes, such as lower β_N due to the higher pressure peaking, reduced radiative cooling from the lower density near the plasma edge would tend to raise divertor heat loads, and impurity density peaking might lead to serious central radiative losses.

Broadening our understanding of density peaking should facilitate planning to deal with these impacts. Tritium fueling and retention are critical issues for ITER and other burning plasma tokamaks, so it is important to learn the isotope dependence of the fuel transport. A pioneering investigation¹⁵ of isotope dependence in burning plasmas with the GYRO Standard Case parameters found slightly more density peaking for tritium than deuterium. Impurity transport sometimes exhibits a charge dependence so it is important to find the expected density peaking for a wide range of species. In some AUG experiments, both experimental and quasilinear theoretical analyses were in qualitative agreement that boron peaking is consistently less than for electrons.⁶ However, one quasilinear study¹⁶ found that peaking increases with impurity charge up to Z ∼ 16 (assuming A = 2Z), but then drops for larger Z.

While analytic quasilinear formulations enable detailed accounting of the importance of different theoretical elements, they involve assumptions that should be validated by comparison with more complete nonlinear simulations. Turbulence saturation is governed by nonlinear processes, and these also determine the relative importance of modes with different k_yρ_s. Quasilinear estimates of particle flux can be particularly unreliable because oppositely directed fluxes^12–14 commonly occur for different k_yρ_s so their relative contributions must be known to obtain the correct net flux. The quasilinear estimates of Ref. 12, for instance, predict a crossover from outward to inward fluxes at k_yρ_s that are significantly higher than what is seen in nonlinear simulations for the same conditions. Nonlinear simulations should closely guide quasilinear estimates, which then enable efficient parametric scans.

We explore whether density peaking is predicted by gyrokinetic turbulence simulations of low collisionality Alcator C-Mod plasmas, and explore the isotope and species dependence of hydrogenic density peaking. Impurity transport is dominated by turbulent processes in L-mode plasmas in C-Mod,^17,18 so we also examine the species dependence of impurity density peaking. These studies are based on nonlinear turbulence simulations because quasilinear techniques can be unreliable unless they are founded on relevant nonlinear simulations,^5,12,19,20 and there are no prior nonlinear studies of isotope and species dependence of low-collisionality density peaking to use as a guide.

While neoclassical transport can be important for impurities in some tokamaks,^21,22 studies for L-mode plasmas in C-Mod find that turbulent transport dominates.^17,18,23 Experimental analysis of impurity transport in low-collisionality H- and I-mode plasmas in C-Mod has not been carried out, but particle confinement of the majority species in these regimes is indistinguishable from L-mode plasmas^24,25 so we have not evaluated neoclassical predictions here.

The density peaking predictions are based on calculations by GYRO, a time-dependent, nonlinear gyrokinetic drift-wave turbulence simulation code.^26–28 GYRO and similar codes have been verified by an extensive series of benchmarks^29–35 with the GS2,^36,37 GEM,³⁸ GENE,^39,40 and PG3EQ (Ref. 41) codes. Verification, based primarily on benchmark calculations from independently developed gyrokinetic codes, assesses how well a code solves the equations on which it is based, while validation assesses whether a code's mathematical model provides an accurate description of the relevant experimental measurements.

GYRO validation studies based on DIII-D plasmas have dealt with ρ_* scaling of heat fluxes in the L- and H-mode regimes,^27,42 and together with density fluctuations in both L-mode^43–45 and H-mode plasmas.⁴⁶ Comparisons of predicted and measured heat fluxes and electron temperature fluctuations,^43,47 as well as the cross-phase between temperature and density fluctuations^48,49 have been compared in L-mode plasmas, including scans^44,45 of elongation and T_e/T_i. Plasmas with strongly modulated T_e gradient have been used in validation studies of the electron heat flux and T_e fluctuations,^50,51 as well as the electron stiffness.⁵² 

In the outer region of DIII-D L-mode plasmas, GYRO simulations often predict lower heat transport and fluctuation levels than are measured,^43,45 and this validation problem is referred to as a “shortfall.” A similar shortfall is evident in GEM simulations⁵³ of DIII-D (but not C-Mod). The trapped gyro-Landau-fluid (TGLF) reduced transport model that is based on GYRO simulations exhibits a large shortfall, too.^45,54 However, a similar shortfall has not been reported for NSTX,⁵⁵ C-Mod,⁵⁶ or AUG,⁵⁷ and GENE simulations of the original DIII-D shortfall discharge produce a substantially higher level of electron and ion heat transport than seen in GYRO and GEM simulations. The source of these code to code differences is under intense study,⁵⁸ but is not yet resolved.

GYRO validation studies based on Alcator C-Mod experiments have compared heat fluxes and density fluctuations in ohmically heated plasmas,^59,60 heat fluxes in RF-heated L-mode plasmas,^56,61,62 impurity transport in L-mode,^17,18,63 and multichannel L-mode studies.^64,65

The experimental conditions used for our turbulence simulations are described in Sec. II, and the basic simulation inputs, resolution choices, and initial simulations are discussed in Sec. III. We show that the mode spectrum of the particle flux is bidirectional in Sec. IV and that the theoretically expected density peaking diminishes at higher collisionality. The initial technique for estimating the density peaking predicted by GYRO is described in Sec. V, and a refined technique is described and used in Sec. VI to predict that hydrogenic species will have a small isotope dependence, and that the degree of impurity density peaking will increase with impurity charge.

A regime of low-collisionality H-modes in Alcator C-mod was discovered serendipitously during an experiment designed to produce type-I ELMs.¹⁰ Type-I ELMs had not previously occurred in C-Mod, but with the JFT-2M plasma boundary shape, ELMs do occur and their frequent expulsion of plasma from the edge and/or the deeper strike point position and altered recycling patterns may be responsible for lowering the average density. The resulting density profile shapes are notably unlike the typical EDA H-mode plasmas in C-Mod, which have very flat density profiles.^10,66 The density profile is not strongly centrally peaked (as it is in C-Mod ITB plasmas), only the outer 50–60% of the minor radius has a significant gradient. The degree of peaking is correlated with the collisionality, confirming the empirical association of density peaking and collisionality originally found in AUG and JET.¹⁰ In later experimental campaigns, density peaking has occurred in H-mode plasmas with the standard C-Mod shape, as well as low-density I-mode plasmas^24,25 and ohmic H-modes.

The turbulence simulations discussed here are based on data measured at 1.0 s during an ELMing low-collisionality H-mode phase of shot 1080213027, formed with the JFT-2M shape. The major macroscopic characteristics are I_p = 850 kA, B_o = 5.4 T, a = 0.22 m, R_o = 0.67 m, κ = 1.4, δ = 0.4, q₉₅ = 3.7, P_rf = 3.5 MW, P_rad = 1.6 MW; $n¯e=2×1020/m3$⁠, Z_eff ∼ 1.7, and n_H/n_D ∼ 0.06. Large sawteeth, with δT_e/T_e ∼ 25%, dominate central transport, so this region is not included in the turbulence simulations.

The modeled plasma has a line average density about 40% lower than an EDA H-mode in C-Mod with similar plasma current and heating power. As shown in Fig. 1, the fits to the Thomson scattering density profile are moderately peaked, as are the two predicted density profiles based on GYRO simulations. The “initial prediction” is based on a matched pair of simulations, one with the experimental density profile and the second with the density gradient reduced by 30%. The “twin prediction” is based on a single simulation with different a/L_n assigned to “twins” of each ion species; details are given in Sec. VI.

A recently developed Gaussian process regression (GPR) technique⁶⁷ has been used to obtain statistically rigorous fits and uncertainty estimates for both the density and the density gradient parameter, a/L_ne. Shown in Fig. 1 as the dotted line (with ±2σ “uncertainty bars”), this new fit is based on 150 ms of Thomson data during the H-mode phase, and possesses a more gradual variation of the gradient than the least-squares spline fit, but with essentially the same variation in density between 0 ≤ r/a ≤ 0.9 as seen in the spline fit.

The ion temperature is collisionally coupled to the electron temperature across the entire profile in high-density EDA H-mode plasmas. This is not the case in low density H-mode and I-mode plasmas in C-Mod where the central T_e can be 50% higher than T_i, or in AUG and JET plasmas where central T_i is usually much higher than T_e. In the discharge studied here, the two temperatures happen to be very similar across most of the profile (although this is not forced by collisional coupling). The results from all three tokamaks demonstrates that neither equal nor unequal temperatures are correlated with density peaking or its absence.

The first generation of GYRO simulations uses measured ion temperatures provided by a full profile diagnostic.^68,69 As shown in Fig. 2, the gradient scale length varies considerably, and this produces neighboring regions of strong and weak turbulence in the simulations. The results presented in Secs. III–VI are based on modified ion temperature profiles (see Fig. 2) constructed from (1) a flattened gradient scale length that smoothly joins the ion temperature measurements in the peripheral and central regions or (2) a minor adjustment to the flattened gradient that produces a better match between the profile shape of the heat flux predicted by GYRO and the experimental analysis by TRANSP. While the modified profiles are plausible, the simulations cannot be used for validation of the gyrokinetic model for turbulent transport because the gradient changes are large and, in any case, it was later recognized that the T_i measurements for the discharge studied here were made with the initial vignetted diagnostic view.⁶⁹ While this vignetting was removed in later campaigns, it is not possible to construct reliable ion temperature profiles for the discharge studied here. We do expect, however, that the simulation results are representative of this type of plasma because the results are qualitatively similar in a number of simulations using different shapes for the ion temperature profile.

The results described in Secs. III–VI are taken from nonlinear turbulence calculations made by the continuum (or Eulerian) gyrokinetic code GYRO,^27,28 which solves the gyrokinetic Vlasov-Maxwell equations. The simulations include non-circular shaping of magnetic flux surfaces through use of the Miller local-equilibrium model.^70,71 Multiple ion species are modeled with the full gyrokinetic equation, but electron dynamics are based on the drift-kinetic equation, including the realistic mass ratio $mD/me=60$⁠, and a Lorentz collision operator is used to model pitch-angle scattering of electrons due to electron-electron and electron-ion collisions.

For the plasma conditions in these simulations, the most unstable modes are driven mostly by the ion temperature gradient (ITG), but non-adiabatic electron effects provide very significant additional drive so the turbulence is best characterized as a hybrid of ITG and TEM, as first described in Refs. 72 and 73. Consequently, the simulations employ non-adiabatic treatments of all ions and the electrons.

The simulations are based on measured profiles of electron temperature and density (the density is modified in some cases noted below), and modified ion temperature profiles shown in Fig. 2 and noted where used. The “flat gradient” ion temperature profile produced a broad region of turbulent transport, but peak total heat flux was over 7 MW so the inverse ion temperature gradient scale length was reduced by ∼20% in some of the simulations to reduce the peak total heat flux to 3 MW. While this magnitude is appropriate for the experimental conditions, the TRANSP power balance analysis finds that 2/3 of the power is carried in the electron channel, but in the GYRO simulations, 2/3 of the power is in the ion channel. Experience with simulations of a number of C-Mod plasmas indicates that matching the TRANSP ion/electron power mix can be achieved only with a large increase in the electron temperature gradient (ETG). There is no reason to expect the two independently calibrated T_e diagnostics have such serious errors so we continued to use the measured electron temperature data.

Computational expense has prevented us from undertaking multi-scale simulations that include short wavelength turbulence, such as ETG modes. For some C-Mod plasmas, the inclusion of ETG modes can raise the electron heat transport and allow matching of both the ion and electron heat fluxes.⁶² In some multi-scale simulations, the presence of short wavelength modes does not qualitatively change the results for ion transport because that is dominated by strongly-driven long-wavelength ITG turbulence.^74–76 However, when the ion drive is sufficiently weak, the inclusion of ETG modes may significantly alter the long wavelength ion heat flux⁶² (and, presumably, the particle flux) so a significant role for ETG modes cannot be ruled out for the conditions examined here without further investigation.

Electromagnetic effects are not expected to be strong for these plasmas because the normalized pressure, β, does not approach the ideal-MHD ballooning limit. Effects of sheared rotation and ExB shear are not included, based upon their modest effect in simulations of similar plasmas with rotation measurements^17,61 and typically low ExB shearing rate compared to the linear growth rates.^56,61–63

The simulations use a “global” domain, which includes radial variations of the magnetic equilibrium's geometric parameters as well as density and temperature, and the gradients of these quantities. While turbulence simulations based on “local” flux-tube geometry (with no radial variations) are often similar to global simulations^30,51 for ρ_* values that are typical for C-Mod, some local simulations^51,52 produce “unphysical streamer-type eddies with radial correlation lengths on the order of the radial domain size,” although the corresponding global simulations have “physically relevant radial correlation lengths.” The full extent of the binormal domain is L_y = 85ρ_s, where the ion sound speed gyroradius, ρ_s is defined as ρ_s = c_s/Ω_ci, $cs=Te/mi$ and Ω_ci = eB/cm_i, while the radial domain encompasses 0.3 ≤ r/a ≤ 0.8, with L_x = 124ρ_s. The 312 radial grid points used produce a radial resolution of dr = 0.40ρ_s.

The simulations described in this paper were run in fixed-gradient mode, with unchanging profiles of density and temperature. The simulations employ fixed (or Dirichlet) boundary conditions with vanishing fluctuations at the outer boundary of two buffer regions, each of width ∼8–15 ρ_s, which lie outside the physical domain that extends over 0.34 ≤ r/a ≤ 0.74. This method permits the fluctuation amplitudes to find their own level within the physical domain by applying the boundary condition at remote locations. The recommended buffer width is 8ρ_s, and a convergence study for multi-scale simulations of ITG/ETG turbulence in a C-Mod L-mode plasma found that buffer zones of 6ρ_s were sufficient.⁶² Insufficient buffer width or damping can allow fluctuations to suppress the temperature gradient over large scales within the physical domain, and consequently produce incorrect estimates of the turbulence level and transport,⁷⁷ but in the simulations reported here, the fluctuations do not cause significant large-scale modifications of the temperature gradients. The turbulent transport fluxes are prevented from making large-scale modifications to the equilibrium profiles (also known as profile relaxation) by an adaptive source model that damps the long-radial-wavelength parts of the n = 0 components of the fluctuations, as described in Ref. 26. Short-scale variations, such as zonal flows are not modified by this procedure.

The parallel spatial grid along the field line is defined by equal intervals of orbit-time,²⁶ with 14 intervals for passing orbits and twice that for trapped orbits; this exceeds the recommended GYRO resolution for standard aspect-ratio tokamaks, and was chosen to tailor the domain decomposition mapping to the available computer architecture. The velocity-space grid is composed of eight energies and eight pitch angles for each direction of parallel velocity for a total of 128 points, the recommended GYRO resolution for standard aspect-ratio tokamaks, to provide well converged results.

The nonlinear electrostatic simulations described here include a total of 16 long-wavelength toroidal modes, the uniformly spaced binormal wavenumbers span 0 ≤ k_yρ_s ≤ 1.1, with Δk_yρ_s = 0.074. A convergence test with 16 toroidal modes, but the maximum k_yρ_s raised to 1.6, had the same average heat fluxes to within half the uncertainty of each flux. The “initial prediction” in Fig. 1 was used as the input density profile, so the absolute magnitude of the average particle fluxes is 0.05 gyroBohm units or smaller in the convergence test and the base case, and each species' particle flux changed by less than 0.01 gyroBohm units. The particle flux k_yρ_s spectra have more modes with inward flux, but the location of the crossover from outward to inward flux moves to slightly higher k_yρ_s, and the magnitude of the low k_yρ_s peak inward flux is increased so the changes offset each other. Variations in density profile were not also carried out (so null-flux roots could not be derived), but raising the maximum k_yρ_s to 1.6 would not be expected to make a significant difference in the predicted density peaking since the changes to the particle flux k_yρ_s spectra are quite small when compared with the changes shown in Fig. 3 which are caused by significant changes in the density gradient. Time integration combines a second-order split implicit/explicit treatment of the linear physics with a fourth order Runge-Kutta treatment of the nonlinear terms; typical time steps are 0.005–0.0029 a/c_s.

Simulation results are typically averaged over several hundred time units (a/c_s) during the steady-state phase of saturated turbulence. Statistical uncertainties arising from the time-averaging are typically much less than 10%, but these are overshadowed by systematic uncertainties arising from inaccurate experimental inputs, particularly for the most important driving gradient, a/L_Ti.

Here, we describe results from turbulence simulations with one ion (deuterium) and electrons that were used to establish that turbulent density peaking is expected in low collisionality C-Mod plasmas, that it grows as collisionality is reduced, and that revealed the bidirectional structure of the k_yρ_s spectrum of the particle flux. The bidirectional shape of the particle flux wavenumber spectrum explains why the initial gyrokinetic study of density peaking in AUG and JET (Ref. 7) found the boundary of the low-collisionality peaking regime occurred at a collisionality that was much lower than in the experiments. Most of the initial C-Mod simulations include only two species, electrons and deuterium, with Z_eff = 1.0, but a few with boron, raising Z_eff to 1.5 or 2, demonstrated that the presence of diluting impurities did not qualitatively change the features of density peaking. The “flat gradient” ion temperature profile shown in Fig. 2 and the measured electron density profile were used in the simulations discussed in this section. With other parameters taken as measured, this produced a peak total heat flux over 7 MW so the inverse ion temperature gradient scale length was reduced by ∼20% in some of the simulations. In spite of the large variations in heat flux across simulations, we found that the density peaking and the qualitative features of the particle-flux mode spectra were essentially unchanged.

The most important result from the initial simulations¹¹ is the recognition that in the low-collisionality density peaking regime, the k_yρ_s spectrum of the particle flux is bidirectional: the flux is outward for low k_yρ_s and inward for higher k_yρ_s modes. As shown in Fig. 3, the location of the crossover point depends on collisionality and density gradient. Compare spectra with the experimental values for these parameters to those with collisionality reduced to a third of the experimental level: with lower collisionality more modes contribute to the inward flux, the magnitude of the inward flux grows at high k_yρ_s, and the net inward flux is much greater. Conversely, for sufficiently high collisionality, no modes have significant inward flux. The collisionless flux is inward even for low k_yρ_s if a/L_n is not large.

The flux's collisionality dependence is quite complex in velocity space since it can reverse sign depending on the value of k_yρ_s.¹² The quasilinear analysis of Angioni et al.¹² shows that the most dramatic collisionality dependence can be traced to ITG modes. Another generic result is that low energy particles generally contribute to inward flux, while higher energy particles are responsible for the outward flux, and collisionality affects the inward flux more than the outward. The most detailed quasilinear analysis of the sources of density peaking¹⁴ particularly emphasizes the contribution to the inward flux by modes with low real frequency, which causes density peaking to be maximized near the ITG-TEM transition. Experimentally, density peaking is strongest in the lowest collisionality C-Mod plasmas, with high power ICRH and high T_e, where linear stability calculations find that TEM drive is slightly dominant or nearly as important as ITG drive, while in most plasmas, ITG drive is dominant. The dependence of the flux on the density gradient has not been examined in as much detail as the collisionality dependence, but the direct Fick's law term, Γ ∝ D∇n, may dominate when the density gradient is not close to the value needed to stabilize ITG turbulence or close to the threshold for destabilizing TEM turbulence (neither of these circumstances applies to the C-Mod plasma studied here).

The dependence on density gradient is similar to that of collisionality: as a/L_n rises, the magnitude of outward flux grows, the crossover point moves to higher k_yρ_s, and the inward flux is diminished for the modes with high k_yρ_s. When either of these controlling parameters are sufficiently large, all the modes contribute only outward flux, and when both are sufficiently small, all modes produce an inward flux. Nonlinear simulations for AUG conditions^12,13 also display these bidirectional features of the k_yρ_s spectrum and the dependences on collisionality and density gradient. As collisionality is lowered, the turbulently driven inward flux will raise the core density and its gradient until the rising gradient offsets the effect of the lower collisionality.

In the absence of a nonlinear simulation based on peaked density plasma conditions, the first quasilinear study of density peaking⁷ was guided by a comparison of quasilinear and nonlinear fluxes for collisionless TEM turbulence.^78,79 The recommended “rule of thumb” was to maximize γ/(k_yρ_s),² which occurred at k_yρ_s ∼ 0.15, but this part of the spectrum is responsible for outward flux in nonlinear turbulence simulations based on peaked-density plasmas in both C-Mod and AUG.^12,13 Consequently, density peaking was found only for extremely low collisionalities “at least one order of magnitude lower” than in the experiments studied. After the C-Mod and subsequent AUG nonlinear simulations were available, revised guidance^12,14 for choosing k_yρ_s led to successful quasilinear predictions of density peaking for JET.¹³ That work used modes that maximize the growth rate, (k_yρ_s ∼ 0.4), and they fall within the range that produces inward flux in the nonlinear turbulence simulations based on peaked-density plasma conditions in both C-Mod and AUG. As shown in Fig. 3, the quasilinear flux at k_yρ_s ∼ 0.3 – 0.4 is a rough guide to the sign of the total flux so this quasilinear guide would probably work for C-Mod plasmas. The flux at significantly lower or higher values of k_yρ_s, however, would not closely track the sign of the total flux of the C-Mod simulations reported here or the AUG spectra in Ref. 12.

Changing the collisionality, density gradient parameter, or ion temperature gradient parameter affects the overall level of heat transport, but the qualitative changes of the mode spectra shown in Fig. 3 are not related to the magnitude of the heat flux, per se. As shown in Fig. 4, the heat flux in all the simulations with 20% lower ion temperature gradient parameter is lower than in those with the higher gradient, but the particle flux spectra in Fig. 3 exhibit qualitatively the same dependence on collisionality and density gradient within both the high and low heat flux groups. As discussed in another context in Ref. 80, the local peaks in Fig. 4 arise from fluctuating zonal flows that do not time average away completely near lowest-order rational q surfaces, and these can be ignored in the present context.

In a plasma with no core fueling like C-Mod, which is highly opaque to neutrals and has only RF auxiliary heating, the density profile shape will achieve its stationary state when the particle flux vanishes everywhere in the source-free region. The stationary condition occurs when the inward flux arising from the pinch is offset by diffusion down the density gradient that the pinch has produced, and in theoretical studies, a convenient local measure of the peaking is a dimensionless ratio of the pinch velocity to the diffusivity, which for null flux obeys the identity $−aV/D=−(a/n)(dn/dr)=a/Ln$⁠, where a is the minor radius, V is the pinch velocity (negative is inward), D is the diffusivity, and L_n is the density gradient scale length. In this work, the local value of a/L_n that nulls the particle flux will constitute the measure of the predicted density peaking. The predicted density profiles are obtained by integrating the null-flux a/L_n profile, and normalizing to a standard boundary density to facilitate comparisons.

One of the advantages of quasilinear studies is the ability to calculate directly the individual diffusivity and pinch terms, while nonlinear simulations simply produce a total particle flux. Finding the density gradient scale length that nulls the particle flux with nonlinear simulations is an iterative process, and the initial procedure is described in this section. A more efficient approach was developed later for the isotope and impurity studies discussed in Sec. VI, but it can be used to determine the density peaking of a single species as well.

The procedure for determining the density profile that is consistent with the simulated turbulent particle transport is based on the expression for the particle flux

recast in a normalized expression below as a linear function of a/L_n:

The particle flux produced by two simulations with different input density gradients will allow us to solve for D and V if the difference in the particle fluxes is caused mainly by the different density gradients rather than differences in D and V in the two simulations. The validity of the decomposition is then confirmed by using the density profile “predicted” by the null-flux solution as the input to a simulation that actually produces a particle flux close to zero. The results from the last simulations used in such an iterative process are shown in Figs. 5–7. The gyroBohm-normalized particle flux is ≪1 everywhere and the null-flux a/L_n differs from the input values by less than 0.1 over most of the domain, with both positive and negative differences (Fig. 6), so the density profile that would produce exactly zero flux is very close to the input density profile. The “tuned gradient” ion temperature profile shown in Fig. 2 was used in the simulations discussed in this section, but a/L_Ti was increased by 10%; the input electron density profile is labeled “initial prediction” in Fig. 1.

If D and V vary slowly with changes in a/L_n at fixed r/a, then two values of a/L_n and their corresponding particle fluxes are sufficient to solve for the root that will produce Γ = 0 at each radial location. The density profile shape implicit in the solution has an arbitrary normalization that we determine by specifying the density at the outer boundary of the domain

where r/a = b at the outer boundary of the domain.

In order to minimize the changes in D and V, it is desirable to apply small changes to the input a/L_n, but then the changes in the output particle fluxes are not large and the flux difference between the two simulations will be more easily affected by the stochastic variations inherent in turbulence simulations, so the estimate of the a/L_n needed to null the flux will be noisy and less accurate. We adopted a 20% reduction of the input a/L_n, which is sufficient to bracket Γ = 0 for r/a > 0.6 (see Fig. 5), so the root there is an interpolation between the two simulations. In the inner half of the domain, the density gradients are lower so the change in flux is much smaller. Nevertheless, the a/L_n root is reasonably well behaved (see Fig. 6); where the fluxes cross and the root diverges (at r/a = 0.41) it nevertheless integrates to only a small jump in density that does not overwhelm the broadly distributed density peaking (see Fig. 7). The small gradients in the inner half of the domain (and the correspondingly small fluxes) prevent the relatively large extrapolation needed to reach the root from affecting the overall peaking, most of which occurs in the outer half of the domain. Although the density gradient in the outermost part of the domain is not an exact match to the input density (see Fig. 1 to compare with two fits to the measured electron density profile), the overall changes in density and its gradient are rather small so the iterations ceased at this point. While the null-flux predicted gradients intersect the fits to the measured density, the predicted radial variation of the gradient is not similar to that of either fit. The unavailability of a reliable ion temperature profile prevents a validation study here, so we can only conclude that the turbulence simulations predict a moderate degree of density peaking that is not greatly different from the measured overall trend but the predicted shape differs locally from the measured density.

Having introduced the two-simulation version of the null-flux method, we next describe a more robust algorithm that is used to determine the dependence of the density peaking on the mix of hydrogenic isotopes and the variation of peaking among impurity species.

The technique used for determining D and V for impurity transport was modified to enable an estimate of the density gradient scale length that produces null flux from a single simulation by using “twinned” species with identical charge and mass, but different density profiles with steeper and shallower gradients than the electron density. For hydrogenic isotope studies, two D and two H species were used; the impurity studies included one D, but two copies each of two impurities (helium and boron, for instance).

The sum of the densities for the two D species in such a simulation has the same shape as the electron density, but their individual density gradient scale lengths are systematically smaller and larger than that of the electrons (see Figures 8(a) and 8(c)). Such a simulation has the same turbulence level and total D flux as a simulation with a single D species and the same total density profile. A single simulation with twinned species is sufficient to determine accurately the a/L_n that nulls the particle flux if D and V are independent of the a/L_n of the background plasma, but further iterations of this method using a sequence of modified density profiles would typically be needed to demonstrate full convergence to the null flux solution.

A turbulence simulation with a new input density gradient may have a significantly different level of turbulence and, crucially, the relative contributions of different ITG and TEM terms may be affected so the prediction for the a/L_n that nulls the flux in the second simulation might be significantly different. The twinned simulations shown here began with a density profile that had been designed to produce null flux (using the method outlined in Sec. V) so it tells us nothing about how rapidly this procedure will converge in general. Nevertheless, the twin technique greatly improved the robustness of the isotope scaling studies, which had begun with the two-simulation approach, and the inevitable small flux differences, due to slightly differing overall turbulence levels from one run to the next, were difficult to disentangle from the flux differences caused by the changes in input density gradient scale lengths. The great strength of the twin method is that the fluxes associated with the two different a/L_n are based on the same turbulent eddies experienced by the twin species simultaneously; so even significant stochastic variations in overall turbulence level tend to have a weak effect on the predicted a/L_n required to null the flux.

With twins of both the D and H species, we have determined the null flux conditions for each species. The ion temperature used in the simulations shown in Figure 8 lies between the “flat gradient” and “tuned gradient” profiles shown in Fig. 2(a), while a/L_Ti is similar to the “tuned gradient” in Fig. 2(b). Figure 8(b) shows that the D particle fluxes for the twins' density profile shapes are quite different so the null flux is obtained by interpolation everywhere. The input a/L_n for each D twin is shown in Fig. 8(c), as well as the a/L_n expected to produce null flux for D and H individually, and the a/L_n of the input electron density. The close similarity of the two null flux solutions and the input background density indicates (1) that the two-simulation procedure followed in Sec. V was successful, and (2) the hydrogenic isotope dependence of density peaking is very weak. The input and null flux profiles for electron density are shown in Fig. 1 as the “initial prediction” and “twin prediction.”

The density peaking is quite similar for D and H, and the small difference is not much affected by varying the isotope mix or the overall turbulence level. The null-flux a/L_n profiles for both D and H are shown in Fig. 9 for H:D mixes of 50:50 (used in Fig. 8) and 40:60. The predicted peaking and small isotope difference are also robust to variations in turbulence level. Other simulations with 10% higher ion temperature gradient had double the turbulent heat transport (see Fig. 10), but the predicted density peaking of each species varied by less than 5%. Simulations with 3.4% boron (and 20% higher a/L_Ti to offset the hydrogenic dilution) also gave very similar peaking for both the hydrogenic species.

The simulations used to estimate the density peaking of impurities include kinetic treatment of electrons, deuterium, and two twin pairs of impurities. The set of impurities—helium, boron, neon, argon, iron, and molybdenum—covers the range of Z for which data are or might become available. The charge of the impurity is taken to be the nuclear charge for all but molybdenum, for which we used Z = 32, which is appropriate for C-Mod core T_e. In most of the runs, we assumed a trace level of impurity because large amounts of anything beyond helium will promptly lead to radiative collapse of the discharge. Some simulations included helium with a density fraction 2.5%, but this had very little impact on the particle fluxes of either the majority D or trace boron present in simulation. The ion temperature and a/L_Ti used in these simulations discussed lies between the “flat gradient” and “tuned gradient” profiles shown in Fig. 2, and the input electron density profile is labeled “initial prediction” in Fig. 1.

The isotope scaling of impurity peaking predicted by GYRO is illustrated in Fig. 11. The dramatic peaking of some species for r/a < 0.45 may be unrepresentative of this regime as it is not seen in the extensive quasilinear work and the turbulence in this part of the domain is very close to threshold. The helium density profile is less peaked than the electron density, and is nearly flat. The boron density peaking is more or less comparable to the electron density, and higher Z species have increasingly more peaked density profiles as far as iron. The molybdenum density is not significantly more peaked, however. The Z dependence up to argon is similar to the quasilinear estimates of Ref. 16 with A = 2Z, but the sustained peaking for iron and molybdenum in Fig. 11 does not reflect the decrease previously found for Z greater than argon nor the decrease found for tungsten with a high ratio of A/Z = 4 (in our molybdenum calculations A/Z = 3). These differences may have less to do with the simplifications of the quasilinear approach than with the sensitive parametric dependence of impurity transport that is produced by the multiple flux contributions of opposing signs that are described in Ref. 16 and references cited therein.

The collisionality dependence of impurity density peaking was investigated in a second simulation with three times the experimental collision rate used initially. As seen in Fig. 12, the density peaking is slightly lower for r/a > 0.5 for all species, but the Z dependence is unchanged: increasing peakedness through iron, but no additional increment for molybdenum. The increased collisionality causes the plasma to move even closer to threshold for r/a < 0.45 and that may be why the medium to high-Z gradients are even steeper there than in Fig. 11.

It is important to remember that in gyroBohm normalized units, ITER will also be close to threshold, so studies of impurity peaking intended to be ITER-relevant should take care to use parameters that produce the expected level of heat transport. The present results should not be assumed even qualitatively valid for ITER conditions, but they do indicate that very strong impurity density peaking is not expected in the outer half of low-collisionality C-Mod plasmas (we have low confidence in the results for smaller radii where the turbulent heat flux is approaching zero). The density peaking of hydrogenic species has not exhibited a similar sensitivity in regions near threshold, but ITER-focused studies, such as Ref. 81, should be done for a range of impurities.

Experimental measurements have shown that boron is not strongly peaked in a number of C-Mod plasma regimes, but none is a close match to the plasma simulated here.⁸² Impurity peaking data are available for laser blow-off experiments in similar discharges, but analysis has not been completed.

Density peaking predicted for low-collisionality plasmas⁴ in ITER and other burning plasma tokamaks would raise fusion power, raise the bootstrap current, and transport edge-fueled D and T to the plasma center. However, density peaking could also produce unwelcome changes, such as lower β_N due to the higher pressure peaking, reduced radiative cooling from the lower density near the plasma edge would tend to raise divertor heat loads, and impurity density peaking might lead to serious central radiative losses. Predictive capability is needed in order to guide hardware design and operational planning as needed.

We have used the turbulence code GYRO to estimate the theoretically expected density peaking in a low-collisionality H-mode plasma in Alcator C-Mod. Modest electron density peaking is predicted by these simulations—and observed in the experiments—and higher collisionality is predicted to reduce peaking, as is observed in higher density EDA H-mode plasmas. A quantitative validation exercise is not possible here because the ion temperature estimate was hampered by a vignetted view,⁶⁹ so we do not have a reliable measurement of the most important source of microinstability drive. There were no fluctuation measurements for the modeled plasmas either so even a qualitative comparison is not possible. In addition to producing the qualitative prediction described above, the simulations are useful in advancing our understanding of density peaking, as detailed next.

The simulations produce a bidirectional k_yρ_s spectrum of the particle flux that necessitates a cautious approach to quasilinear estimates of the particle flux, because it is essential to use the correct relative weighting of the contributions from turbulent modes of differing length scales. Prior theoretical analyses^12,14 also indicate that even for a fixed k_yρ_s, there is very significant cancellation of inward and outward fluxes arising from different mechanisms. This means that quasilinear estimates should be well grounded by nonlinear simulations, and the first such work involving density peaking of multiple ion species is reported here.

We find that density peaking is expected for all ion species in a typical low-collisionality Alcator C-Mod plasma. The expected isotope dependence of density peaking for hydrogenic species is too weak to be detectable with current diagnostic technology. Impurity density peaking is predicted to increase with impurity Z (at least to Z ∼ 25), and is probably strong enough to be measured for Z greater than ∼ 10. Experimentally, very weak density peaking is found for tungsten for r/a < 0.6 in low-collisionality H-mode plasmas in C-Mod,⁸³ but peaking similar to that shown in Fig. 11 is seen for r/a > 0.6. Additional experimental data on transport of high-Z impurities do exist for I-mode and H-mode plasmas in the low-collisionality regime, but they have not yet been analyzed. Results on high Z impurity peaking are of some concern in burning plasmas so the analysis, and perhaps further experiments, should be carried out.

In general, predictions of density peaking agree with experimental trends,^6,12–14 and the experimental scalings as well as initial gyrokinetic calculations for ITER indicate that it will have density peaking of $ne(0)/〈ne〉∼1.5$ with ∼30% enhancement of fusion power.⁴ However, the character of the turbulence could be different in ITER and other reactors because burning plasmas are unlike the plasmas used to validate theoretical predictions of density peaking: electron heating will be dominant, high β will be required, and with high temperatures, the gyroBohm normalized fluxes will be small so the turbulence-driving gradients will be just above threshold. In present tokamaks, strong electron heating is in some cases experimentally associated with strong peaking and this is theoretically expected.⁶ However, electromagnetic effects may reduce adiabaticity and the dominance of the trapped electron pinch,⁵ although new simulations tailored for ITER suggest that this reduction in peaking will be smaller than it is for AUG.⁸⁴ Further nonlinear simulations with ITER-relevant heat fluxes and other parameters are needed to establish fully what we should expect in ITER.

We thank C. Angioni for enlightening discussions, D. R. Ernst for his least-squares spline-fitting software, and M. A. Chilenski for providing his Gaussian process regression fitting software.⁶⁷ This work was supported by U.S. Department of Energy Contract Nos. DE-AC02-09CH11466, DE-FC02-99ER54512, and DE-FG02-95ER54309, and used data obtained from the Alcator C-Mod tokamak, a DOE Office of Science user facility. This research used resources of the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. Use of parallel computer clusters at PPPL is also gratefully acknowledged. Any opinion, finding, conclusion, or recommendation expressed in this material are those of the author and do not necessarily reflect the views of the National Science Foundation.