Transport and confinement physics: Chapter 2 of the special issue: on the path to tokamak burning plasma operation

Two approaches have been employed to obtain further understanding of density peaking and its dependencies on various plasma parameters: regression over large databases and analysis of selected discharges close to ITER parameters. Special attention has been devoted firstly to confirm the role of collisionality in density peaking and secondly, to quantifying the role of the core NBI fueling source in contributing to density peaking in various conditions.

Several studies using the regression method over a large database using data from individual machines have been reported [4, 6–8]. These single machine density peaking scaling studies indicated that collisionality (often in these papers estimated as $\nu_{\mathrm{eff}}$) has statistically the strongest correlation with density peaking while the NBI particle source, temperature ratio, magnetic shear and Greenwald fraction $n_\mathrm G$ have somewhat smaller, but still significant correlation with the density peaking. The need to eliminate or reduce some of the parameters in the density peaking scaling led to creating multi-machine density peaking databases within the ITPA TC working group [12, 13]. In the JET-AUG combined density peaking scaling database [12], $\nu_{\mathrm{eff}}$ was found to be statistically the most relevant and important parameter in multi-variable regressions to describe density peaking. The contribution of the beam fueling alone was found insufficient to describe the observed variation of density peaking. The multi-machine scaling showed that the observed density peaking can only be explained by the existence of an anomalous inward pinch in low collisionality H-mode plasmas. This experimental multi-machine scaling on the role of collisionality in contributing to density peaking is consistent with the particle transport theory described in section 2.2. The derived scaling for density peaking, which includes $\nu_{\mathrm{eff}}$ in the regression variables, predicts a peaked density profile for ITER. The estimated peaking factor is of the order of $n_\mathrm e(\rho_{\mathrm{pol}} = 0.2)/\langle n_\mathrm e \rangle = 1.45$ where $\langle n_\mathrm e \rangle$ is the volume averaged density. In JET-AUG density peaking database, density peaking was also strongly correlated with $n_\mathrm e/n_\mathrm G$, the ratio of the plasma density to the Greenwald density limit. Based only on two machines, it was not possible to eliminate the latter parameter as one of the most important scaling variables. Without knowing the relative importance between $\nu_{\mathrm{eff}}$ and $n_\mathrm e/n_\mathrm G$, it is not possible to predict unambiguously density peaking for ITER. The addition of C-Mod data allowed the correlation to be broken, strengthening the case for $\nu_{\mathrm{eff}}$ as the controlling parameter. Data from the three devices strongly suggest that density peaking in ELMy H-mode plasmas in ITER will be in the range $n_\mathrm e(\rho_{\mathrm{pol}} = 0.2)/\langle n_\mathrm e \rangle = 1.4$–1.6, as illustrated in figure 3. The relative importance between $\nu_\mathrm{eff}$ and $n_\mathrm e/n_\mathrm G$ as the principal scaling parameter for density peaking is illustrated as the difference between figures 3(a) and (b).

Zoom In Zoom Out Reset image size

Further studies using similar approaches were performed as independent complements to the previous studies on density peaking in the JET tokamak. These focused on low collisionality data as the most representative of the reactor conditions [34]. The study confirms that $\nu^*$ is the most important parameter governing density peaking in H-mode, followed by the NBI particle flux and the $T_\mathrm i/T_\mathrm e$ temperature ratio. Moreover, similar studies on electron density peaking were performed in JT-60U tokamak and LHD helical plasmas at low collisionality [35]. The electron density peaking increased with decreasing collisionality in ELMy H-mode plasmas of JT-60U when the collisionality at half the minor radius was in the collisionless regime. The collisionality dependence of density profiles in LHD plasmas was similar to that in JT-60U plasmas in the same collisionality regime when neoclassical transport was reduced by geometrical optimization. JT-60U data with respect to density peaking versus $\nu^*$ fits well in a multi-machine database as shown in [9].

There are other approaches to study density peaking in particular with respect to the possible role of the NBI particle source. In JET, an analysis of a restricted number of JET plasmas with sufficiently ITER-similar parameters (low NBI, ITER like low $\nu^*$ and q₉₅) [20] was performed. It showed a clear correlation between density profile and particle flux. As a consequence, particle sources must be taken into account when extrapolating the density profiles to ITER or other next-step devices. Transport simulations of similar JET plasmas by adding the particle source due to NBI and an anomalous pinch proportional to $D \nabla q/q$ have been performed [36]. These simulations gave results that are compatible with the experimental measurements, however, the effects of the two terms are comparable and difficult to separate and compare in magnitude. The presence of an anomalous pinch can only be clearly identified in shots heated mainly by ion cyclotron resonance heating (ICRH) where the NBI particle source is negligible and, therefore, further experiment to quantify the magnitude of NBI fueling was needed.

Consequently, the role of NBI fueling was studied on JET with a dedicated experiment. This is needed as the analysis is complicated by the fact that often $T_\mathrm i/T_\mathrm e$, $\nu^*$ and the NBI particle source are correlated, and simultaneously, theyhave an opposite influence on the density peaking [34]. As a conseqeunce, this correlation between density peaking and the source may cause such a bias in the database type of analysis in such a way that the role of NBI cannot be resolved precisely. An H-mode experiment with a comparison of 8 MW of NBI heating/fueling with an identity discharge but with 8 MW of ICRH heating was performed in JET [18]. In the plasma with NBI fueling, density peaking was increased by a factor of two. This experimental result where the NBI fueling increases the density peaking by a factor of two was confirmed both with gyro-kinetic simulations and with full-radius transport modeling. On the other hand, on AUG in an experiment with plasma density and auxiliary power levels scanned to cover a wide parameter space indicated that the density peaking is mainly sustained by turbulence, with only a minor role for the NBI fueling source [19]. These at first sight contradicting results from AUG and JET with respect to the role of NBI fueling can be explained by the different turbulence properties and partly by the machine size. When the turbulence intensity and particle diffusion coefficients are large, NBI fueling has a minor role while in the opposite case it may even dominate. Both on JET and AUG, the NBI fueling is quite centrally located which will not be the case for example in ITER. However, there is not yet a full characterization on the plasma regimes in which the NBI fueling has a significant role in density peaking and in which it does not have that.

Furthermore, electron density peaking has been found to increase in electron cyclotron resonance heating (ECRH) heated plasmas on AUG with respect to the corresponding NBI ones [37]. The higher density peaking with increased central electron heating is explained by a concurrent reduction of the outward convection and a simultaneous increase in the inward thermo-diffusion. Further for AUG, density peaking was found to correlate with ITG-TEM turbulence, the highest peaking observed in the vicinity of ITG-TEM threshold [27, 38]. Similar conclusion was reached on DIII-D with NBI versus ECRH heating, i.e. the maximum density peaking takes place at ITG-TEM threshold, and then rapidly decreases when moving deeper into the TEM regime [39, 40]. Concerning the effect of rotation on density peaking on DIII-D, a strong increase in the inward particle pinch and consequent increase in density peaking in the counter-torque injected plasma is reported [40] In ITER, due to large volume with respect to NBI fueling, the role of NBI fueling in contributing to density peaking is not believed to be dominant and in most operational scenarios it will not be even particularly significant.

Major progress has been achieved in measuring local particle transport coefficients in the last decade, thanks to the improvements in profile reflectometry and High-Resolution Thomson Scattering (HRTS) diagnostics, such as those ones upgraded in JET [41, 42]. Local particle transport has been analyzed to quantify the magnitudes of particle diffusion, particle convection and thereby the NBI fueling in contributing to density peaking. The particle flux can be expressed

$$\begin{equation} \Gamma = - D \nabla n + nV,\end{equation} \tag{ 3 }$$

where D is the particle diffusivity and V is the convection velocity. They are defined separately for electrons $D_{\mathrm{e}}$ and ions $D_{\mathrm{i}}$. A typical notation is that the negative value of $V_{\mathrm{e}}$ or $V_{\mathrm{i}}$ indicates an inward particle pinch. The local particle transport experiments and analyses are described separately for electron particle transport and ion particle transport below.

2.3.2.1. Electron particle transport.

Plasma core electron particle transport has been studied by executing numerous dimensionless collisionality $\nu^*$ scans while simultaneously keeping otherdimensionless profiles fixed in many different plasma scenarios on JET and DIII-D [14, 15]. These experiments are the first ones where electron particle transport coefficients (diffusion and convection) in H-mode were measured on JET with high resolution diagnostics. A low-frequency gas puff modulation technique has been developed and optimized in JET to study core particle transport coefficients as well as core and edge particle sources [43]. The local electron density response to the gas puff modulation was measured with a multi-band density profile reflectometry which iscapable of both good spatial and temporal resolution, and complemented with HRTS diagnostics.

The diffusive and convective electron particle transport coefficients have been measured in JET H-mode plasmas where $\nu^*$ has been scanned over 0.1–0.5 by keeping other dimensionless plasma parameters $\rho^*, \beta_N, q$ and $T_\mathrm e/T_\mathrm i$ constant. Similar three-point dimensionless collisionality scan was performed in L-mode in JET. The electron particle transport coefficients from both the H-mode and L-mode dimensionless $\nu^*$ scans are presented in figure 4. This also demonstrates the existence of the inward particle pinch from the experimental side in both L- and H-mode. The gyro-kinetic code GENE and TGLF and EDWM transport models have also been used to calculate diffusion and pinch for the same $\nu^*$ scans [44–46]. The simulation results showed good match with the experimentally measured electron particle transport coefficients. The particle transport coefficients are significantly larger in L-mode plasma than in the H-mode plasma, giving rise to the fact that the NBI fueling plays a larger role in H-mode plasmas in contributing to density peaking, and a smaller role in L-mode plasmas. The significant role of NBI fueling in JET plasmas was already discussed in section 2.3.1. Another density modulation method, high frequency (up to 20 Hz) strike-point-sweeping, was also developed to study particle transport and edge fueling properties in JET [43]. This paper includes also the newly developed workflow to estimate the particle transport coefficients from a modulated density data beyond the previously used Takenaga’s formulas [47].

Zoom In Zoom Out Reset image size

A similar three-point $\nu^*$ scan in H-mode was also performed in DIII-D [15]. Similarly to JET, $\nu^*$ was varied by a factor 5, while keeping $\rho^*$, β_N, q and $T_\mathrm e/T_\mathrm i$ fixed. Moreover, a 3 Hz gas puff modulation was applied to modulate the density profile and extract the particle transport coefficients. The perturbative particle transport analysis indicated that the increase in density peaking at small $\nu^*$ is associated with an increase in the inward particle pinch and not an increase in core fueling. Furthermore, the experimental dimensionless scan confirmed that there is an increase in density peaking towards lower $\nu^*$, linked to a shift in the turbulence regime from ITG towards TEM. The existence of the electron particle pinch with high-quality diagnostics was also reported earlier on DIII-D [39, 40].

The electron particle diffusion coefficient and electron convection velocity in ECRH heated plasmas have been evaluated using a density modulation experiment in LHD [48]. The values of the perturbative electron particle diffusion coefficient, $D_{\mathrm{e,mod}}$, were significantly higher than those obtained by neoclassical prediction. On the other hand, the perturbative electron particle convection coefficients, $V_{\mathrm{e,mod}}$, were comparable with the neoclassical estimations. Turbulence measurement at ion-scale was compared with the measured $D_{\mathrm{e,mod}}$ and $V_{\mathrm{e,mod}}$, suggesting a clear link between the two. This low anomalous pinch of the same order as the neo-classical pinch seems like a specific feature in stellarator particle transport not observed in tokamaks.

2.3.2.2. Ion particle transport.

The theoretical prediction for the isotope scaling of particle transport was discussed in section 2.2, concluding that for otherwise similar parameters except the isotope, the particle pinch decreases with increasing mass of the hydrogen isotope. This leads to the theoretical result that the density gradient at zero particle flux is largest in hydrogen, and progressively decreases in deuterium and tritium plasmas, with a difference among the isotopes which increases with increasing collisionality [23].

A dimensionless ion mass scaling experiment between the hydrogen and deuterium plasma in JET was performed to study confinement as well as heat and particle transport [49]. Experimentally, no isotope scaling in particle transport was found between hydrogen and deuterium at matched $\rho^*_i$, $\nu^*$ β_N, q and $T_\mathrm e/T_\mathrm i$. This is confirmed in TGLF simulations with JINTRAC transport code where an agreement between the experimental and simulated density profiles was achieved between the dimensionlessly matched hydrogen and deuterium plasmas. The collisionality $\nu^*$ was small in these experiments and therefore, we do not expect to see a collisionality dependence in the particle transport in this regime, in accordance with theory [23] in section 2.2. At larger collisionalitity regime, theory predicts an isotope dependence on collisonality [23].

In the recent JET DT campaign, a similar, dimensionless ion isotope mass scaling experiment in L-mode between pure deuterium and pure tritium plasmas with matched $\rho^*$, $\nu^*$, β_n, q and $T_\mathrm e/T_\mathrm i$ was also performed in JET [50]. Moreover in recent JET DT campaign, a dimensional isotope scaling experiment between hydrogen, deuterium and tritium was performed [51]. In the dimensionless isotope mass scaling experiment, the isotope mass dependence in particle transport was found to be negligible, supported also by the perturbative particle transport analysis with gas puff modulation [50]. The same conclusion, i.e. negligible isotope scaling in particle transport and density peaking was found in H-mode [51]. Local GENE simulations in the L-mode pulses indicated, in agreement with the experimental results, a negligible effect on particle transport. On the other hand, in TGLF-SAT2 simulations, the effective turbulent particle diffusivity shows a very large increase in the tritium discharge outside of mid-radius, which is due to the change of the isotope mass. Similar observations, i.e. too flat density profiles were predicted for the H-mode discharges outside the mid-radius. This unfavorable isotope scaling of particle transport from D to T seen in the TGLF-SAT2 simulations on the outer radii at $0.5 \unicode{x2A7D} \rho \unicode{x2A7D} 0.8$ is not observed in the experiment. Further work is needed to validate the particle transport model in TGLF-SAT2 to obtain the experimentally observed isotope scaling, when moving from deuterium plasma to tritium plasma, and this work is already ongoing [52].

Isotope scaling of particle transport was also studied in LHD and CHS stellarators [53–55]. A reduction, favoring the deuterium plasma, in particle transport and turbulence characteristics between the hydrogen and deuterium plasmas was found in the core region at $\rho \unicode{x2A7D} 0.8$ in high collisionality regime on LHD. Gas puff modulation experiments in the high collisionality regime also showed correspondingly smaller particle diffusion coefficients in the deuterium plasma with respect to the hydrogen plasma. On the other hand, no such reduction in particle transport and turbulence was observed in low collisionality regime. On CHS, the particle confinement was improved in the deuterium dominant plasma with respect to the hydrogen plasma. The gas puff modulation experiments showed the quantitative difference in the particle transport coefficients, i.e. a lower particle diffusivity and a larger inward pinch velocity were observed in the deuterium plasma with respect to the hydrogen one.

2.4.1.1. Pellet ablation.

2.4.1.2. Plasmoid drift.

As described above the pellet deposition is a result of two simultaneous processes: pellet ablation and redeposition of pellet material. Such description defines the structure of any pellet deposition code. There are several leading codes that simulate pellet deposition process: PELDEPD2 [78], HPI2 [79] and PAM [80–82]. These codes differ in detail of calculations of such parameters as pellet lifetime, the size of plasmoids emanating from pellet and their trajectory in ambient plasma. In this respect recent PAM code can predict the pellet ablation and deposition very quickly while maintaining good accuracy. Finally pellet deposition codes also differ into which transport code the pellet model is implemented. It has to be noted that conventional transport codes describe the particle transport after pellet is ablated and homogenized on magnetic surfaces. The time of changeover from pellet deposition to conventional particle transport approximately coincides with the end of pellet evaporation processes and thus is conveniently handled by plug in of pellet deposition module into transport codes.

After a pellet is deposited, the plasma density quickly homogenizes along the flux surfaces. Such a post-pellet density profile takes the form of bump located at the plasma periphery. In ITER the maximum of the perturbation is located at $r_\mathrm{abl}/a \lt 0.82$ and its amplitude is about 30% [83] as seen in figure 5. Such a density perturbation creates two distinct zones: (1) inner zone with inverted density gradient and (2) outer zone with increased density gradient. These two zones broadly coincide with zones of inward and outward particle flux, respectively. The overall pellet fueling efficiency then depends on the balance between inward and outward particle fluxes in the pellet deposition zone. For example, when the inward flux dominates and the outward flux is small, then the fueling efficiency will be high. In the opposite case, when the inward flux is small and the pellet material is transported outwards, the fueling efficiency of the plasma will be low. The real situation is more complicated as inward and outward particle fluxes act simultaneously and also change in time during the pellet cycle.

Zoom In Zoom Out Reset image size

2.4.3.1. Inward flux (inner pellet deposition zone).

2.4.3.2. Inward propagation of isotope mix.

2.4.3.3. Outward flux (outer pellet deposition zone).

Deposition of pellets in the plasma periphery creates a zone with increased (negative) density gradient ($\partial n/\partial r \lt 0$). On ITER this zone spans from normalized minor radius $r/a \sim 0.8$ of up to the separatrix. Outward particle flux in this zone determines the pellet fueling throughput needed to keep the plasma density at the required value. In other words, it is the particle transport in this zone which controls the pellet fueling efficiency and ultimately the reactor burn-up ratio [88]. The later is the most critical parameter of fusion power plant as it determines the amount of tritium on site.

Particle transport in the outer pellet deposition zone is complex. Firstly, it has to be remembered that the outward flux is dynamic and changes during the pellet cycle. Secondly the outer pellet deposition zone (i.e. zone with negative density gradient) covers both the pedestal and the pedestal top regions so that the particle transport is not expected to be described by single physics model. For practical reasons, one can treat the outer flux as a sum of two fluxes: (1) flux due to micro-turbulence and (2) flux resulting from ELM mitigation or suppression. These are now discussed in next paragraphs.

The modes due to micro-turbulence can be studied by a conventional gyro-kinetic model as discussed above for inner pellet deposition zone. In a simulation for the AUG tokamak [87] it is found that the dominant mode is driven by an increased negative density gradient due to the pellets, and corresponds to the usual TEM instability. This mode produces large diffusive particle flux directed outwards, and becomes more unstable with increasing trapped particle fraction and with decreasing collisionality. This means that this mode is expected to be relevant at reactor conditions.

Concerning the flux due to ELM suppression/mitigation, there is now a consensus that, in fusion power plants, ELMs should be mitigated or completely suppressed in order to protect the divertor from high power bursts. Several ELM control methods are considered to achieve this, such as application of RMPs, pellet ELM pacing, and utilizing natural MHD modes. Among these methods only one regime has been subject to pellet fueling so far, namely the application of RMPs. The second method, testing pellet fueling on a plasma with pellet ELM pacing, would require two pellet injectors and such a system is not yet available in present tokamaks. Nevertheless a pellet injection system, including two pellet injectors, has been installed on EAST for exploring ELM pacing performance on EAST simultaneously with pellet fueling. Finally, pellet fueling of plasmas with ELM suppression/mitigation by natural MHD modes has not yet been tested.

Pellet fueling of plasma with RMP ELM mitigation was tested on MAST [97], EAST [98, 99], AUG [100] and DIII-D [101, 102]. A common feature of ELM mitigation/suppression is the so called density pump-out i.e. increased outward particle flux caused by RMP fields. These experiments showed that pump-out can be refueled by pellets. The pellet particle flux $\Phi_\mathrm{pel}$ required for refueling of density pump-out on AUG is [103]

$$\begin{equation} \Phi_\mathrm{pel}\sim0.07P_\mathrm{aux}/T_\mathrm{i,ped}\end{equation} \tag{ 4 }$$

where $P_\mathrm{aux}$ is the heating power and $T_\mathrm{i,ped}$ is the pedestal ion temperature. When the density pump-out is refueled by pellets, the plasma transitions from ELM suppression to ELMy regime, but a significant degree of ELM mitigation is preserved. The pellet particle flux required will, therefore, depend on which level of mitigation/suppression plasma has to be controlled.

Phenomenology of the transition from ELM suppression to an ELMy regime was studied in detail in AUG [100] using long trains of shallow pellets. It is found that each pellet triggers an ELM-like event which is followed by a transient (20 ms) ELM free H-mode, and then the plasma returns to the pre-pellet phase with ELM suppression. When the plasma density gradually increases (spontaneously or by increasing pellet rate) then ELMs start to occur during the post-pellet H-mode phase with increasing frequency until ELM suppression is completely lost in the pellet cycle. This behavior is in line with existence of density limit for ELM suppression and is also observed with gas fueling.

The physics basis for particle transport in the outer pellet deposition zone under condition of ELM suppression/mitigation is under development. Based on experimental observations the main elements of such model should be (1) the particle transport in presence of RMPs and (2) role of the steep density gradient due to pellets [104]. To identify the character of overall particle transport it would be useful to distinguish between convective or conductive losses. In MAST [97] post-pellet ELM losses seem to be convective as the pellet material is removed by ELM outwards even from inner pellet deposition zone.

Pioneering results on helium density profiles and transport were obtained in TFTR [105–107] and DIII-D [108]. First measurements of helium transport coefficients [105] revealed that diffusion was at least one order of magnitude higher than neoclassical at the edge and decreased towards the plasma center. Similar radial variation was observed for the inward directed convection, with centrally peaked helium density profiles. The inward convection was measured to be stronger than the neoclassical predictions in TFTR supershots and the helium diffusivity was observed to increase with increasing ion heat transport [106]. The confinement time of helium ash produced by D–T fusion reactions in TFTR was observed to be dominated by the edge pumping rates and up to one order of magnitude longer than the energy confinement time [107], still compatible with ignition in a reactor [116]. At DIII–D, an extensive analysis of helium density profile shapes in deuterium plasmas in various regimes showed that these have similar shapes to the profiles of the main plasma electron density [108]. Moreover, helium diffusivities were also measured which were of the same size as the effective heat conductivity $D_{\mathrm{He}} / \chi_{\mathrm{eff}} \sim 1$ [108] and, therefore, higher than the neoclassical expectations. Research in JT-60U [122, 123] also found helium diffusivities which were significantly larger than the neoclassical predictions, with the additional result that no observation of central accumulation of light impurities like helium and carbon was obtained even in internal transport barriers which featured ion heat transport at the neoclassical level. In contrast, an impurity with higher charge, such as argon, was observed to accumulate. The unfavorable behavior of impurities with increasing charge in the presence of internal transport barriers was also observed at JET [124], consistent with the predictions of neoclassical transport theory.

A reduction of the argon accumulation, corresponding to a flattening of the argon density profile, was found to occur as a consequence of the application of central ECRH [123]. The leading mechanism producing the flattening was identified to be connected with a reduction of the peaking of the main plasma density profile, also produced by the central electron heating in the high $\beta_\mathrm{p}$ transport barriers of JT-60U. At the same time, experiments in ASDEX Upgrade in normal H–mode plasmas also demonstrated the beneficial role of central electron heating in reducing the peaking of the impurity density profiles [125]. Silicon laser ablation experiments revealed that central ECRH significantly increased the Si central diffusivity to values which were clearly above the neoclassical level, while, concomitantly, reducing the inward convection or even reversing it to an outward convection [125]. This body of results provided the basis for one of the main control tools to avoid central accumulation of heavy impurities in the operation of present tokamaks with tungsten plasma facing components [126], namely, the application of central wave heating. As it will be described below in more detail, in these conditions the electron heating fraction plays the dominant role.

At the same time, these experimental evidences clearly showed that predictions of impurity transport could not be based on collisional transport only, and this motivated the investigation and development of a more comprehensive theoretical description of turbulent impurity transport [29, 127–131]. The main mechanisms of turbulent impurity transport have direct analogues to those which produce the transport of the main plasma species. In addition to the diagonal diffusive component, connected to the density gradient of the impurity itself, off–diagonal contributions are also present. Pure convective terms are produced by the compression of the E×B drift in toroidal geometry and by the impurity parallel dynamics [128]. Moreover, the presence of the temperature gradient of the impurity causes an additional off-diagonal contribution, usually dubbed thermo-diffusion [128–130]. A more complete description also includes an effect which was identified only later [132] and which provides an additional convective component, proportional to the gradient of the toroidal rotation, dubbed roto–diffusion. A comprehensive review of these mechanisms can be found in section 3.3 of [113] and is not reported here. It is, however, important for what follows to recall some of the general properties of turbulent impurity transport, which are also connected with several aspects of the comparisons with the experimental results. The ratio of the turbulent impurity diffusion coefficient to the turbulent heat conductivities is theoretically found not to be constant, but in fact to feature a strong non–monotonic dependence on the electron to ion heat flux ratio. This is produced by the different types of turbulence, from ITG to TEM, and the corresponding impact of the relevant resonance, determined by the mode frequency, on the different transport channels [31]. This effect is particularly strong for highly charged impurities, which are predicted to reach values of diffusivities that are significantly larger than the effective turbulent heat conductivity, in conditions when electron and ion heat fluxes become comparable. In contrast, for large dominant ion or electron heat fluxes, the ratio of the impurity diffusivity to the heat conductivity becomes significantly smaller than one, resulting in a variation of this ratio by more than one order of magnitude depending on the value of the electron to ion heat flux ratio. Off–diagonal components also exhibit the property that their direction can change sign depending on the type of turbulence [27]. In particular both thermo–diffusion and roto–diffusion are directed outward in the general conditions of dominant ITG turbulence, and oppose the usually inward directed pure convective term. In contrast, they can be directed inward in strong TEM turbulence. Moreover, in conditions of comparable electron and ion heat fluxes, also subdominant modes can play an important role in determining the impurity flux, as confirmed by comparisons with nonlinear results [32].

Concomitant with the first theoretical works on the properties of turbulent impurity transport, a growing number of experimental observations [133, 134] allowed the first dedicated comparisons between theoretical predictions and experimental results. These first applied still limited theoretical models [135–137], but already revealed some elements of general, qualitative agreement, particularly in the overall consistency that no evidence of impurity accumulation in conditions of dominant turbulent transport was obtained, in agreement with general expectations from theory. However, also elements of disagreement were identified, which motivated further research.

An increase of the central diffusivity of nickel with increasing electron temperature logarithmic gradient was reported in Tore Supra [138, 139]. This increase was demonstrated to be of turbulent origin, and to be consistent with theoretical quasi–linear calculations of impurity diffusivity. Two complementary works for the validation of turbulent impurity transport theory were performed in ASDEX Upgrade and Alcator C–Mod. The first [37] was dedicated to the comparison of quasi–linear gyrokinetic predictions with a set of intrinsic boron density profiles measured in H–mode plasmas [140] with NBI heating and different levels of central ECRH. Here an experimental trend, which is largely confirmed by future works, is identified. The boron density profiles are increasingly peaked with increasing central ECRH, while still exhibiting a peaking that is lower than that of the electron density profile. The electron density profile also peaks in the confinement region with increasing ECRH power in these conditions, consistent with theoretical predictions of electron particle transport in ITG turbulence, see discussion in section 2.2. In contrast, with NBI heating only, the boron density profile becomes flat or even hollow around mid–radius. These experimental trends are at least qualitatively consistent with the combined outward roles of thermo–diffusion and roto–diffusion in ITG turbulence, which increase with decreasing ECRH heating. However, as described later, the pursuit of quantitative agreement with the increasingly accurate measurements of the observed hollowness of the impurity density profiles has been further investigated in the ITPA activity. In fact, this aspect still challenges the theoretical predictions, despite several attempts and recent additional improvements in the completeness of the theoretical description. The second work [141] compares, for the first time, global non–linear gyrokinetic simulations with the gyrokinetic code GYRO [142], matching the heat fluxes, with the separate measurements of the diffusion and convection coefficients of a calcium impurity, injected in Alcator C–Mod L–mode ICRH plasmas with dedicated laser ablation experiments [143]. The results provide a demonstration, for the first time, that, within the experimental uncertainties, gyrokinetic simulations can reproduce experimentally estimated values of the impurity transport coefficients, with large diffusion and large inward convection.

These promising results motivated more specific and extended research, aiming at improving the quality of the measurements, extending the investigated domain of parameters and thereby more strongly constraining the theoretical validation and informing on the reliability of our current prediction capabilities of the impurity density profiles predominantly determined by turbulent transport. The impact of experimental uncertainties were comprehensively considered [144], and linear and non–linear gyrokinetic simulations were found to reproduce the trend of hollow boron density profiles with strong NBI heating and strong plasma rotation, but an under-prediction of the level of hollowness was regularly reported. Impurity transport modeling activity at JET consistently documented the inability of the gyrokinetic modeling to reproduce the observed hollowness of the intrinsic carbon density profiles [145, 146]. In these studies also the dependence of the impurity peaking factor with increasing impurity charge was investigated, and found to be qualitatively reproduced by the theoretical predictions, with an increase of the peaking factor with increasing impurity charge for Z $\unicode{x2A7D}$ 20 and a saturation to constant values for larger values of Z, much below the neoclassical predictions. This is consistent with the theoretical expectation of turbulent impurity transport which becomes independent of the impurity charge at high charge.

These considerations in connection with the disagreement on the prediction of hollow impurity density profiles has further motivated the joint research under the ITPA. Analyses and modeling of both JET C–Wall [147] intrinsic carbon profiles and JET–ILW [148] profiles of extrinsic helium, berillium, carbon, nitrogen and neon confirmed that the gyrokinetic simulations were predicting weakly peaked impurity density profiles in conditions where measured profiles were flat or hollow. Moreover, while the experimental results exhibit a significant variation in impurity density peaking, with helium as peaked as the electron density profile, and flat to slightly hollow carbon and nitrogen around mid–radius, the theoretical results are characterized by a more limited variation of the density peaking.

This has also motivated new experimental studies on intrinsic impurities at ASDEX Upgrade [149], with an extensive comparison of boron and helium density profiles, as well as on injected impurities with transient transport experiments, separately determining diffusion and convection in DIII–D [150] and ASDEX Upgrade [151, 152]. Consistent with the specific profile measurements of different impurities in JET–ILW [148], in [149] a large database of helium density profiles is observed to have a density peaking which mainly follows the peaking of the electron density, whereas the corresponding profiles of boron are generally flatter than the electron density. Cases with hollow density profiles were observed for both helium and boron. The comparison with quasi–linear gyrokinetic calculations of turbulent impurity transport at mid–radius is shown in figure 6. The main parametric trends are correctly reproduced by the gyrokinetic results. However, the helium normalized logarithmic density gradient $R/L_{n\mathrm{He}}$ is systematically underpredicted in the low NBI power cases, whereas, in contrast, the boron normalized logarithmic density gradient $R/L_{n\mathrm B}$ is usually underpredicted for the low NBI power cases, and is largely overpredicted for the high NBI power cases, in particular at low ECRH power, where the observed boron hollowness is not reproduced.

Zoom In Zoom Out Reset image size

These extensive datasets of boron and helium profiles were also made possible by upgrades and improvements in the CXRS diagnostics, for light impurities [153] and helium in particular [111]. In [153], the importance of the additional emission from the beam halo population for the interpretation of the CXRS spectra has been highlighted. In [111], a kinetic forward model for the helium plume [109] emission has been developed.

In [150, 151], where the experimental analysis of the injected impurities transients allows the separate determination of diffusion and convection, measurements are compared with both quasi–linear gyrofluid and quasi–linear and nonlinear gyrokinetic simulations respectively. In [150], dedicated to high confinement regimes at DIII–D without edge localized modes (ELMs), agreement with quasi–linear gyrofluid predictions is mainly identified in the convective term, not in the diffusion. In [151], a 0.5 MW ECRH and 5 MW NBI heated H-mode plasma at ASDEX Upgrade, with additional time modulation of 1 MW ICRH power to produce a modulated edge boron source from the antenna limiters, reports experimentally determined boron diffusion and inward convection coefficients both quantitatively consistent with quasi–linear and heat flux matched nonlinear gyrokinetic predictions within experimental uncertainties. It is worth observing that this case falls in the domain of parameters where an agreement between gyrokinetic results and observations can be expected from the results in [149].

The systematic analyses of databases of intrinsic carbon and boron density profiles, respectively from JET C-Wall [147] and ASDEX Upgrade [149] H-mode plasmas, allowed the identification of the main parameter domains in which the disagreement between theory and carbon and boron experimental profiles was large. In JET the disagreement is found to be particularly large at high collisionality, while in ASDEX Upgrade the disagreement is found to be particularly large at high ion temperature logarithmic gradient and rotation gradients [149], correlating with large NBI powers. The comparison between the data of the two devices, performed within the ITPA joint activity, showed that the ASDEX Upgrade cases with higher ion temperature and rotation gradients were close to the JET parameter space investigated in [147] and confirmed the overprediction of impurity density peaking found in that previous work. These results motivated two lines of theoretical investigation, which were performed within the ITPA group. The first, related to the high collisionality, lead to the study of the direct impact of a neoclassical equilibrium on the turbulent impurity transport [154]. The second, connected with the correlation of the disagreement with high NBI power, focused on the impact of fast ions [155]. The first study found only minor effects, which become negligible for light impurities like boron or carbon at core plasma conditions. The second study identified an important additional ingredient, the beam ion density, which was missing in the modeling activity performed so far, and which leads to a strong reduction or even reversal of the convection at high values of beam ion density and gradients, due to a combination of effects involving both turbulent and neoclassical transport.

In Alcator C–Mod separate diffusive and convective transport coefficients of calcium have been experimentally inferred with a fully Bayesian approach and the corresponding measured density profiles [156]. The results have been compared with theoretical predictions obtained from both quasi–linear gyro–fluid and non–linear gyrokinetic simulations [157], and a good agreement has been found in general, particularly for the diffusion coefficient of all examined cases and for conditions with inward impurity convection, as shown in figure 7. At the same time, clear experimental evidence of outward impurity convection and hollow calcium density profile is obtained in an ICRH EDA H-mode plasma (figure 2, bottom–right), which cannot be reproduced by the theoretical predictions. As compared to the other examined cases, in L–mode and I–mode, this EDA H–mode case features lower current, higher collisionality, lower logarithmic density gradient, higher logarithmic electron temperature gradient, and comparable logarithmic ITG. The results provide a strong indication that, while the impact of beam ions, as identified in [155], is certainly an important physical ingredient to be included in the theoretical modeling, at the same time it cannot be considered the only missing piece leading to the solution of the disagreement between observations and predictions in conditions of outward impurity convection.

Zoom In Zoom Out Reset image size

From the impurity diagnostic standpoint, clear evidence has mounted that impurity transport studies approaching the outer part of the plasma, also well inside of the pedestal, cannot ignore the impact of charge exchange between impurities and background (deuterium) neutrals [157–159]. This has important ramifications for spectroscopic measurements. As a result of these reactions, the ionization equilibrium of impurities can be radically modified, with higher charge states recombining to lower ones. This has been shown to be a fundamental ingredient for experimental analysis and interpretation of spectroscopic studies on ASDEX Upgrade [158], DIII–D [159], and even C-Mod [157], where plasma densities were much higher and neutral penetration comparably lower. Charge exchange with NBI neutrals, often used for diagnostic purposes, has also been quantified and observed to play a more limited role in modifying the ionization balance, while possibly still leading to misinterpretation of spectroscopic data in the core, especially for partially stripped ions. The consideration of the role of edge neutrals is challenging since their densities are poorly known. In the C-Mod [157] and DIII-D [159] works, neutral densities where estimated based on SOLPS-ITER runs using the EIRENE Monte Carlo neutral code, demonstrating a new element of integration between core and edge models. As the community expands its activities towards validation of pedestal and scrape–off layer (SOL) impurity transport, these efforts gain clear importance.

At ASDEX Upgrade, further developing the experimental technique which allows the separate identification of the diffusive and convective components of boron transport [151], a large database of boron transport coefficients has been assembled, allowing a practically unprecedented analysis of the experimental properties of light impurity transport in terms of separate diffusion and convection coefficients as a function of a variety of plasma parameters [160], in particular including a very large variation in NBI and ECRH powers, and a limited variation in edge safety factor. The normalized ion temperature logarithmic gradient $R/L_{Ti}$ is identified as the strongest ordering parameter for both boron diffusion and convection. It is also shown that a large value of $R/L_{Ti}$ is a necessary condition to obtain outward convection and hollow density profiles of boron. A critically important result, which confirms and extends previous knowledge [37, 125, 140], is that, while even large variations of the NBI power have relatively limited impact on the impurity transport, the inclusion of low levels of ECRH can dramatically increase both the diffusive and the convective coefficients, also leading to a moderate increase of the boron density peaking. Similar to the results on Alcator C–Mod, the comparison of these experimental results to a combination of neoclassical and quasi–linear gyrokinetic results shows a relatively good agreement on the diffusion coefficient and on the convection coefficient when this is inward directed. In particular, a non–monotonic dependence of the impurity diffusivity to heat conductivity ratio on the electron to ion heat flux ratio is demonstrated by the experimental analysis, qualitatively consistent with the theoretical predictions [31]. In contrast, even with the inclusion of the impact of beam ions, the outward convection measured in NBI dominated plasmas, featuring strong rotation and ITGs, is not reproduced by the theoretical predictions. A strong variation of the impurity diffusivity to heat conductivity ratio as a function of the electron to ion heat flux has been also observed in DIII-D experiments of aluminium transport [161], demonstrating a commonality in this generic property of turbulent impurity transport. A comparison of the experimental results, as well as of the related gyrokinetic calculations, obtained in these recent works dedicated to the separate experimental determination of the impurity diffusion and convection coefficients is presented in figure 8. The ratios of the impurity diffusion coefficients to heat conductivities for aluminium and boron in DIII-D and ASDEX Upgrade, respectively, are plotted as a function of the electron to ion heat flux ratio. In the left plot, from [161], the DIII-D experimental results are compared with nonlinear electrostatic gyrokinetic simulations. In these simulations the ion heat flux was matched by modifying the ITG in input. The experimental and theoretical aluminium diffusivity to heat conductivity ratios are plotted as a function of the electron to ion heat flux ratio from power balance and from simulations respectively. In the right plot, from [160], the ASDEX Upgrade experimental results are compared to quasi–linear gyrokinetic results. These have been performed with nominal input parameters from profile measurements only, without attempting any matching of the electron to ion heat flux ratio. In this case, both the experimental and the theoretical results of the boron diffusivity to heat conductivity ratio are plotted as a function of the experimental electron to ion heat flux ratio from power balance. Current work, still within the ITPA, is now dedicated to nonlinear electromagnetic gyrokinetic simulations of a subset of these data.

Zoom In Zoom Out Reset image size

An interesting result has recently been obtained within a statistical test–particle approach by the explicit inclusion of the polarization drift in the gyro-center test–particle motion [162]. Within this approach, the direction of the impact of the polarization drift on the impurity transport is determined by the mode propagation frequency and it is directed outward for ITG modes and inward for TEM modes. Therefore, it provides an additional contribution that, in particular, allows the ASDEX Upgrade experimental observation of hollow impurity density profiles with large NBI heating to be reproduced. In the conventional gyrokinetic formulation, the polarization effect is treated as a shielding term in the gyrokinetic Poisson equation, provided by the so–called polarization density [163]. Alternative formulations of the gyrokinetic system of equations which explicitly include the polarization drift have been derived [164, 165]. However, so far the effect of the polarization drift on the impurity transport has not been computed within a consistent gyrokinetic formalism. Moreover, an additional term in the impurity particle flux, which is produced by an externally applied toroidal momentum torque, typically by NBI, and has been largely neglected in the last decades, has been rederived very recently [166], confirming past results [167–169]. This term is directed outward with co-current injected torque and, therefore, also contributes to the explanation of the observation of hollow light impurity density profiles with high co-current NBI. It is caused by the interaction between parallel impurity inertia, poloidal flow damping and friction, and by the link between parallel motion and perpendicular particle flux caused by the poloidal stress tensor. It arises from the impact of the turbulent toroidal momentum flow on the radial neoclassical particle flux and, as such, is not included in separate descriptions of turbulent and collisional impurity transport, and its effect has been ignored so far in recent comparisons between theoretical predictions and experimental impurity density profiles. Its effect is estimated to be large enough to provide an observable modification of the predicted profile shape of light impurities, potentially leading to the observed level of hollowness in conditions of large NBI heating power [160]. In contrast it should be expected to become relatively small in comparison to the neoclassical convection of highly charged and heavy impurities, particularly in the concomitant presence of strong toroidal rotation, typical of high NBI heating powers.

Investigations of impurity transport were also performed in spherical tokamaks. The dominance of neoclassical transport for externally injected neon was first reported in NSTX H-mode plasmas with NBI heating [170]. The transport of intrinsic carbon and lithium impurities was analyzed at NSTX [171], consistently revealing a dominance of neoclassical transport effects. In particular, the analysis shows that carbon tends to accumulate during the ELM–free phases (induced by lithium coating). The observed accumulation is explained by neoclassical transport in connection with modifications of the ion temperature and density profiles, whereas the high carbon concentration enhances the neoclassical lithium diffusivities, limiting the lithium concentration to very low values. One of the few recent studies providing experimentally inferred helium transport coefficients has been performed on MAST [110], by helium gas puff transient transport analysis. In the outer part of the L–mode plasmas, helium transport is dominated by turbulence, providing large diffusion and inward convection. In contrast, in H–mode the diffusivity is found at the predicted neoclassical levels, without evidence of a pinch. The turbulent pinch in L–mode is consistent with the gyrokinetic predictions for collisionless TEM, in connection with steep density gradients, while ITG are stabilized by the strong E × B shear. Building on these helium transport results, also the charge dependencies of neoclassical and turbulent light impurity transport have been experimentally investigated in MAST [172]. While carbon and nitrogen are observed to be weakly screened in the core, with an outward convection at mid–radius, helium is peaked over the entire plasma radius. Over a large central portion of the plasma, the transport is at the neoclassical level, whereas it is turbulent in the outer region, outside mid–radius, where gyrokinetic predictions confirmed the dominant role of TEM. Here the modeling, combining neoclassical with quasi–linear turbulent transport, provides reasonable estimates of the impurity transport coefficients and the impurity charge dependence. The reversal from inward to outward convection observed around mid–radius when moving from the periphery to the center of the plasma for both carbon and nitrogen is found to be consistent with the neoclassical predictions of banana–plateau transport. Neoclassical transport can play a dominant role in this extended radial region for light impurities in a spherical tokamak, characterized by strong rotation, stabilizing ITG turbulence, and with TEM only unstable outside mid–radius, due to flat density profiles in the core. These conditions with dominance of neoclassical transport over a large portion of the radial domain appear to be specific to spherical tokamaks like MAST and NSTX, whereas in conventional aspect ratio tokamaks turbulent transport dominates over neoclassical transport for light impurities over a broader radial region approaching the plasma center. This implies that an explanation of the observed hollow impurity density profiles in the confinement region of conventional aspect ratio tokamaks has to be obtained from the turbulent transport predictions.

One theoretical element, that is not included in that very recent review (as it became clear even more recently), is related to the impact of plasma rotation on neoclassical impurity transport in the low collisionality limit. This aspect has been identified in very recent calculations of neoclassical tungsten transport with the drift–kinetic code NEO [174–176] applied to plasma parameters achieved in the KSTAR tokamak [177], which are characterized by very high rotation at low collisionality. It is shown that in this limit of the parameters, the relative effect of the screening produced by the ITG [119, 178] is increased, leading to more favorable conditions for the avoidance of central impurity accumulation. This behavior, which is characteristic of the impurity transport in the banana regime [179], was in fact already present in the results of the first study on the impact of rotation on impurity transport [180], although this specific favorable property was not highlighted in that pioneering paper. Recent experimental results at JET–ILW [181, 182] have entered this regime, where the collisionality is low enough to reach conditions in which the impact of the plasma toroidal rotation increases the relative strength of the temperature screening.

This behavior is also consistent with the modeling results of the transient phase of an ASDEX Upgrade plasma with very high energy confinement, low collisionality and very high rotation [183]. This regime can be expected to be achieved at the low collisionalities of a burning plasma, in case of sufficiently large rotation, and could be particularly interesting for low aspect ratio devices, which are often characterized by strong plasma rotation. In contrast, at high collisionality, with highly charged impurities which typically are in the Pfirsch–Schlüter regime, toroidal rotation not only strongly increases the neoclassical transport [184–188], but can also reduce the relative strength of the temperature screening [186, 189, 190] more easily leading to central accumulation of the impurities. This effect is clearly exacerbated in the operation with NBI heating, as already documented in first observations [191, 192], since NBI produces the highly unfavorable combination of a central particle source and an external torque, which both have a negative impact on the behavior of heavy, highly charged impurities like tungsten. This unfavorable effect has been clearly demonstrated on JET–ILW and ASDEX Upgrade [183, 193–195], and has been very recently confirmed also on EAST [196]. The role of plasma rotation in inducing poloidal asymmetries of the tungsten density distribution and in modifying the tungsten impurity flux has been experimentally verified in ASDEX Upgrade, consistent with the theoretical predictions [197]. The mitigation of this effect, which unfavorably impacts the operation of present tokamak plasmas in the most common high collisionality regime of heavy impurities, critically relies on the modification of the plasma profiles in order to reduce or reverse the inward convection, typically by reducing the ion density gradient and by increasing the ITG, and/or on the relative increase of the turbulent transport, which has to offset the neoclassical convection, by central wave heating and/or by reducing the plasma rotation and/or the poloidal asymmetry. These elements clearly explain the interest in understanding and realistically predicting the physics and the magnitude of the turbulent diffusion, already presented in section 3.2, given the critical role that it can have to avoid heavy impurity accumulation. At the same time, the strong impact of the inhomogeneity of the poloidal distribution of heavy impurities on the neoclassical perpendicular transport has renewed the interest in the study of the parallel transport of heavy impurities.

Here below, we briefly review the main aspects of these physical mechanisms, particularly those that have been considered within the ITPA.

There are three main effects which can break the homogeneous distribution of the impurity density on the magnetic flux surface. The first is the plasma toroidal rotation by the centrifugal force [180, 184–188, 198], the second is the impurity–ion friction [185, 188, 199–202], the third is the presence of particle species featuring significant temperature anisotropies, like resonant particle species in the presence of ICRH, which also develop non homogeneous poloidal distributions [203–209] and impact the impurity distribution by quasi–neutrality. Also the usually limited temperature anisotropy developed by NBI ions can have a non–negligible impact on the poloidal distribution of heavy impurities [197, 206]. Depending on the poloidal asymmetry which is generated, the neoclassical transport of the impurities, and in particular the Pfirsch-Schlüter transport, can be strongly increased, for instance in the case of strong low field side localization as produced by the centrifugal force, or decreased, for low HFS asymmetry which can offset the effect of the poloidal asymmetry of the magnetic field [186].

While the impact of toroidal rotation on the neoclassical transport has been known for a long time [180, 184, 185], the realization of its strong impact and the requirement of its inclusion in the modeling of experiments has been highlighted only recently, particularly within the ITPA [183, 189, 193], thanks to the modeling possibilities offered by the NEO drift–kinetic code. Figure 4 shows the impact of rotation on the NEO predictions for the diffusion and the convection terms of W for three time slices of a JET–ILW discharge, from [189]. Not only the diffusion and the convection are strongly increased, figure 9(a), but also the convection to diffusion ratio is modified, figure 9(b), and the convection becomes more inward, due to a reduction of the strength of the temperature screening component produced by the centrifugal effects. Considering the edge region of the plasma, including the H-mode pedestal, parallel friction [185, 199, 202] and centrifugal effects from poloidal rotation [210] can become important. These effects have been also shown to produce poloidal asymmetries in B2SOLPS5.2 simulations [211], with increased impurity density on the HFS in the region of the ETB, consistent with observations in AUG and MAST.

Zoom In Zoom Out Reset image size

The unfavorable properties of collisional transport imply that operation with highly charged metallic elements, such as tungsten or molybdenum, in the plasma facing components requires some precautions in order to avoid accumulation [212–214]. Two concomitant main lines have to be pursued for safe operation. The first element is related to the impurity concentration at the edge, which has to be kept small. This can be obtained by limiting the impurity sources from the wall, as well as by sufficient transport in the pedestal region. The latter can be produced by regular ELMs [215, 216] or by favorable pedestal conditions, in particular with strong ITGs [217]. ELMs are observed and predicted to produce a net tungsten outward transport in present conditions with inward neoclassical convection and peaked impurity density profiles in the pedestal. At the same time, in the case of weak main ion density gradients in the pedestal, and consequent outward neoclassical convection and hollow tungsten density profiles, ELMs are predicted to produce an opposite effect and increase the tungsten density in the core [216]. The impurity transport in the H–mode pedestal region is experimentally found to be mainly caused by neoclassical convection [215], which strongly increases the impurity confinement time in H–mode plasmas when the main ion density gradient inward component dominates the neoclassical convection.

In these conditions, a sufficiently high ELM frequency is required to limit the edge tungsten concentration, which typically imposes operational constraints to the minimal level of allowed gas puff. High gas puff levels can have detrimental effects on the energy confinement [218–223]. Non–stationary ELM–free conditions are prone to an increase of the tungsten concentration at the edge and eventually to a radiation collapse. The same property of relatively long impurity confinement times characterizes stationary high confinement regimes with small or no ELMs that maintain a H-mode density pedestal, featuring an ETB with reduced particle turbulent transport, such as the EDA H–mode regime [224, 225], as well as the quiescent H–mode regime and regimes characterized by active ELM-suppression [150]. However, while particle transport, and therefore also impurity transport, is in general reduced in the radial region of the density pedestal, specific properties of the different conditions which characterize the various high confinement regimes of stationary operation, with small or no ELMs, allow certain regimes to be less prone to, or even free from, problems connected to W accumulation. This is particularly the case of the conditions of the quasi-continuous exhaust regime, with small type II ELMs, [226], due to the high density in the SOL and consequent cold far SOL and reduced tungsten sources, in combination with possible increased transport close to the separatrix. This is also the case of conditions of ELM suppression with externally applied magnetic perturbations [227, 228], because of the usually lower edge density and reduced density pedestal and the possible increased transport in connection with the peripherally modified magnetic topology. Stationary ELM–free regimes, like the EDA H–mode [224, 225, 229, 230] are often characterized by the presence of quasi–coherent fluctuations at the edge. These keep the operational conditions below those which would lead to an ELM trigger and can also have favorable consequences on the impurity behavior, leading to a reduction of the impurity content in the plasma [224, 229, 230] with respect to similar plasma conditions featuring ELMs. A stronger reduction of the impurity confinement time with increasing heating power as compared to the reduction of the energy confinement time has been observed in the EDA H–mode regime in ASDEX Upgrade. The general property of relatively long impurity confinement times which characterizes the regimes with a strong density pedestal consistently also implies that, in contrast, a reduced impurity confinement time is present in regimes which are not characterized by a strong density gradient at the edge. This is the case of the L–mode or the I–mode, as demonstrated by a systematic comparison in Alcator C–Mod [225]. Here the increased particle transport as well as the reduced neoclassical inward convection, due to a reduced main ion density gradient at the edge, lead to significantly shorter impurity confinement times, where in particular the I–mode regime breaks the otherwise usual correlation between impurity confinement and energy confinement times.

A second element which allows the avoidance of heavy impurity accumulation is related to the impurity profile shape in the core. The impurity density profile has to feature a limited peaking in the core, particularly without developing strong central gradients close to the magnetic axis. A central particle source, as produced by operation with NBI heating, has an unfavorable impact by increasing the peaking of the main ion density profile close to the magnetic axis. In contrast, operation with wave heating only is clearly much less prone to accumulation, because of the combined beneficial effects of absence of a central particle source and of an external torque, as shown in both C–Mod [231] and more recently at WEST [232].

Central wave heating also has very beneficial effects in the case of operation with NBI heating, due to a combination of effects, which can favorably act concomitantly. A first element is provided by the increase of the impurity turbulent diffusion [125, 160, 161, 233]. The increase of central turbulent transport can also reduce the peaking of the density profile close to the magnetic axis, more strongly balancing the impact of the NBI particle source. An additional and related element, also provided by the increase of the central turbulent transport due to the central wave heating, is connected with the increase of central toroidal viscosity and the reduction of the plasma toroidal rotation [234], which can also be produced by a modification of the turbulent residual stress. These effects, mainly connected with an increase of central turbulent transport, can be considered to be common to both ECRH and ICRH. There are, however, also some differences between the effects of ECRH and ICRH. First, it has to be mentioned that ICRH with W–coated antennas has been observed to increase the tungsten sputtering, due to radio-frequency (RF) sheaths driven by parallel RF electric fields, an effect which has been strongly reduced on ASDEX Upgrade by offsetting the effect of the electric fields with a three–strap antenna [235]. Second, ICRH produces highly energetic ions with strong temperature anisotropy, providing an additional source of tungsten poloidal asymmetry [203–209].

Fast ions can also contribute to the total temperature screening [183], an effect which strongly depends on the actual peaking of the effective temperature profile of the fast ions. The predicted strength of this effect is observed to depend on the impact of finite orbit width effects in the modeling of the fast ion distribution [195]. Central ICRH and ECRH have the tendency to increase the overall turbulent transport level. This can reduce the impact of the central particle source, particularly in case NBI heating is applied, or of the Ware pinch, in conditions in which this is significant, close to the magnetic axis. Therefore, central ECRH and ICRH can reduce the central peaking of the main plasma density [191, 236]. Finally, ICRH produces both electron and ion heating, whereas ECRH only heats the electrons. While ion heating can be efficient in increasing the ITG and increasing the neoclassical screening, depending on the plasma heating conditions, the electron heating can have a stronger effect in increasing the turbulent diffusion [31]. A comparison on the efficiency of ICRH and ECRH in preventing tungsten accumulation in ASDEX Upgrade has shown that it is the electron heating fraction which orders the peaking of the tungsten concentration, providing experimental indication that, at least in the conditions of those experiments, it is the electron heating fraction that has to be maximized for a more efficient avoidance of tungsten accumulation [194]. It has to be underlined that while on–axis ECRH and ICRH are widely documented to have these beneficial effects in avoiding central impurity accumulation, the application of too off–axis ECRH or off–axis ICRH usually does not allow one to reach the same type of effects and thereby does not lead to the avoidance of central impurity accumulation [125, 207, 237]. Recently it was reported that in EAST the application of LHW also helps to avoid the central tungsten accumulation in dominant NBI heating plasmas [214, 238]. Integrated modeling is performed like in [233], and the results suggest an increase in the turbulent diffusion by increasing the electron temperature and a reduced inward neoclassical convection by decreasing the toroidal rotation.

The increased understanding of both the neoclassical and the turbulent transport processes in impurity transport, and particularly the clear identification of the critical ingredients leading to neoclassical central accumulation of heavy impurities, allows the integrated modeling of the impurity behavior to be performed, making use of a combination of theory–based models. Recent state–of–the–art examples of this activity, which apply completely theoretical models for the description of neoclassical and turbulent transport, are [195, 239]. The latter [239] is based on a fully integrated core–pedestal modeling framework, including a model for the pedestal structure and the application of neural–network accelerated models, also for light impurity turbulent transport. It makes use of the TGLF–NN turbulent transport model, also trained for impurity transport, coupling the transport solver TGYRO [240] with the STRAHL impurity transport code [241]. It compares predictions of this highly integrated modeling workflow with DIII-D experimental results and also includes predictions of the impact of increasing concentration of neon on the ITER baseline scenario. The former integrated modeling approach [195] adopts the JINTRAC [242] integrated transport simulator and also models the transport of heavy impurities like tungsten and nickel for comparison with JET experimental results and to optimize tungsten control with ICRH. It further develops a previous multi–channel integrated modeling effort [243] and applies the new modeling framework to the prediction of the time evolution of the tungsten behavior in JET. It combines the comprehensive description of neoclassical impurity transport provided by the NEO code [176], also applied in [239], with the turbulent transport, including impurities, of the QuaLiKiz [244] quasi–linear transport model. The turbulent transport model includes the impact of poloidal asymmetries [245, 246]. These transport simulations allow comprehensive descriptions in an integrated modeling approach of the physical effects which have been separately identified in previous dedicated investigations with the drift–kinetic and gyro–kinetic modeling of specific experimental observations. An example of the results obtained with this comprehensive integrated modeling approach is presented in figure 10. Following the approach developed in [239], more recent works have been dedicated to the prediction of tungsten impurity and radiated power density profiles in EAST plasmas [196, 233, 247] respectively considering an improved high–performance fully non–inductive discharge, a hybrid scenario discharge and a long–pulse steady–state H–mode discharge.

Zoom In Zoom Out Reset image size

At the same time, the computational requirements of the drift–kinetic code NEO as a coupled module in the framework of transport modeling are heavy and have motivated the development of reduced analytical models, which can closely reproduce the NEO results. An example is the FACIT model [179, 188], which has been also extended to include the impact of rotation at low collisionality [179].

A even more recent application [248] has been dedicated to the modeling of a pair of JET plasmas in which an increase of confinement is observed with neon seeding in comparison to the unseeded case. This includes the modeling of the neon density profiles by combining the QuaLiKiz and NCLASS for turbulent and neoclassical transport and also includes a comparison between QuaLiKiz neon impurity transport predictions with the results of the GENE gyrokinetic code.

Although not directly connected with actual transport mechanisms produced by collisions or turbulence, additional important effects on the impurity behavior are caused by the occurrence of MHD instabilities. These can be mainly classified in two groups. The first group is given by the impact of modes connected with the presence of the q = 1 surface, in particular sawtooth relaxations as well as ideal or resistive deformations of the magnetic equilibrium in the innermost region of the plasma, affecting the impurity behavior in the proximity of the magnetic axis. The second group is connected with the occurrence of neoclassical tearing modes (NTMs). Therefore, these are effects which impact the impurity behavior in a radial domain which is more usually located around mid–radius. Sawteeth, as well as fishbones, have been historically recognized to have a beneficial effect against central impurity accumulation, producing a flattening of the central impurity density profile and consequently a more off–axis redistribution of the profile of the radiated power density [192, 249, 250]. Consistent observations of this phenomenology have been also reported and analyzed more recently [251]. Moreover, deformations of the magnetic equilibrium in the proximity of the magnetic axis, as produced by the presence of long–living, saturated internal kink modes at the q = 1 surface, can significantly modify the impurity behavior, particularly through an interplay with central ECRH, leading to experimentally measured centrally hollow impurity density profiles [252, 253]. These central impurity holes cannot be explained by the usual combination of neoclassical and turbulent transport effects in the presence of an axisymmetric magnetic equilibrium [194, 253] and still await a physical explanation, likely connected with 3D effects of this helical magnetic core on the transport. The occurrence of NTMs is observed to impact the impurity behavior [193, 254] and experimental evidence is obtained that the tungsten accumulation can be initiated or accelerated by the appearance of the modes [254].

The threshold value of $R/L_T$ for turbulence onset is identified experimentally by reaching the condition of heat flux close to 0 inside a given radius, i.e. by localizing external power off-axis, or alternatively by performing a scan of heat flux by changing the proportion of on- vs off-axis heating and calculating the $R/L_T$ intercept for 0 flux, whilst the slope yields the stiffness. This implies that neoclassical ion heat transport is negligible, which is a good assumption for tokamaks outside the magnetic axis region [240]. To aid with the interpretation of the measurements against local gyrokinetics, it is useful to normalize to the natural dimensionless units of the system, i.e. gyroBohm units, e.g $q_\mathrm{i,e}^{gB} = q_\mathrm{i,e} (eBL^2)/ (nT^2\rho_\mathrm{ref})$, with n the reference density, T the reference temperature, L a reference length, and B the reference magnetic field. $\rho_\mathrm{ref}$ is the reference Larmor radius $\rho_\mathrm{ref} = (mT)^{1/2}/eB$, with m the reference mass, most often taken as the main ion mass, and e the electron charge. The specific reference values (e.g. $T_\mathrm i$ or $T_\mathrm e$ for reference temperature, major or minor radius for reference length, etc) may vary for different authors or codes, and care must be taken when comparing. If the heat flux scan is experimentally repeated with different values of one dimensionless plasma parameter (to be kept constant during each heat flux scan) then the threshold dimensionless parametric dependencies can be measured. This requires very careful dedicated transport experiments, which were performed on JET [295] and AUG [296]. An example of expected parametric trends can be seen in an analytical formula from [300], which in the flat density limit predicts dependencies of ion threshold on $T_\mathrm e/T_\mathrm i$ and $\hat{s}/q$ as

$$\begin{equation} R/L_{Ti}^\mathrm {crit} = 4/3 \left(1+T_\mathrm i/T_\mathrm e\right) \left(1+2 \hat{s}/q\right)\end{equation} \tag{ 5 }$$

where q is the safety factor and $\hat{s}$ the magnetic shear. More exact absolute predictions and multivariate parametric trends are found through carrying out gyrokinetic calculations, e.g. as in [301]. Figure 10 in [295] shows the ion threshold measurements at different $T_\mathrm e/T_\mathrm i$, in good agreement with linear predictions by the gyrokinetic code GKW [30], whilst figure 11 in [295] shows the dependence on $\hat{s}/q$, also found in good agreement with the gyrokinetic code GS2 [302]. Also figure 7 in [296] supports the $T_\mathrm e/T_\mathrm i$ dependence of the ion threshold predictions.

Particular attention has been given both in experiments and modeling to the stabilizing effect of the E × B flow shear. This introduces a shift in the effective ion threshold, as discussed e.g. in [270]. Early simplified treatment of this effect [303, 304] was based on the replacement in a transport model of the maximum linear growth rate $\gamma_\mathrm{max}$ by an effective growth rate $\gamma_\mathrm{eff} = \gamma_\mathrm{max} - \alpha_\mathrm E\gamma_\mathrm E$ where

$$\begin{equation} \gamma_\mathrm E = \left|\frac{RB_\mathrm p}{B_\mathrm T^2}\frac{\mathrm d}{\mathrm dr}\frac{E_\mathrm r}{RB_\mathrm p}\right|\end{equation} \tag{ 6 }$$

with $B_\mathrm p$ and $B_\mathrm T$ being the toroidal and poloidal fields, $E_\mathrm r$ the radial electric field and $\alpha_\mathrm E$ a constant. This ‘quench rule’ often provides a significant degree of E × B stabilization, although this can be reduced by kinetic electrons, TEM drive, and destabilizing PVG which depends on the magnetic geometry [305]. A general quench rule model would demand a parameterization of $\alpha_\mathrm E$ taking this into account. A more nuanced approach revolves around the finding that the E × B flow shear induces a shift in the peak of the radial wavenumber spectrum of electric potential fluctuations, and a reduction in their amplitudes. The occurrence of this shift has been recently measured directly in a dedicated experiment [306]. Such a spectral shift can be modeled analytically [307] and is now the paradigm for the E × B implementation in the TGLF model, with a generally reduced effect of stabilization compared to the previous implementation, in better agreement with experiments. In QuaLiKiz both parallel and perpendicular flow shears are taken into account through a sheared Doppler shift of the mode frequency, a modification to the diamagnetic frequency, and through inducing a shift in the ansatz Gaussian in the eigenfunction calculation, which leads to symmetry breaking and momentum transport [308].

In higher fidelity gyrokinetic codes such as CGYRO [277, 309] and GENE [310, 311], the implementation of the E × B shear can be taken into account through radial wavevector evolution methods. GENE uses a time-dependent wavevector shift scheme [312], which has been recently corrected to avoid spurious nonlinear couplings [313], whereas CGYRO uses a continuous wavevector advection scheme [277]. The impact of E × B shear has been shown to be negligible in nonlinear gyrokinetic simulations, for low-flux regimes with strong ITG stabilization due to e.m. effects [314]. For the same settings but taking the electrostatic (ES) limit, E × B stabilization is predicted. While conjectured to be related to predicted decreases in auto-correlation time and increased parallel mode structure in the e.m. case, the reason for this observation is still under investigation.

Measurements of ion or electron stiffness are obtained by performing a scan of the ion or electron heat flux by changing the proportion of on- vs off-axis heating and calculating the slope of the $q_\mathrm{i,e}^{gB}$ vs $R/L_{T_\mathrm{i,e}}$ curve. Alternatively, power modulation experiments allow an independent determination of stiffness, extracting from the heat pulse propagation the local derivative of the $q_{i,e}^{gB}$ vs $R/L_{T_{i,e}}$ curve, as described in detail in [292]. The first measurements of ion stiffness were performed in JET L-mode plasmas, thanks to well time and space resolved CX diagnostics for $T_\mathrm i$ and to the availability of ICRH power in the $(^3He)-D$ scheme, providing localized ion heating with short slowing down time [293–295]. This allowed both the heat flux scan and the modulation techniques to be used. On AUG, the heat flux scan was made by using on- and off-axis NBI [296], without modulation. In addition, ICRH in $(^3He)-D$ was applied, but only centrally to study the effects on local transport without heat flux scans or modulation [315]. On DIII-D, different levels of centrally localized NBI power with co- or balanced injection were applied in H-mode plasmas to study the effects on local transport with and without rotation, but no heat flux scans or modulation were performed [297, 298].

Figure 11(a) shows the results of the JET ion heat flux scan. A significant reduction of core ion stiffness at high NBI power was observed, leading to very peaked $T_\mathrm i$ profiles in the inner half of the plasma (figure 11(b)). The local slope of the $q_i^{gB}$ curve from modulation (indicated by the two segments on the encircled points) confirms the stiffness reduction. This effect was found to be enhanced by low magnetic shear ($\hat{s}$). Although in JET, with NBI heating, power and rotation are inherently coupled, this effect cannot be ascribed to increased rotation, which is mainly producing a threshold up-shift, partly compensated by the destabilization due to the PVG, as discussed in [316, 317]. Instead, an explanation was found in terms of nonlinear electromagnetic (e.m.) stabilization linked to thermal and suprathermal pressure [316, 317]; the associated complex physics is discussed in detail in [299] and will be summarized in section 4.2.3. These effects have been the object of intense theoretical studies in the last decade. Additional interesting JET results were the repetition of the ion stiffness measurements in Hydrogen, finding a much lower reduction of ion stiffness [318]. This was found due to differences in the fast ion content between H and D plasmas, inducing lower nonlinear e.m.stabilization in H, thus providing an explanation for the observed anti-gyroBohm dependence of the core transport in these discharges. A similar effect was found in H verses D H-mode plasmas in AUG [319]. Instead, in AUG L-mode plasmas with pure electron heating and no fast ions, ion stiffness was found independent of mass, with mass entering only in the energy exchange term between electrons and ions [320].

Zoom In Zoom Out Reset image size

The experiments on AUG and DIII-D support the same type of phenomenology. In AUG [296], moving NBI power from off-axis to on-axis leads to a non-monotonic $q_i^{gB}$ vs $R/L_{T_i}$ trend, interpreted as an increased ion stabilization for central NBI power. This seems less effective for higher $T_\mathrm e/T_\mathrm i$. Also on AUG, central ICRH in $(^3He)-D$ minority scheme leads to very peaked T_i profiles in the inner region [315], very similar to those observed on JET (figure 11(b)). On DIII-D the role of e.m. stabilization and of rotation were clearly identified by a NBI power scan with co- and balanced injection, the results of which are shown in figure 12 [297]. Here the shift between co- and balanced NBI points is indicative of the rotation effect, but in both sets we see the trend of decreasing $q_i^{gB}$ with increasing $R/LT_\mathrm i$ seen in JET and AUG. We remind that in figure 12, being a power scan for central heating and not a heat flux scan, we only see the tips of the curves shown in figure 11. These observations have been explained by gyrokinetic modeling as due to nonlinear e.m. stabilization associated with increasing thermal pressure, with little contribution from fast ions [298], which were instead a dominant contribution in the cases of JET and AUG experiments.

Zoom In Zoom Out Reset image size

Besides these specific transport experiments, a general experimental observation is that in the core of Hybrid plasmas, where $\hat{s}$ is low and pressure is high, very high $T_\mathrm i$ peaking is observed for very low $q_i^{gB}$, which is suggestive of very low ion stiffness, e.g. in JET [314, 321, 322] and AUG [323]. This has also been ascribed by gyrokinetic simulations to the same physics of nonlinear e.m. stabilization, with a small effect of E × B shear, which is instead more dominant in the external part of the plasma, where the nonlinear e.m. stabilization is reduced due to increased $\hat{s}$ and reduced thermal and suprathermal pressure [317, 321]. However, e.m. stabilization can still play a (diminished) stabilizing role also at outer radii, according to gyrokinetic simulations [299].

Finally, all the above experimental results on ion stiffness reduction concerned plasmas with moderate energy (∼100 keV) fast ions, as provided by present-day NBI and ICRH systems. An obvious question arises whether such stabilization effects would be present also in a reactor characterized by fast ions in the MeV range. Recent JET experiments using ICRH in a three-ion scheme (see [324]) revealed that stabilization could take place also in the presence of fast ions in the MeV range [325]. In spite of a variety of Toroidal Alfvèn Eigenmodes (TAEs) with different toroidal mode numbers excited in the central region, the ion heat transport was found improved with respect to a discharge with similar power and low energy NBI fast ions only. This encouraging observation, nevertheless, requires further verification and will, hopefully, be corroborated by new results from the JET 2021 D-T campaign, which will enable the investigation, in detail, of the effects of α particles on thermal transport, enhancing understanding towards the extrapolation to burning plasmas.

As discussed in [299], a number of linear mechanisms that stabilize ITG turbulence are known. Main ion dilution by fast or impurity particles [326, 327] and Shafranov shift stabilization [328–330] are well known basic electrostatic (e.s.) mechanisms. A newly discovered resonant linear e.s. mechanism is described in [283], based on wave-particle resonance between fast ions and the ITG mode when the fast ion magnetic drift frequency approaches the linear ITG frequency. The effect is more pronounced for ICRH than NBI fast ions, with an optimal fast ion temperature that depends on parameters such as the fast ion charge and vanishes for highly energetic fast particles [331]. This mechanism was found effective on JET with $(^3He)-D$ ICRH [332], and recent experiments on AUG succeeded in producing very peaked $T_\mathrm i$ profiles optimizing ICRH in $(H)-D$ minority [285]. Global GENE simulations in the inner half of the plasma show a narrow region $\rho_\mathrm{tor} = 0.2$–0.3 with full turbulence suppression leading to high T_i gradients, with localized shearing layers at its margins [285].

Linear e.m. effects through coupling of ITG to shear Alfvèn modes lead to ITG stabilization with increasing β, already without fast particles [333, 334], but significantly enhanced by the presence of fast particles increasing the drive of the shear Alfvèn modes [317, 335]. These mechanisms were, however, shown in [317] to be insufficient to explain the JET observations of figure 11. To match experiment, it was key to include in gyrokinetic simulations with GENE both fast particles and nonlinear e.m. effects [316, 317]. Figure 13 shows that most of the stabilization leading to the high $T_\mathrm i$ peaking of the JET high NBI case of figure 11(b) is due to the fast ion nonlinear e.m. effect. Also for the AUG results, both with NBI and $(^3He)-D$ ICRH, nonlinear e.m. stabilization by fast ions was found to be the dominant mechanism [296, 315]. In JET high power Hybrid plasmas a contribution to the ITG stabilization was found due to fast ions [314], with a non-negligible effect also in ITER simulations [336]. The physics behind such stabilization was investigated in detail in [337] using parameters of the JET high power case of figure 11(b). In figure 14(a), a stronger reduction of $q_i^{gB}$ with increasing β_e is found when including fast ions with respect to the case of pure thermal β_e. Two phases for the stabilization process are recognized, as shown in figure 14(b) for the case with highest β_e. In phase I, fast ions provide marginally linearly stable TAEs, which are then further nonlinearly excited by a cross-scale energy redistribution by the ITG modes, leading to a first degree of stabilization. In phase II, the fast ion modes transfer energy nonlinearly to ZFs, whose increase provides a major stabilization of the ITG turbulence. This in turn reduces the nonlinear TAE drive until an equilibrium is reached at lower transport level. We note that this mechanism concerns only ITG driven heat flux and has no impact on the electron heat flux driven by TEMs, for which ZF saturation is not key [338–340].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Fast ions are not a necessary ingredient for nonlinear e.m. stabilization, although they greatly enhance its effect. For the DIII-D experiments shown in figure 12, CGYRO simulations indicated that the main mechanism for ITG flux reduction was nonlinear e.m. stabilization due to thermal β. In figure 15(a) one can see a very large effect of nonlinear e.m. stabilization and a very small effect of rotational shear. Fast ions in this case add a minimal correction to the large thermal β stabilization, as discussed in [298]. Also in JET and AUG Hybrid scenarios, the dominant stabilization is by the thermal β, with fast ions adding an additional flux reduction (see figure 15(b) and all cases illustrated in [314, 321, 323, 341, 342]). In this case, the underlying mechanism is different from the one due to fast ions described above. Early nonlinear gyrokinetic simulations with Cyclone Base Case already showed a significant decrease of transport at finite-β, much larger than in linear calculations [343, 344] and quasilinear estimates [345]. Such increased nonlinear e.m. stabilization in simulations is associated with a relative increase in ZF coupling efficiency. In [345] it was shown that while there is no change to the fundamental saturation mechanism of three-wave nonlinear interactions (triplet coupling of stable ITG, unstable ITG and ZF) which lead to the rise of ZFs, the efficiency of this nonlinear coupling increases in the finite-β case. Analytically, it was possible to relate the coupling efficiency to a triplet correlation time τ [345, 346]. While an exact calculation of τ, including stable modes, is complex, a proxy for τ based on couplings between two unstable ITG modes and a ZF was suggested, which increases with β. This could be used to account for these effects in quasilinear estimates. In addition, also for the thermal β stabilization a similar mechanism as the one described for the fast ion case could be envisaged, but in which the role of TAEs could be played for example by Kinetic Ballooning Modes. Overall, a lot of progress has been made in the understanding of ITG nonlinear e.m. stabilization, and its importance is well demonstrated. Further studies of the relative role and parameter sensitivities of the various pathways for more effective saturation via ZFs in e.m. turbulence would be useful, including the understanding of why in some cases the manifestation of e.m. stabilization involves a threshold upshift rather than a stiffness reduction, as shown in [341].

Zoom In Zoom Out Reset image size

Even more work is still needed for the understanding of stabilization mechanisms when the energy of the fast ions is in the MeV range. Gyrokinetic flux-tube simulations of the JET 3-ion experiments mentioned in section 4.2.2 [325] showed linearly that ITGs are dominant and TAEs are unstable. A scan of the fast ion normalized pressure gradient, $R/L_{p_\mathrm{FD}}$, showed a strong stabilizing effect of the fast ions on the thermal ion diffusivities, leading to agreement with power balance, when the TAEs are near marginal stability, persisting also when TAEs are fully destabilized. Such stabilization is due to the formation of intense poloidally oriented and radially sheared ZF structures in the e.s. potential, with a complex interplay between the fast-ion driven TAEs and the nonlinearly generated ZFs. The experimental electron heat flux was more difficult to recover by simulations due to growth of the part driven by magnetic flutter. Nevertheless, by reducing within uncertainties $R/L_{T_\mathrm e}$ and $R/L_{p_\mathrm{FD}}$, conditions were found where TAEs were still unstable and both ion and electron heat transport matched power balance. These results, albeit requiring further verification, are nevertheless encouraging for future reactors. However, since for MeV-scale fast ions the validity of the local limit is questionable, ultimately, highly demanding global flux-driven simulations are necessary to both validate present-day local gradient-driven simulations, as well as shed light on the precise mechanisms at play, to assess whether these fast-ion related ITG stabilization mechanisms remain relevant in burning plasmas.

Although section 9 describes in detail the state of the art of quasilinear models for integrated scenario simulations, in this section we focus attention specifically on the capabilities of reduced models to reproduce ion critical gradient and stiffness measured by the above described experiments, and on the stand-alone comparison for identical input of reduced models with gyrokinetic simulations, which contain a more complete physics description. In fact, a large effort was devoted within ITPA T$\&$C to such validation of two reduced models, TGLF and QuaLiKiz, contributing to the release of a number of improved versions. TGLF now features 4 versions of the saturation rule: SAT0 [287], local in k_y, SAT1 [33], featuring a nonlinear upshift of the critical threshold and multi-scale interactions, SAT1-geo which added to SAT1 an improved geometrical treatment, and SAT2 [288, 347] with an improved collision model and a fit of the geometry factors to 3D CGYRO spectra. QuaLiKiz besides the original version [29, 246] now features a version with improved collisional operator [348].

The agreement of TGLF and QuaLiKiz on the critical gradient (which is determined by linear physics) is in general very good with respect to GENE and experiments, both in absolute value and in trends with $T_\mathrm e/T_\mathrm i$ and $\hat{s}/q$. Examples are reported in [299] (figures 36 and 37), in [332] and in [322]. One notable exception has been reported in [349], where the effect of $Z_{\mathrm{eff}}$ by light impurity injection was found to cause an excessive threshold increase in QuaLiKiz with respect to GENE, TGLF and experiment. The treatment of the E × B effects on thresholds in both models has already been discussed in 4.2.1 and is now adequately matching gyrokinetics and experiments.

Much more complex is the correct reproduction of ion stiffness by reduced models, since that involves nonlinear physics and all the effects described in section 4.2.3, many of which are not included in the saturation rules of TGLF or QuaLiKiz. Both models include the linear e.s. effects, i.e. dilution, α stabilization and also the resonant effect, although the latter requires including fast ions as a kinetic species and improved numerical resolution, not often feasible in integrated modeling. However, neither model includes the nonlinear e.m. effects. TGLF includes the linear e.m. effects whilst QuaLiKiz is purely e.s.. Therefore, we expect that both models overestimate ITG transport when nonlinear e.m. stabilization (from thermal and/or suprathermal pressure) is a main player.

For low power cases without significant nonlinear e.m. stabilization, comparison of TGLF, QuaLiKiz, GENE and experiment indicated that early TGLF versions (SAT0 and SAT1) had too low levels of ion stiffness [350], which triggered a revision leading to the new TGLF versions SAT1-geo and SAT2 with higher ion stiffness, more in agreement with GENE and experimental data. Due to their low inherent ion stiffness, SAT0 and SAT1 in high power cases with significant e.m. stabilization were in fact matching experiment without the need of including any nonlinear e.m. effect, at variance with GENE showing that they are crucial to match experiment. In other words, they were reproducing experiments well for the wrong reasons. This indicates that the predictive capabilities of these models are unreliable. The two recent versions SAT1-geo and SAT2 instead have a more correct ion stiffness level, but since they are still missing nonlinear e.m. effects, they reproduce present high power experiments worse. However, this is a necessary step toward achieving a complete model, since it leaves margins to add the missing stabilizing effects. QuaLiKiz ion stiffness is found in good match with GENE, putting the model on similar grounds for ion heat transport as the new TGLF versions. The QuaLiKiz version with improved collisional operator in fact has mainly changed the electron heat transport predictions (see [299]), as well as predictions of density peaking leading to improvements in comparisons against JET experiments.

Such outcomes of stand-alone validation simulations are in fact confirmed by integrated modeling. In [272] it is discussed how the two models reproduce well H-mode baseline scenarios, but tend to fail when applied to advanced scenarios, with a significant under-prediction of T_i and density peaking in the inner half-radius, where the missing nonlinear physics is important due to high β and a broad region of low $\hat{s}$. As a provisional remedy for integrated modeling of such scenarios, ad-hoc models have been developed for both QuaLiKiz and TGLF, as discussed in [299]. For QuaLiKiz in JINTRAC [242], the input $R/L_{Ti}$ is decreased by the local ratio of thermal to total pressure [94, 195], since in NBI and ICRH heated discharges the degree of e.m. stabilization correlates strongly with this ratio. This proved successful in modeling high power JET D discharges for subsequent extrapolation to DT. For TGLF in ASTRA [351], a local increase of the concentration of an artificial impurity, which stabilizes ITG turbulence, proportional to the fast ion pressure and dependent on $\hat{s}$, was applied successfully to simulations of advanced scenarios in AUG [352]. It is clear, however, that the development of a theory-based e.m. stabilization model for quasilinear transport models should be a high priority research item for the near future. As discussed in [299], initial attempts to develop a proxy for the triplet coupling efficiency, related to ITG saturation with ZFs, were reported in [353]. This proxy can provide an approximate correction to quasilinear model saturation rules. However, the impact of fast-ions on the proxy and the role of subcritical Shear Alfvén modes has not yet been investigated. At the moment, such models have not been introduced into either TGLF or QuaLiKiz.

While the core electron heat transport has been historically interpreted as mostly due to ITG or TEM micro-instabilities [358, 360, 368], evidences are piling up indicating that ETGs could also play a key role in some cases. Early theoretical indications were given of the importance of the ETGs in setting the $T_\mathrm e$ peaking under particular experimental conditions for cases with dominant electron heating [355, 369], due to the formation of radially extended structures (‘streamers’) that could lead to a very effective cross-field electron heat transport. We recall that ETGs carry negligible ion heat and particle transport. But only in the last decade stronger arguments have been obtained which seem to support this picture. These arguments are based on the careful comparison of the experimental $q_\mathrm e$ and its slope versus $R/L_{Te}$ with the results of nonlinear gyrokinetic simulations, as well as with reduced quasilinear models. Two main strategies are adopted: the computationally cheapest consists of comparing the experiment with ion-scale gyrokinetic simulations ($k_{\theta}\rho_\mathrm i \lesssim 1$, where k_θ is the poloidal wave number and $\rho_\mathrm i$ is the ion Larmor radius), testing if the ion-scale contributions to $q_\mathrm e$ are sufficient to explain the experiment. If there is a gap between the ion-scale gyrokinetic $q_\mathrm e$ and the experimental value, beyond reasonable uncertainty propagation of the simulation input parameters, then ETGs could be invoked as an explanation for the missing flux. The second strategy consists in performing multi-scale gyrokinetic simulations ($k_{\theta}\rho_\mathrm e \lesssim 1$, where $\rho_\mathrm e$ is the electron Larmor radius), covering both ion and electron scales, and directly comparing the results with the experiment.

The first approach, despite the drawback of only giving a necessary but not sufficient condition for the ETGs to contribute to the experimental $q_\mathrm e$, is extremely computationally faster than the second, since a full-physics nonlinear multi-scale gyrokinetic simulation of experimentally relevant cases nowadays still usually requires tens of millions CPU hours to run. The first strategy is used for example in [361, 370]. In [370] this method shows that the ETGs could play a role for Alcator C-Mod L-mode discharges where the low-k TEMs are not dominant and instead the ITG are the stronger ion-scale instability. Since $q_\mathrm e$ is underpredicted by ion-scale GYRO gyrokinetic runs for those cases, while the experimental $q_\mathrm i$ is explained, the ETGs are invoked as possibly responsible for the missing $q_\mathrm e$. In [361], JET pulses are considered. In figure 14 of [361], the comparison of gyrokinetic GKW nonlinear ion-scale $q_\mathrm e$ vs $R/L_{Te}$ scans with the experiment is presented at mid-radius, showing that GKW is not able to explain the larger experimental fluxes, for $R/L_{Te} \gt 9$. These indicates a possible role for ETGs in setting the $T_\mathrm e$ core profiles.

The second strategy is used instead in [362, 364, 371, 372]. In [362, 371], the same experimental Alcator C-Mod pulse that had been considered in [370] is analyzed, performing (in [362]) $R/L_{Ti}$ and $R/L_{Te}$ scans of both ion-scale, electron-scale and multi-scale nonlinear GYRO simulations, showing that multi-scale simulations are required to match experimental ion and electron heat fluxes and electron profile stiffness. In particular, this kind of analysis allows the observation of a rich variety of cross-scale interactions between ion and electron scales, which are evidenced by the comparison of the different kinds of gyrokinetic runs, covering different regions of the k_θ spectrum. In figure 16 (figure 3 from [362]), representing the results of three multi-scale simulations with a scan in $a/L_{Ti}$ at fixed $a/L_{Te}$, it is shown that effects of both ion scales on electron scales and electron scales on ion scales are present. In figure 16(B) the reduction of the efficiency of ZF shear generation of the ITG-dominant ion scales with decreasing $a/L_{Ti}$ lets the ETG streamers develop, leading to an increase of the electron-scale fraction of $q_\mathrm e$ (green line) with decreasing $a/L_{Ti}$, with a maximum effect that is obtained when the ITGs are marginally unstable. This is consistent with what was observed in [373, 374]. On the other hand, a back reaction of electron scales on ion scales is observed, with enhanced ion heat flux (blue line in figure 16(A)) and ion fraction of $q_\mathrm e$ (yellow line in figure 16(B)) with respect to the result of ion-scale simulations (red line). Figure 4 in [362] shows that this second effect is enhanced with increasing ETG drive $R/L_{Te}$, and that only accounting for ETGs it is possible to recover the experimental electron stiffness. Indeed, high-k ETGs damp low-k ZFs, and the reduced ZFs result in enhanced ITGs, especially for marginally unstable ITGs. Both effects have been previously observed in [375], where a slightly different mechanism was deduced for ETG quenching by ion scales, with ETG streamers suppressed by ITG eddies instead of ZFs.

Zoom In Zoom Out Reset image size

An additional strong dependence of these cross-scale couplings on the E × B shearing rate $\gamma_{E\times B}$ is evidenced in [372], where multi-scale GYRO runs are performed for a DIII-D ITER baseline scenario with dominant electron heating and low input torque. Whilst noticeably observing the same kind of cross-scale effects of [362] for values of γ_E that are smaller than the experimental one, at experimental γ_E the GYRO multi-scale run predicts a complete stabilization of ion scales, due to the formation of a state dominated by ZFs, for values of the logarithmic gradients that are not consistent with ITG stability if only the ion scales are considered. This shows that through modifications of γ_E (with increasing γ_E for this case) the electron-scale$\rightarrow$ion-scale cross-scale effect can be reversed compared to [362] and the ETGs can suppress the low-k turbulence instead of enhancing it (see figure 8 in [372]). In particular, for an intermediate γ_E value, the experimental fluxes can be matched with a similar proportion of $q_\mathrm e$ coming from ion and electron scales (see figure 9 in [372]), indicating that the electron heat transport for this ITER-relevant case could be very intrinsically multi-scale in nature. Following these analyses, in [364] the results of a nonlinear multi-scale GENE [310] simulation are shown for a representative pulse extracted from JET C-wall dedicated electron heat transport experiments in L-mode plasmas, where different levels of NBI and ICRH power (on- versus off-axis, using mode conversion to provide pure electron heating) were applied, to vary the ETGs linear threshold through the τ parameter. Indeed, for those shots, the dependence of the critical $R/L_{Te}$ for the onset of the electron heat transport (derived from experimental data using a critical gradient model [376]) on τ is found consistent with the ETG threshold (equation (8)). The results of the nonlinear analysis indicate that the ETGs are needed to explain the experimental $q_\mathrm e$, but the ETG flux provided by the multi-scale GENE simulation is still not sufficient to match the experimental q_e (see figure 6 in [364]). This ‘multi-scale’ shortfall could be due to the fact that the impurities have not been accounted for in the runs for lack of computational resources, and accounting for them could have resulted in more ITG stabilization with weaker stabilizing cross-scale effect of ITGs on ETGs. However, this work represents a first indication of a possible role of ETGs in JET plasmas.

Inspired by these discoveries, new dedicated experiments have been performed during the last years on the TCV [366], AUG [296] and JET [365] tokamaks. The results obtained on the three machines are collected and compared in [367]. The same experimental framework has been considered for the different machines. Steady state electron heat flux scans ($q_\mathrm e^{gB}$ vs $R/L_{Te}$), obtained by varying the electron heating radial deposition, have been performed to estimate the electron temperature stiffness and test the possible presence of an ETG wall, also varying the ion/electron heating power ratio since it impacts the ETG linear threshold. When possible (TCV and AUG), an independent determination of the stiffness has been pursued by a perturbative analysis (based on ECRH modulation or on the study of the propagation of the sawtooth). The experimental results indicate that a high electron stiffness consistent with ETGs is compatible with the cases with $T_\mathrm e\sim T_\mathrm i$ and sufficiently high $R/L_{Te}$, in agreement with the linear theory. These promising cases have been studied with linear, nonlinear ion-scale, and multi-scale gyrokinetic flux-tube runs with the GENE code, at mid-radius. Impurities have been neglected for AUG for lack of computational time and the E × B shearing has only been considered for the TCV cases with NBI (negligible rotation at the radius of analysis for JET and much smaller γ_E for AUG with respect to TCV). Also the fast ions have only been accounted for when simulating the TCV cases with NBI, while they have been neglected for AUG (from NBI) and JET (from NBI and ICRH), due to the fast ion density fraction $n_\mathrm{FI}/n_\mathrm e$, which is more than 10 times smaller for AUG and JET than for TCV. Moreover, only ion-scale simulations have been performed for the TCV cases, given the importance of possible global effects in the smaller size TCV tokamak (a full-physics global multi-scale run was not computationally feasible). The ion-scale gyrokinetic runs indicate that the AUG and JET ion scales are ITG-dominant, while for TCV the TEM are the main unstable mode. Moreover, it is not possible to match both the $q_\mathrm e$ levels and the electron stiffness with the ion-scale simulations. Only the flux levels can be matched for JET by varying the parameters within error bars, for AUG $q_\mathrm e$ is slightly underestimated, while for TCV $q_\mathrm e$ is strongly underestimated because of a synergic effect of fast ions and E × B shearing (both due to NBI), stabilizing the TEMs. As a consequence, ETGs could be invoked for all the TCV, AUG and JET cases as responsible for the disagreement. Finally, the multi-scale runs (AUG and JET) give the following picture: for AUG explain both the experimental flux levels and stiffness as possibly coming from ETGs (see figure 17(a)), where the results of the GENE multi-scale runs (black stars: neglecting impurities, red stars: taking into account impurities, lumped into an effective carbon impurity with experimental $Z_\mathrm{eff} = 1.4$) are compared with the experimental stiffness, represented in magenta, and with the experimental steady state scan, represented by circles. The results of the AUG GENE ion-scale and multi-scale simulations with impurities at $R/L_{Te} = 11$ are new. The corresponding q_e spectra, compared with the results of the simulations with $Z_\mathrm{eff} = 1$, are shown in figure 18. A stabilizing ${\sim}{-}20\%$ effect of impurities is found on q_e, approximately equally divided between ion-scales ($k_y \lt 1$) stabilization and electron-scales ($k_y \gt 1$) stabilization. These results, however, have to be taken with caution, since a relatively small statistics ($\Delta t c_s/a\sim 44$, with c_s the ion sound speed and a the edge unnormalized toroidal radius). For JET, the gyrokinetic fluxes do not explain both the experimental flux levels and stiffness (see figure 17(b)), where GENE multi-scale is shown in blue and the steady state scan data in black), with ETGs that seem to play a minor role.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The stabilizing effect of collisions and electromagnetic effects, which are stronger for the AUG case, could be invoked to partially explain the larger role of ETGs that is observed for AUG. Finally, the high electron stiffness that is observed for the JET case, made impossible a fine tuning of the ITG drive $R/L_{Ti}$ within error bars close to ITG marginal stability. One then cannot exclude the possibility of finding an optimal value of $R/L_{Ti}$ that could allow to match the experimental $q_\mathrm i$ for smaller $R/L_{Ti}$ than the one needed to match $q_\mathrm i$ with ion-scale runs, due to electron scales enhancing the ITGs through cross-scale effects. This would allow the ETGs to have a larger impact on $q_\mathrm e$.

In order to further investigate the impact of the impurities on the AUG results, and also test their effect on the JET results, without performing computationally unfeasible gyrokinetic multi-scale $R/L_{Te}$ fine scans, TGLF quasilinear simulations have been run. The results, collected in figure 17, indicate that while for AUG it is possible to explain the experimental data with or without considering the impurities, for JET it is only possible to get close to explaining the experimental data when neglecting heavy impurities, which heavily impact the position of the ETG wall. In particular, for the AUG case, the TGLF results agree with GENE multi-scale ones predicting a moderate impact of impurities on the ETGs, and in particular on the position of the ETG wall. TGLF strongly overestimates the stabilizing effect of impurities at fixed $R/L_{Te}$ with respect to GENE, and it seems to also overestimate the $R/L_{Te}$ upshifting of the ETG wall when adding impurities, even if this is difficult to quantify, considering that only a single GENE multi-scale simulation with impurities is available ($R/L_{Te} = 11$).

Up to this point, the picture that emerges is that ETGs could impact $q_\mathrm e$ when $T_\mathrm e\sim T_\mathrm i$ and $R/L_{Te}$ is sufficiently high. The recent work [377], consisting in GKV gyrokinetic multi-scale simulations of a high electron temperature tokamak plasma with electrons, deuterium and tritium fuel with helium ash, mimicking future experiments on burning plasma, indicates that ETGs could impact the electron heat transport up to $T_\mathrm e/T_\mathrm i\sim3$. ETGs, through cross-scale interactions, could suppress near marginal TEMs, while the shearing of ETGs by TEM turbulent eddies can suppress ETGs (see figure 4 in [377]). This is different from what was observed in [378], where TEM-driven ZFs were at work to suppress ETGs, and similar to [375], where ETGs are quenched by ITG eddies.

Finally, since in addition to ITG and TEM modes, micro-tearing modes (MTM) can be destabilized with finite β at ion scales, it is worth mentioning two works that analyze cross-scale interactions between ETGs and MTMs. The first [379] shows that MTMs are suppressed by ETG turbulence, since the fine-scale E × B flows of ETGs distort the localized current-sheet structures of MTMs. The second [380] shows that MTMs and ETGs can coexist and interact in a tokamak pedestal.

Simplified criteria, based on the results of linear gyrokinetic simulations, are often invoked to predict the nonlinear behavior of ETGs, guiding the interpretation of the experiments and the planning of computationally heavier multi-scale nonlinear runs. The two most commonly used heuristic criteria can be applied once the linear spectra of the growth rates $\gamma(k_y)$ of the most unstable mode are known. The first states that when the ratio $\gamma_\mathrm{ETG}^{\mathrm{max}}/\gamma_\mathrm{ITG}^{\mathrm{max}}$ of the maximum linear growth rate in the electron-scale and ion-scale k-ranges is larger than $\sqrt{m_\mathrm i/m_\mathrm e}$ (here referred to as ‘criterion 1’ for brevity), significant levels of electron heat flux are driven by ETG turbulence [310, 370]. The second, based on a ZF mixing model [381], states that ETGs should nonlinearly impact $q_\mathrm e$ when $(\gamma/k_y)_\mathrm{ETG}^{\mathrm{max}}/(\gamma/k_y)_\mathrm{ITG}^{\mathrm{max}} \gt 1$ (‘criterion 2’). The second criterion is also recovered by [382], with a different model based on accounting for the effect of parallel-to-the-field E × B shear. The latter criterion reduces to the former in the approximation $k_\mathrm{ETG}^{\mathrm{max}}\sim\rho_{Te}^{-1}$, $k_\mathrm{ITG}^{\mathrm{max}}\sim\rho_{Ti}^{-1}$, with $\rho_{Te,i}$ the electron and ion thermal Larmor radii. While in most of the cases the predictions of these criteria are confirmed a posteriori by the results of nonlinear multi-scale gyrokinetic simulations (e.g. [296, 362, 364, 367, 370] (AUG)), for some recent studies this does not always hold, and even if a simplified criterion is invoked to predict strong ETGs, sometimes the gyrokinetic multi-scale runs do not confirm the picture, predicting weaker or stable ETGs (e.g. [299, 365, 367, 383] (JET)). In [383], studying a DIII-D L-mode plasma in the near-edge, it is found that even if criterion 1 is satisfied and criterion 2 is marginally satisfied by the ITG-dominant ion scales and ETG-dominant electron scales (see figure 6 in [383]), nonlinear multi-scale GENE runs show little impact of ETGs on q_e (see figure 10 in [383]). The same behavior is recovered for [365, 367] (JET), as can be seen comparing figure 2(c) with figure 5(b) and figure 6(b) in [367]. Finally, in [299], a high-performance hybrid H-mode scenario is investigated, comparing predictive transport modeling with gyrokinetic GENE and GKW simulations at $\rho_\mathrm{tor} = 0.65$. The results at fixed radius indicate that, even if criterion 2 is consistently met (see figure 19(top), adapted from figure 11 in [299]), the ETG contribution to $q_\mathrm e$ is negligible (see figure 19(top), adapted from figure 13 in [299]). There are few caveats in applying this gyrokinetic result to the direct interpretation of the experiment, such as the simplified $s-\alpha$ geometry adopted in the runs, the $Z_\mathrm{eff} = 1$ assumption, and the impossibility to perform a fine $R/L_{Ti}$ scan down to ITG marginal stability to spot possible electron-scale $\rightarrow$ ion-scale cross-scale effects allowing ETGs to emerge, with multiscale effects raising q_i up to power-balance consistent values. However, this is a clear case where the simplified rules of thumb based on the linear eigenvalue spectra are not applicable. The authors thus emphasize that ‘there are hidden variables at play which complicate the existence criteria of significant ETG flux beyond simple rules-of-thumb’, and they also point out that ‘exceeding the nonlinear threshold for ETG turbulence does not necessarily correlate to significant ETG flux, which also requires the existence of significant ETG streamers—elongated structures’, which are absent in this case.

Zoom In Zoom Out Reset image size

The existence of reliable reduced models that include cross-scale effects is essential in view of allowing integrated modeling of scenarios where the transport is manifestly multi-scale in nature, due to the impossibility of coupling computationally heavy multi-scale gyrokinetic modules to transport codes. Sometimes, the implementation of the cross-scale effects in reduced quasilinear models is based on the simple criteria that have been shown above. Indeed, in QuaLiKiz only the quenching of ETGs due to ITG ZFs is included, through the $(\gamma/k_y)_\mathrm{ETG}^{\mathrm{max}}/(\gamma/k_y)_\mathrm{ITG}^{\mathrm{max}} \gt 1$ criterion (ETG are turned on when the criterion is approximately met). TGLF, on the other hand, employs a more refined zonal-flow mixing model [33]. However, the inverse electron scale $\rightarrow$ ion scale effects like the enhancement of ITGs due to ETGs [362, 375](or quenching, depending on γ_E, see [372]) are still not retained by the present quasilinear codes (including QuaLiKiz and TGLF). This could be, for example, one of the reasons that prevent TGLF from obtaining agreement with the experiment in [365, 367] (JET) (see figure 17(b)). In order to validate the simple criterion 2 versus the more refined TGLF multi-scale model, an effort has been pursued in [385] to compare TGLF (used within the VITALS [386] framework) predictions with experiment for twelve discharges from the Alcator C-Mod and ASDEX Upgrade tokamaks. A multi-channel comparison (matching experimental $q_\mathrm e$, $q_\mathrm i$, the percentage electron temperature fluctuation $\delta T_\mathrm e/T_\mathrm e$ and the perturbative electron heat diffusivity $\chi_\mathrm e^\mathrm {pert} = -(1/n_\mathrm e)\mathrm dq_\mathrm e/\mathrm d(\nabla T_\mathrm e\cdot \hat{r})$, varying $R/L_{Te}$, $R/L_{Ti}$, $R/L_{n}$ and $Z_\mathrm{eff}$ within experimental error bars) is performed, testing when ion scales are sufficient to explain the experiment and comparing the outcome with the predictions of two rules of thumb: the criterion 2 and a different rule that is met when at ion scales ITG and TEM are balanced (this criterion reads ${\mathrm{exp}}[-((\gamma/k_y)_\mathrm{TEM}^{\mathrm{max}}/(\gamma/k_y)_\mathrm{ITG}^{\mathrm{max}}-1)^2/0.09] \gt 0.5$ and will be here referred to as ‘criterion 3’ for brevity). The results show that when all the channels are considered, the multi-scale TGLF model is always able to explain the experiments, and that sometimes when retaining only ion scales and just looking at matching $q_\mathrm e$ and $q_\mathrm i$, it is possible to explain the experiment by fortuitous agreement, since $\delta T_\mathrm e/T_\mathrm e$ and/or $\chi_\mathrm e^\mathrm{pert}$ are not matched (see table 2 in [385]). Regarding the validation of the rules of thumb, it is found that, when both criteria 2 and 3 are fulfilled, the ion scales are not sufficient to explain the experiment and the transport is multi-scale in nature (see figure 2 in [385]). In particular, the fulfillment of criterion 3 means that at ion scale ITG and TEM are balanced, and just a small multi-scale interaction of ITGs or TEMs with ETGs could transform the overall character of the transport, unbalancing ITGs and TEMs, and impacting the fluxes, fluctuations and diffusivities. A remark on these results has to be made since, as it has been pointed out at the end of the previous paragraph, some cases where criterion 1 is met, but ETGs are shown not to impact $q_\mathrm e$ in gyrokinetic multi-scale simulations, have been found. This can be partially due to the fact that the multi-scale model that is implemented in TGLF is simplified with respect to the gyrokinetic formalism, that is potentially able to include all the physics that matters. Moreover, the parameter region that has been explored in [385], as well as the amount of retained physics, is limited (e.g. in [299] it is pointed out that all the cases in [385] are in a regime with $q_\mathrm e \gt q_\mathrm i$ and TEMs that fill a spectral gap between ion and electron scales, while in [299] a regime with $q_\mathrm e \lt q_\mathrm i$ and no spectral gap is investigated). However, more multi-channel analyses of this kind should be pursued, also comparing with gyrokinetic runs for selected cases, in order to improve the understanding of these phenomena.

The local delta-f approximation of the gyrokinetic equation retains only the lowest order terms in a $\rho_*$ expansion ($\rho_* = m_\mathrm i V_\mathrm{th}/L_{\perp} e B_\mathrm t$), where $V_\mathrm{th}$ is the thermal velocity, $L_{\perp}$ is the characteristic perpendicular scale length, generally represented by the minor radius a in the plasma core, and $B_\mathrm t$ is the toroidal magnetic field. Since rotation is a signed quantity, the symmetry properties of the gyrokinetic equation imply that the radial flux of parallel momentum due to the fluctuating E × B drift is zero, unless the parallel symmetry is broken [27, 416]. That is to say, when integrating over a flux surface, positive and negative velocities are equally probably and balance each other out, unless an asymmetry is present that favors one over the other. A systematic examination of the gyrokinetic equation and parallel symmetry breaking mechanisms [356, 422] yields exactly four contributions to the toroidal momentum flux that are lowest order in $\rho_*$. These four contributions, together with the next order flux in $\rho_*$, $\Pi_*$, can be formally written as:

$$\begin{align} \frac{\Pi_\phi}{n_\mathrm im_\mathrm iR_0} = -\chi_\parallel v^{^{\prime}}_\parallel - s_j\chi_E \gamma_{E} + R_0V_\phi \Omega_\Phi + C_\mathrm{FS} + \frac{\Pi_*}{n_\mathrm im_\mathrm iR_0}\end{align} \tag{ 10 }$$

where $v^{^{\prime}}_\parallel$ is the flux-surface-averaged parallel flow shear, $\gamma_E = \frac{1}{B}\frac{\partial^2\Phi}{\partial r^2}$ is the E × B flow shear and s_j the sign of the plasma current with respect to the positive toroidal direction. In the limit consistent with the local delta-f framework, i.e. to lowest order in $\rho_*$, the plasma flow is purely toroidal [119]. As a consequence, the parallel flow shear and the E × B shearing rate are coupled and proportional to the toroidal flow shear: $v^{^{\prime}}_\parallel\sim\frac{RB_t}{B}\frac{\partial \Omega_\Phi}{\partial r}$ and $\gamma_E\sim -s_j\frac{RB_p}{B}\frac{\partial \Omega_\Phi}{\partial r}$. Their contributions to the momentum flux can thus be grouped into a diffusive flux with $\chi_\phi = \chi_\parallel \frac{B_t}{B}-\chi_E\frac{B_p}{B}$. The next two contributions are the Coriolis pinch $R_0V_\phi$ and the residual stress driven by any up–down asymmetry in the magnetic equilibrium, $C_\mathrm{FS}$. The four lowest order contributions drive a turbulent momentum flux that scales with the square of the fluctuating fields and is, therefore, proportional to $\rho_*{}^2$ in the local delta-f approximation.

To next order in $\rho_*$, a number of additional contributions arise. Similarly to heat transport, these next order contributions can change the magnitude of the diffusive and convective momentum fluxes in smaller, low $\rho_*$, devices. In larger machines, however, they are not expected to significantly impact either the diffusion or convection. In contrast, for the residual stress, these next order $\rho_*$ contributions can be large and strongly impact the intrinsic rotation. This is because strong local gradients in kinetic profiles and turbulence properties can develop, for example in edge pedestals, that can cause higher order $\rho_*$ terms to be of the same magnitude or even larger than the up–down asymmetry residual stress. In that case, residual stress contributions of higher order in $\rho_*$ can play a determining role in setting the intrinsic rotation profile. In addition, the $\rho_*$ scaling of $\Pi_*$ is not yet firmly established and it can potentially compete with lowest order contributions even at low $\rho_*$ values. In the following, we review the four lowest order contributions to the toroidal momentum flux and two next order residual stress contributions deemed to be particularly important.

The first work that showed that the ion heat diffusivity, $\chi_\mathrm i$, and momentum diffusivity driven by the parallel flow shear, $\chi_\parallel$, were expected to be similar was based on a simple model of ITG turbulence [404]. This provided an explanation for the experimental similarity between the rotation and ion temperature profiles. The result can be expressed more succinctly as a Prandtl number of order unity, where the Prandtl number is defined as the ratio of the ion-momentum to heat diffusivities, $\textrm{Pr} = \chi_{\phi}/\chi_\mathrm i\sim\chi_\parallel/\chi_\mathrm i$. This result was reinforced by a later study using a linear gyrokinetic model [356], which also showed that, for ITG turbulence, $\chi_\parallel/\chi_\mathrm i$ is of order one with only a limited variation as a function of plasma parameters. This was, itself, explored further using nonlinear simulations with adiabatic electrons in [423, 424] and including kinetic electrons in [425–427], all indicating that $\chi_\parallel/\chi_\mathrm i$ is of order unity.

In Strintzi et al [428], gyrokinetic calculations of the momentum diffusivity were compared to those of a fluid model that reproduced the gyrokinetic model as long as the terms associated with the particle drift in a curved magnetic geometry were maintained. This work also used linear gyrkokinetic simulations of ITG turbulence with kinetic electrons to show that $\chi_\parallel/\chi_\mathrm i$ varies between only 0.7 and 1.2 over a wide range of plasma parameters. In addition, $\chi_\parallel/\chi_\mathrm i$ was predicted to depend only weakly upon plasma β_e, magnetic shear, and $R/L_{T_\mathrm i}$. The behavior of the momentum diffusivity was also explored for tokamak plasmas dominated by TEM turbulence in [429] via quasilinear simulations and, in contrast to ITG turbulence, the momentum diffusivity was found to be larger than the ion heat conductivity ($\chi_\parallel/\chi_\mathrm i \gt$1) at low normalized electron density gradient, $R/L_{ne}$. With more peaked density profiles, the ratio was again predicted to be of order one, similar to the ITG results.

A diffusive momentum flux driven by E × B flow shear was first identified in [430] and further explored in [431–434]. For pure toroidal rotation and positive magnetic shear, it always tends to oppose the diffusive momentum flux driven by the parallel flow shear, i.e. $\chi_E \gt 0$, making the total momentum diffusivity, $\chi_\phi = \chi_\parallel \frac{B_\mathrm t}{B} - \chi_E\frac{B_\mathrm p}{B}$, lower and thereby reducing the Prandtl number. For small toroidal rotation shear, typical of intrinsic rotation, a Prandlt number as low as 0.5 can be obtained when including the χ_E contribution [433]. For negative magnetic shear, the situation reverses with $\chi_E \lt 0$ [433]. In contrast to $\chi_\parallel/\chi_\mathrm i$, which is mostly independent of $v_\parallel^{^{\prime}}$, $\chi_E/\chi_\mathrm i$ decreases with increasing γ_E making its contribution to the Prandtl number smaller for higher toroidal rotation shear [425, 433].

The experimental results on the magnitude of the Prandtl number are presented in section 5.3.1. However, as the Prandtl number is not expected to vary significantly, at least for ITG turbulence, it is difficult to provide robust experimental validation of the theory. The differences observed in the gyrokinetic predictions between ITG and TEM plasmas, on the other hand, provide a potential avenue by which some validation becomes possible and should be actively pursued by experiment. Complete validation of our theoretical picture will require self-consistent modeling of the rotation profile (see section 5.2.5) including not only the diffusive term, but also contributions from the convective and residual stress momentum fluxes, presented in sections 5.2.3 and 5.2.4, respectively.

The convective momentum flux usually refers to the total convective flux present in the plasma and includes two terms: a radial particle flux that carries some momentum and a true convective momentum flux, better known as the Coriolis pinch. The identification of mechanisms capable of generating a strong inward convection was particularly important in completing our theoretical picture of momentum transport and helping to explain the experimental observations. The role of the centrifugal force on momentum transport has also been considered, but is expected, under most circumstances, to be small [435]. However, for high beta plasmas in NSTX, where MTMs, hybrid ITGs, and kinetic ballooning modes (KBMs) were found to be the dominant instabilities, the momentum flux from the centrifugal force was significantly larger than that from the Coriolis force [436], showing that there are circumstances in which this term cannot be neglected. The formalism behind the centrifugal momentum flux is not discussed here, but can be found in [435].

The physics behind the Coriolis momentum flux is discussed in detail in a number of publications [27, 416, 435, 437–441] and is only briefly summarized here. The Coriolis pinch can be understood by considering a tokamak plasma with a rigid body rotation given by $v_{\phi} = R\Omega_{\phi}$. A particle moving parallel to magnetic field lines has a velocity component perpendicular to the toroidal direction and, therefore, feels a Coriolis force, $\vec F_C = 2m_\mathrm iv_{\parallel}\mathbf{b}\times \boldsymbol{\Omega_{\phi}}$ with b the unit vector in the direction of the magnetic field. As with any force within a magnetic field, this generates a perpendicular particle drift that enters into the gyrokinetic equation similarly to other drifts. This drift is vertical at the LFS mid-plane and adds to the magnetic (curvature and $\nabla B$) drift for $v_\parallel \gt 0$ (or reduces it for $v_\parallel \lt 0$). Therefore, for a given fluid density perturbation, the part that flows parallel to the field will drift vertically faster than that anti-parallel to the field. This means that the fluid element with $v_\parallel \gt 0$ will have a phase shift with respect to a potential perturbation that is different from the phase shift of the fluid element with $v_\parallel \lt 0$. This corresponds to a different radial transport by $\mathbf{E}\times\mathbf{B}$ convection of the two fluid elements and a net transport of momentum, which is inward in most cases and increases with the density gradient. The parallel velocity dependence of the Coriolis drift generates a coupling between density and parallel velocity perturbations. This is the underlying reason for the dependence of the Coriolis pinch on the density gradient. Incidentally, this exact same coupling is also responsible for the generation of roto-diffusion: a convective term relevant to impurity particle transport, see section 3.2.3.

The net momentum transport generated by the Coriolis pinch depends on the parallel mode structure [441]. If the mode structure is assumed to be independent of the toroidal rotation, as in simple fluid models, the predicted Coriolis pinch is maximized resulting in a strong dependence of the pinch on the normalized electron density gradient, $R/L_{n_\mathrm e}$. If the impact of toroidal rotation on the mode structure is retained, in the limit of adiabatic electrons, the mode structure adjusts so that the net momentum driven by the Coriolis pinch is zero when averaged over a flux surface. In the presence of trapped electrons, or other effects that make the electrons less adiabatic, the compensation becomes imperfect, resulting in a finite, non-negligible pinch, explaining why the strength of the Coriolis pinch depends upon the fraction of trapped particles. Gyro-kinetic simulations that include the impact of the parallel mode structure also show a strong, albeit weaker, dependence on $R/L_{n_\mathrm e}$ as well as dependences on the magnetic shear, safety factor, and ITG [437] as shown in figure 20.

Zoom In Zoom Out Reset image size

Most gyrokinetic simulations that examine the Coriolis pinch considered ITG dominated plasmas. The existence and behavior of the Coriolis pinch in TEM dominated plasmas was explored in [429], which showed that a finite pinch is predicted under most conditions in the presence of TEM turbulence. For both ITG and TEM turbulence, in standard aspect ratio tokamaks at moderate plasma β, the resulting convection is inward, meaning it acts to strengthen the rotation. However, for high-β and low aspect ratio, as found in spherical tokamak plasmas, sign reversals of the convection have been obtained within gyrokinetic simulations [442, 443].

Due to the Coriolis momentum pinch, the rotation in the plasma core is strongly dependent on the magnitude of the rotation at the plasma edge. This is true to such a degree that the edge rotation is expected to be critical for providing strong rotation in the core of ITER and other large tokamak devices. This is discussed more in section 5.3.1.

Experimental measurements, over a wide range of machine and plasma conditions, have demonstrated the presence of intrinsic rotation in fusion plasmas [406, 408, 413, 444, 445]. In particular, the observation of rotation profiles that exhibit a finite rotation gradient at zero rotation, even when considering shifts of the profiles of the order of a few km s⁻¹ to account for possible differences between main and impurity ion rotation, cannot be explained by any mechanism other than residual stress. This showed that significant rotation, and rotation shear, can be established in the core of tokamak plasmas even if the applied NBI torque and the Coriolis momentum flux are small. In consequence, the transport community and the ITPA TC group dedicated significant effort to understanding turbulence driven intrinsic rotation, which led to a rapid development of the corresponding momentum transport theory.

In this section we review pieces of the residual stress that are presently believed to be most relevant for the comparison with present day experiments and important for predicting rotation in future, larger machines. These include:

• residual stress arising from up–down asymmetry in the magnetic equilibrium. This is the only residual stress term at lowest order in $\rho_*$.
• residual stress arising from first order contributions in $\rho_*$ to the flow shear.
• residual stress arising from profile shearing/intensity gradient effects.

For a comprehensive overview of all the identified residual stress mechanisms, see recent theoretical reviews [417, 418].

A residual stress induced by an asymmetry of the flux surface itself ($C_\mathrm{FS}$, see equation (10)) was described in [446, 447] using linear simulations and with nonlinear simulations in [448]. This momentum flux was demonstrated experimentally in [449, 450]. The sign of the up–down asymmetry residual stress depends upon the signs of the toroidal magnetic field and toroidal plasma current, and on the direction of the up–down asymmetry of the flux surface. All of these dependences were validated experimentally [450]. The dependence upon the flux surface shape is subtle and involves competing contributions from the local magnetic shear and magnetic drifts [447]. The main effect can, nevertheless, be captured by considering the tilt of the flux surface main axis with respect to the vertical direction. If the flux surface is tilted in the direction opposite to the ion $\mathbf{B}\times\nabla B$ drift (for instance, a flux surface with a positive top triangularity and a negative bottom triangularity for a clockwise toroidal magnetic field as viewed from above), $C_\mathrm{FS}$ will tend to enhance the toroidal rotation peaking in the co-current direction. For ITER, and tokamak plasmas in general, this mechanism is expected to be significant only at the plasma edge, where the flux surface asymmetry is strongest. However, through the Coriolis momentum pinch, edge generated intrinsic rotation can be extremely important in determining the core rotation profile. Therefore, this generation mechanism, together with others that are expected to be edge localized [451, 452], may be important when combining edge and core physics in integrated models for rotation profile predictions.

A second contribution to the residual stress arises from the deviation of the plasma flow from a purely toroidal and species independent E × B flow [453, 454]. The relationship between the toroidal flow shear, $\frac{\partial \Omega_\phi}{\partial r}$, and the E × B flow shear γ_E is given by [119, 455]:

$$\begin{equation} \gamma_E = -s_j\frac{1}{|\nabla r|}\frac{RB_p}{B}\frac{\partial \Omega_\phi}{\partial r} + \gamma_E^*\end{equation} \tag{ 11 }$$

with

$$\begin{equation} \gamma_E^* = s_j\frac{1}{|\nabla r|}\frac{RB_p}{B}\frac{\partial}{\partial r}\left \langle \hat{u}_\theta \frac{B_t}{R}+\frac{1}{Z_sen_s}\frac{\partial r}{\partial \Psi}\frac{\partial p_s}{\partial r}\right \rangle .\end{equation} \tag{ 12 }$$

Here, $\hat{u}_{\theta} = (\mathbf{v}\cdot\nabla\theta) / (\mathbf{B}\cdot\nabla\theta)$ is the ratio of the poloidal rotation to the poloidal magnetic field and p_s is the species pressure. For the main ion species, neglecting radial derivatives related to the magnetic equilibrium, but retaining radial derivatives of the kinetic profiles, and using the neoclassical poloidal flow for $\hat{u}_\theta$, one obtains:

$$\begin{equation} \gamma_E^*\sim -\rho_*v_\textrm{thi}\frac{\partial}{\partial r}\left(k_iR/L_{Ti}\right) +\rho_*v_\textrm{thi} \frac{\partial}{\partial r}R/L_{pi}\end{equation} \tag{ 13 }$$

where $v_\textrm{thi} = \sqrt{T_\mathrm i/m_\mathrm i}$ is the ion thermal velocity, $k_\mathrm i$ the neoclassical flow coefficient ($k_\mathrm i\sim 1$ in the banana regime), $R/L_{T_\mathrm i}$ the normalized ITG and $R/L_{p_\mathrm i}$ the normalized ion pressure gradient.

The two terms contributing to $\gamma_E^*$ are due to the poloidal rotation and to the diamagnetic flow, respectively. They scale with $\rho_*$, but can become comparable, or even exceed, the toroidal rotation contribution to γ_E in the presence of steep profiles such as those encountered in the edge pedestal. They enter the parallel and E × B flow shear contributions to the momentum flux as:

$$\begin{align} \left(\frac{\Pi_\phi}{n_im_iR_0}\right)_{\parallel+E\times B} = -\left(\chi_\parallel\frac{B_\mathrm t}{B}-\chi_E\frac{B_\mathrm p}{B}\right)\frac{\partial \Omega_\phi}{\partial r} - s_j\chi_E\gamma_E^*.\end{align} \tag{ 14 }$$

As discussed in section 5.2.2, for positive magnetic shear, the diffusive flux driven by E × B flow shear is positive, $\chi_E \gt 0$. In this case, and with $\gamma_E^* \gt 0$, an inward flux of co-current momentum is generated that tends to enhance the peaking of the rotation profile in the co-current direction.

Note that, in the literature, E × B residual stress often refers to the momentum flux related to γ_E and, therefore, includes a part of the diffusive flux that is lowest order in $\rho_*$. This is the first term on the right hand side of equation (11). In the following, we will make the distinction between the total E × B momentum flux, related to γ_E and the E × B residual stress that only includes the component related to $\gamma_E^*$, i.e. the second term on the right hand side of equation (11), which is first order in $\rho_*$. The generation of a momentum flux by the E × B shear is closely linked to the tilting of the turbulent eddies away from the low-field side mid-plane (θ₀ = 0). The resulting stress is proportional to the tilt angle, θ₀, which is determined by the amplitude and sign of γ_E. In ballooning coordinates, the tilt of turbulent eddies translates to a shift in the radial wave-vector spectrum (k_x). In [307, 456] a model capturing the k_x spectrum shift induced by E × B shear was developed and proved to accurately capture the E × B momentum flux in quasilinear models (together with the turbulence stabilization by E × B shear).

Together with the E × B flow shear, other mechanisms related to the radial variation of background profiles, such as profile shearing [434, 457], including first and second derivatives contributions, and the turbulence intensity gradient [458–460] can generate a residual stress momentum flux. These mechanisms formally require a radially global treatment of the gyrokinetic equation and involve next order terms in $\rho_*$ that are neglected in the flux-tube approximation.

The impact of the radial shear in the turbulence propagation frequency, connected to the radial variation of the density and temperature across the simulation domain (first derivatives), has been explored in the linear regime in [457]. The eddy tilting generated by the mode frequency shear was shown to induce a residual stress proportional to the tilt angle θ₀. The magnitude of the residual stress is higher for ITG modes than for TEMs and changes sign at the ITG/TEM transition. This component of residual stress peaks the rotation profile in the co-current direction for ITG modes and in the counter-current direction for TEM modes.

The impact of second derivatives of the density and temperature profiles was addressed in nonlinear simulations in [434]. The resulting picture is very similar to that of the E × B residual stress and the resulting momentum flux can be formally written as

$$\begin{equation} \left(\frac{\Pi_\phi}{n_\mathrm im_\mathrm iR_0}\right)_{2^{nd} der.} = - s_j\chi_n\gamma_n - s_j\chi_T\gamma_T\end{equation} \tag{ 15 }$$

with $\gamma_n\sim\rho_*v_\textrm{thi}\frac{\partial}{\partial r}R/L_n$ and $\gamma_T\sim\rho_*v_\textrm{thi}\frac{\partial}{\partial r}R/L_T$ the second radial derivatives of the density and temperature profiles, respectively. For the reference case explored in [434], $\chi_n\sim\chi_E \gt 0$ and $\chi_T\sim-\chi_E$ is observed, meaning that positive γ_n, as typically encountered in the edge region, will drive an inward flux of co-current momentum whereas positive γ_T (also routinely observed in the edge) will drive an outward flux of positive momentum. For density driven TEM turbulence, $\chi_n\sim\chi_E \gt 0$ was still observed, but $\chi_T\sim0$. Similarly to the E × B residual stress, there is a saturation of the driven flux due to χ_n and χ_T decreasing towards zero at high values of γ_n and γ_T. The dependence of profile shearing residual stress (from first and second derivatives) with respect to the sign of the magnetic shear has not been explored to date. Due to the sign dependencies of the first and second derivatives and the relative importance of the gradients of the different background profiles, it is difficult to determine the net sign of the profile shearing contributions with a simple rule of thumb. Rather, global modeling is required to assess the total contribution.

Modeling momentum transport is extremely challenging due to the wide variety of contributing mechanisms whose relative importance, over a wide and multi-dimensional parameter space, is unclear from simple models. In addition, the momentum channel displays strong inter-dependences upon a variety of plasma parameters, including the equilibrium profiles and their gradients, requiring a self-consistent treatment of these coupled dependences. In practice, gaining quantitative traction requires a gyro-kinetic or gyro-fluid approach, with gyro-kinetics providing a more general and appropriate framework. It has been shown that gyrokinetics, as implemented in simulation codes, conserves toroidal momentum [401, 461–463], a prerequisite in the study of momentum transport. A wide variety of local and global gyrokinetic codes has been applied to the study of momentum transport. Local simulations, in addition to being easier computationally, enable specific effects to be independently tested, providing a direct comparison of their relative importance. This is more difficult to do in global simulations, which simultaneously incorporate all of the residual stress mechanisms described in section 5.2.4. However, by implementing diagnostics within the global simulations, it is possible to identify those residual stress mechanisms responsible for the observed momentum fluxes.

Independent of the type of modeling, experimental validation is very challenging as experiments have only succeeded in separating diffusive, convective, and residual stress momentum fluxes within a limited parameter space, often providing only sparse information against which models can be validated. One exception to this are experiments pioneered at DIII-D in which the externally applied torque is tailored in order to zero the rotation profile across a large portion of the minor radius of the machine. Here, the diffusive and convective fluxes are, by definition, zero and the intrinsic torque profile inferred from experiment provides a clear validation quantity [445, 464]. These experiments were modeled by GYRO simulations that produced similar levels of intrinsic torque as seen in the experiment through a combination of E × B shear and profile shearing contributions [434]. However, the uncertainties were large due to strong sensitivities to the profile gradients. An important result from this work was that the residual stress from profile shearing and E × B shear can be of similar magnitude. Therefore, depending on their respective signs, can either compete or complement one another, resulting in either a cancellation or enhancement of the net residual stress. In general, the strong sensitivities seen within such global simulations place very strong requirements on the quality of the experimental data needed for such validation exercises.

Another research avenue are experiments focusing on intrinsic rotation, i.e plasmas without externally applied torque. In these scenarios, in steady state, the momentum fluxes must balance (both in experiment and simulation) to produce the measured rotation gradients, providing a clear validation criteria. This can be seen by rearranging equation (9) as

$$\begin{equation} u^{^{\prime}} = -\frac{\Gamma_\mathrm{RS}}{\chi_{\phi}}-\frac{R_0V_{\phi}u}{\chi_{\phi}},\end{equation} \tag{ 16 }$$

where a normalization to the ion thermal velocity, $v_\mathrm{th,i} = \sqrt{2T_\mathrm i/m_\mathrm i}$, has been introduced and the equation written in terms of the normalized rotation gradient, $u^{^{\prime}} = \left(-R_0^2/v_\mathrm{th,i}\right) \mathrm d\Omega_{\phi}/\mathrm dr$, normalized rotation $u = \Omega_{\phi}R_0/V_\mathrm{th,i}$, and normalized residual stress, $\Gamma_\mathrm{RS} = \Pi_\mathrm{RS}/\left(n_\mathrm i m_\mathrm i v_\mathrm{th,i}\right)$. The residual stress and convective momentum fluxes from simulations (normalized by the momentum diffusivity) determine the normalized rotation gradient in the simulations. These can be directly compared to the u′ measured experimentally, which are often observed to be of order 1 [38, 234, 465].

Although local gyrokinetic models do not contain all of the physics needed to correctly describe global effects such as profile shearing, see section 5.2.4, they can still be useful in probing the residual stress mechanisms that are consistent with the experimental data and in demonstrating the degree to which the global effects are reflected in the local models. In calculations with local models, the tilt angle in the ballooning transform or, correspondingly, the radial wave number, is usually set to zero, because this usually corresponds to the condition with maximum growth rate. However, a finite tilt angle, θ_i, can be applied to the mode structure in the simulations, which generates a finite residual stress from the eigen-function asymmetry, see for example [457]. Ideally, a tilt angle would be taken from an experimental measurement or a global simulation, however, neither of these are regularly available. Within local, quasilinear simulations, the residual stress is linear in the applied tilt angle such that a prediction of the rotation gradient for an arbitrary value of θ₀ is possible as

$$\begin{equation} u^{^{\prime}}_p = -\frac{\Gamma_\mathrm{RS}}{\chi_{\phi}}\frac{\theta_0}{\theta_i}-\frac{R_0V_{\phi}u}{\chi_{\phi}}.\end{equation} \tag{ 17 }$$

Hence, when comparing quasilinear results to experiment, the value of θ₀ can be found that best reproduces the measured gradient. This approach was applied to databases of AUG intrinsic rotation observations in [38, 465] and showed that different tilt angles were required to describe the rotation profiles observed in plasmas with dominant ITG versus dominant TEM instabilities. Moreover, by fixing the applied tilt angle to that which best described a given dataset on average, it is possible to explore the dependences of the generated residual stress as a function of other input parameters. This procedure revealed a clear dependence of the produced residual stress on the normalized electron density gradient, $R/L_{ne}$, consistent with the experimental data, see section 5.3.4.2, and consistent with the results of the nonlinear global modeling, which later identified the density curvature as a critical parameter [466]. These works speculated that the experimental data can be qualitatively explained by a combination of residual stress from equilibrium E × B flow to explain hollow rotation profiles, and a mechanism dependent on the sign of the mode frequency, to explain co-$I_\mathrm p$ directed rotation profiles present in ITG dominated plasmas. Interestingly, these results are also qualitatively consistent with later global modeling of DIII-D intrinsic rotation profiles [467], which identified zonal E × B flow as the dominant driver of hollow rotation profiles, while the co-$I_\mathrm p$ directed intrinsic rotation profiles in TEM could be explained with turbulent diffusion, without the need to invoke a residual stress contribution. This is discussed further in section 5.3.4.2.

The large values of u′ observed in experiment are challenging for simulations, particularly within local models, as demonstrated in Hornsby et al [468]. Therein, a large database of intrinsic rotation observations from AUG L-mode plasmas [38] was modeled by nonlinear and quasilinear GKW [30] simulations. These included the Coriolis pinch [437], the residual stress from up–down asymmetries [446], and neoclassical equilibrium distribution function effects [417, 469], obtained by coupling GKW to the neoclassical transport code NEO [174, 176]. This tested whether the asymmetry in the gyrokinetic equation produced by neoclassical and diamagnetic flows [469, 470] could reproduce the AUG intrinsic rotation data behavior and, in particular, if it was large enough to drive the level of rotation gradient, u′, seen in experiments. Previous work using this model found success in predicting the sign of the intrinsic rotation in MAST [471], but neither its magnitude nor its gradient. The work by Hornsby et al [468] is important as it shows quantitatively that the momentum fluxes from all of the considered mechanisms in the local model are of similar magnitude and that the summed flux leads to significantly smaller gradients than observed in the experiment.

In contrast to local simulations, global, delta-f simulations of residual stress were able to obtain normalized rotation gradients of the same magnitude as experiments [466, 472]. In [472] global, nonlinear simulations with GTS are able to produce residual stress profiles, which, when integrated outside of GTS to obtain intrinsic rotation profiles, quantitatively reproduce the experimentally observed hollow rotation profile, see figure 21. Similar results were obtained in [466] using global, delta-f simulations with GKW. Here, the experimental hollow rotation profiles in AUG were quantitatively reproduced, self-consistently, within the simulation. It is interesting to note that these results were obtained with two very different, but complementary, approaches; GKW is a Eulerian Vlasov code while GTS is a particle in cell code. The GKW approach additionally highlighted the importance of including kinetic electrons, as the residual stress changes sign compared to using adiabatic electrons. This effect was first seen in full-f, gyrokinetic simulations using a hybrid kinetic electron model [473]. There, the sign of the intrinsic rotation gradient, and consequently the sign of the intrinsic rotation, changed sign when switching from adiabatic to kinetic electrons and sizable intrinsic flows up to 0.15 of $v_\mathrm{th,i}$ were obtained.

Zoom In Zoom Out Reset image size

The work with global GKW [466] also demonstrated a strong influence of the density profile curvature on the residual stress, which is sufficient on its own to explain the correlation between the core intrinsic rotation and electron density profiles observed in AUG experiments [38, 465]. In fact, using the ultra fast sweeping reflectometer [474], a linear correlation between the first and second derivatives of the density profile was established in the experiment [466], connecting the experimental and simulated observations. The correlation from the modeling can be seen in figure 22 and provides a clear demonstration of the presence of a density-curvature-driven, profile-shearing, momentum flux, as described in equation (15) in section 5.2.4.

Zoom In Zoom Out Reset image size

This work also points out, however, that the curvature of the density profile cannot be the only key parameter, as significant variation in the residual stress is also observed for the same density curvature, i.e. within the AUG database, profiles with differently directed rotation gradients can be identified for the same value of $R/L_{ne}$. This is very reminiscent of results from Alcator C-Mod, where strong changes in the intrinsic rotation have been observed under conditions of constant electron density and temperature profiles (and gradients) within the uncertainties of the measurement [475, 476]. Such observations are not in contradiction with the AUG results. Analysis of individual profiles is complementary to the database approach. Whereas large collections of observations over a wide parameter range enable broad correlations to be identified, individual profile comparisons demonstrate that additional drive terms are required to explain the observations. In addition, such examples present challenging targets for global modeling to understand why the rotation reversals occur in what appears to be a very narrow region of parameter space. Conversely, the sensitivities to first and second derivatives seen in the simulations pose a challenge for experimental validation. While measurements of kinetic profiles and their gradients have improved greatly in recent decades, reliable measurements of their curvatures remain, with few exceptions, beyond the capabilities of most diagnostic systems.

It is important to note that the simulations done with GKW in [466] did not distinguish the profile shearing residual stress produced by zonal E × B flow from the total profile shearing residual stress. Rather, both contributions are included in what is termed ‘profile shearing’ in this work. Hence, these results are likely consistent with other works that specifically identify the zonal E × B flow as a key component. For example, the GTS simulations of the intrinsic rotation in a DIII-D LOC-SOC transition [467] show that the zonal E × B flow shear correlates most strongly with the residual stress in the case of hollow rotation profiles and in [477], again modeling DIII-D experimental data with GTS, the zonal E × B flow and turbulence intensity gradient are identified as the dominant symmetry breaking mechanisms. In contrast, turbulence intensity gradient was not found to be important in [466].

Global, delta-f gyrokinetic codes are now able to reproduce residual stress momentum fluxes in the core of tokamak plasmas that correspond to intrinsic rotation profile gradients that are in quantitative agreement with experimental measurements. This represents a significant, and important, step-forward in our predictive modeling capability. However, due to strong sensitivities to small changes in the background plasma conditions, such successes do not yet translate to an easily applicable, predictive capability, particularly when extrapolating to the entirely new parameter spaces that will exist in future machines. Within global simulations, diagnostics must be implemented to identify sources of residual stress. Within the literature, different approaches to do this have been employed and, while no strong inconsistencies between the different works performed so far have been identified, the implementation of the same diagnostics within all codes would allow more stringent and conclusive assessments of the relative roles and dominance of the different mechanisms.

An important consequence of the $\rho_*$ ordering in the gyro-kinetic equation (see 5.2.1) is the predicted scaling of these momentum flux contributions to larger devices, such as ITER, which will operate at significantly lower $\rho_*$ than in present machines due to the larger machine size, and higher magnetic fields. It is possible that terms that are higher order in $\rho_*$, although potentially critical to explain observations in present day experiments, may not be dominant in affecting the rotation in profile in those conditions. Therefore, it becomes necessary to clarify the $\rho_*$ dependences. This can be approached experimentally, see section 5.3.4, as well as via gyro-kinetic simulations. This is critical in determining which sources of momentum will be negligible, or considerable, in future machines. This question was examined non-linearly with GYRO [434] and both linearly [457] and non-linearly with GKW [466, 478]. The GYRO simulations showed a linear scaling of the residual stress to diffusive flux ratio with the level of profile shearing, γ_n and γ_T, at low values of shearing and then a saturation when the level of shearing increases. Since γ_n and γ_T scale with $\rho_*$ this would suggest a residual stress to diffusive flux ratio proportional to $\rho_*$ at low $\rho_*$ values. The linear model produced a weaker scaling of $\rho_*^{1/3}$. In later nonlinear work with GKW in the adiabatic electrons approximation, a linear scaling of the residual stress to diffusive flux ratio with $\rho_*$ was obtained for low $\rho_*$ values [478], in agreement with the GYRO results [434]. However, subsequent nonlinear GKW simulations with kinetic electrons and based on AUG experiments revealed a lack of scaling with $\rho_*$ [466]. Together, these works suggest the $\rho_*$ scaling of residual stress to diffusive flux ratio ranges from $\rho_*^0$ to $\rho_*^1$ and may be parameter dependent. It is, therefore, not clear at present whether higher order $\rho_*$ contributions to the residual stress will be smaller or not in future devices. Smaller residual stress contributions in future machines, would result in the rotation profile being dominated by a combination of diffusive and convective transport and the edge boundary condition, which would become the key to accurate predictions of the rotation profile. Experimental attempts to illuminate the $\rho_*$ dependence of the residual stress are discussed further in section 5.3.4.

A great deal of progress has also been made in modeling plasma rotation via global full-f simulations. As already discussed, local delta-f models retain the lowest order terms in the $\rho_*$ expansion, and describe only turbulent transport induced by steady background profiles without finite $\rho_*$ effects. Global delta-f models also describe only turbulent transport under steady background profiles, but with finite $\rho_*$ effects. On the other hand, full-f simulations compute the gyrokinetic equation as it is without such a scale separation, and describe both neoclassical and turbulent transport and the resulting evolution of background profiles including zonal and radial electric fields. These features are important, in particular, for analyzing edge plasmas and transient plasma responses. The toroidal angular momentum conservation in full-f models is given as [401, 461], $\partial_t{\cal L}_\phi+\partial_\psi(\Pi_{\phi,\mathrm{turb}}+\Pi_{\phi,\mathrm{neo}})+\partial_\psi T_\phi = {\cal J}$, where ${\cal L}_\phi$ is the gyrocenter angular momentum, $\Pi_{\phi,\mathrm{turb}}$ and $\Pi_{\phi,\mathrm{neo}}$ are turbulent and neoclassical components of the Reynolds stress, T_φ is the field term stress, and ${\cal J}$ is the J × B torque due to the radial current. Here, the field term stress is determined by the phase shift between the perturbed density and fluctuating toroidal electric fields, and is also interpreted as the polarization stress [479] or the generalized Maxwell stress [461]. The toroidal angular momentum conservation was verified in several gyrokinetic simulations based on full-f models [461, 480–482]. The accuracy of momentum transport and radial electric fields was also verified by cross-code comparisons showing quantitative agreement in radial electric field profiles and intrinsic rotation profiles with different numerical approaches [279], and by computing the third-order gyrokinetics, which showed that higher order corrections on the quasi-neutrality condition are negligible for ion scale turbulence [481].

It is noted that Parra and Catto [483] raised a serious concern about the accuracy of momentum transport and radial electric fields in full-f models. They argued that in order to determine the correct transport of momentum, one needs to compute higher order momentum transport up to ${\cal O}(\rho_*^3)$ using third-order gyrokinetics, while most gyrokinetic codes are based on first-order gyrokinetics [484]. In this argument, it is assumed that the lowest order momentum transport is canceled by a microscopic symmetry of the local gyrokinetic equation in up–down symmetric tokamaks, and that the remaining momentum transport is ${\cal O}(\rho_*^3)$ under the low flow ordering, $v_\parallel/V_\mathrm{th}\sim{\cal O}(\rho_*)$, and assuming a gyro-Bohm scaling of the momentum diffusion. These assumptions are valid in the $\rho_*\to 0$ limit (local flux-tube model), but at moderate $\rho_*$, the momentum transport does not necessarily scale as ${\cal O}(\rho_*^3)$. For example, in a quasilinear theory [457], a weak $\rho_*$ scaling was shown for the profile shear stress $\Pi_*/\chi_\phi\propto \rho_*^{1/3}$. In fact, global delta-f simulations in circular equilibrium geometry for $1/\rho_*\unicode{x2A7D} 220$ showed that the profile shear stress scales as $\Pi_*/(n_\mathrm im_\mathrm iR_0 \chi_\phi)\sim v_\parallel^{^{\prime}} \propto \rho_*^{0}$, and thus is comparable to the lowest order momentum transport [466]. Full-f gyrokinetic simulations using the third-order gyrokinetics [485] showed that higher order corrections have a negligible impact on momentum transport and the resulting intrinsic rotation induced by ion scale turbulence in circular equilibrium geometry at $1/\rho_*\sim 150$ [481]. Here, the third order gyrokinetics [485] was derived by using the standard first order guiding-center transform in $\epsilon_B\sim \rho_\mathrm i/L_B\sim \rho_\mathrm i/R$ and the third order gyro-center transform in $\epsilon\sim \rho_\mathrm i/L_n\sim \rho_\mathrm i/a$, where $\rho_\mathrm i$ is the ion gyro-radius, L_n is the scale length of the density gradient, and L_B is the scale length of the magnetic field inhomogeneity. This formulation may be valid for large aspect ratio tokamaks, in which ${\cal O}(\epsilon_B^2)\sim{\cal O}(\epsilon^2\epsilon_A^2)$ correction terms can be neglected under the inverse aspect ratio expansion in $\epsilon_A = a/R\ll 1$. However, depending on the choice of ordering, other higher order terms can be derived [486]. Their impact on momentum transport in gyrokinetic simulations remains to be examined. Another critical issue is $\rho_*$ scaling of various symmetry breaking mechanisms. If the lowest order momentum transport is cancelled in future large machines with small $\rho_*$, higher order momentum transport can be important. This issue remains an area of ongoing research.

Fundamental physical mechanisms of momentum transport have been investigated with full-f models, and the following properties were disclosed. Firstly, avalanche like non-local heat transport due to interactions between temperature profiles and turbulent fluctuations is universally observed in full-f simulations, and similar non-local transport features were observed also for momentum transport [487, 488]. Secondly, in addition to the turbulent Reynolds stress, which was discussed based on delta-f models in the previous sections (see equation (9)), other players contribute to the toroidal angular momentum conservation in full-f models. In ITG turbulence with adiabatic electrons, turbulent and neoclassical contributions of the Reynolds stress are dominant and tend to cancel with each other [461, 480]. This feature was further examined by full-f simulations over an energy confinement time, where the toroidal angular momentum conservation in the steady state showed that the turbulent Reynolds stress due to avalanche like non-local transport cannot cancel only by itself, but the remaining turbulent Reynolds stress is cancelled by the neoclassical Reynolds stress, which is induced by the magnetic drift, and is largely enhanced by non-axisymmetric turbulent fluctuations [489]. The field term stress also makes a finite contribution, and is enhanced by kinetic electrons, which modify the phase relation between the perturbed density and fluctuating toroidal electric fields, leading to the sign flip of intrinsic rotation [473]. In steady plasmas, the J × B torque is negligible in the plasma core. However, in the presence of transient particle transport, e.g. in modulation heating experiments, the J × B torque due to the ion radial current can also make a significant contribution [490]. Thirdly, most full-f simulations showed the development of intrinsic rotation from zero initial rotation without external momentum input, and the turbulent Reynolds stress was shown to correlate with the E × B shear [255, 488, 491] and the turbulent intensity gradient [488]. The profile shear stress theory [457] was also examined in full-f simulations, and it was shown that the turbulent Reynolds stress correlates with the tilt angle θ₀, which depends on the profile shear and the E × B flow [489]. According to the higher order ballooning theory [492], θ₀ is determined by the sum of the profile shear and the E × B shear. In addition, global simulations show that θ₀ is affected not only by the E × B shear, but also by the magnitude of E × B flow itself [489, 493]. Although the profile shear stress theory [457] was proposed without taking account of the E × B flow, it enables a unified treatment of these effects through the tilt angle θ₀.

Numerical experiments on the intrinsic rotation and the residual stress driven by ITG turbulence with adiabatic electrons were conducted with GYSELA [494], XGC1 [495], and GT5D [255]. In [461, 480], $\rho_*$ scan numerical experiments were conducted with GYSELA at far above marginal $R/L_{Ti}\sim 10$ for $\rho_* = 1/512\sim 1/64$, where the divergence of the Reynolds stress scaled as $\partial_{\psi}\Pi_\phi/(Rv_\mathrm{th,i}^2/a)\propto\rho^{*0.7}$. On the other hand, in [491], $\rho_*$ scan numerical experiments were conducted with GT5D at near marginal $R/L_{Ti}\sim 6$ for $\rho_* = 1/225\sim 1/100$, and the $\rho_*$ dependence of the Mach number of intrinsic rotation was shown to be weak, while in [496], this $\rho_*$ scan was extended up to $\rho_* = 1/450$, which showed lower Mach number at smaller $\rho_*$. In [491], by using an adaptive momentum source, which imposes a nearly zero rotation condition, the residual stress was shown to correlate with the E × B shear. By comparing simulations with and without an adaptive momentum source, it was shown that the intrinsic rotation enhances the E × B shear and the resulting turbulence suppression. In [488], numerical experiments with XGC1 using the kinetic profiles from DIII-D experiments [497] showed the development of intrinsic rotation in the co-current direction, and a nearly zero rotation condition was produced by imposing an external torque in the counter-current direction, which confirmed the intrinsic torque proportional to the normalized temperature gradient $R/L_{Ti}$. This relation is consistent with the correlation between the intrinsic rotation and the temperature gradient in the experiment [498, 499]. In [487], an impact of scrape-off-layer flows on the core intrinsic rotation was investigated with GYSELA, and it was shown that large-scale, avalanche-like events transport both heat and momentum, and scrape-off layer flows can modify the core intrinsic rotation profile in a large fraction of the radial domain.

In [500], the inward momentum transport phenomenon at the edge region of an H-mode plasma on DIII-D [501] was analyzed with XGC1, and it was shown that the edge rotation inside the separatrix is generated mostly from neoclassical physics such as the X-point orbit loss and the finite radial orbit effect [502], where turbulent fluctuations enhance the ion orbit loss by modifying the velocity space distribution near the loss cone. In [503], the development of intrinsic rotation in the presence of magnetic ripple was investigated with GYSELA simulations and a synergistic impact of the magnetic ripple on the turbulent Reynolds stress was identified, by which the ripple indirectly effects the Reynolds stress through an enhancement of the E × B shear. An impact of TEM turbulence on the intrinsic rotation was also investigated with full-f simulations including kinetic electrons. In [490], an electron heating modulation numerical experiment was conducted with GT5D, and it was shown that transitions from ITG turbulence to TEM turbulence induced by electron heating generate density peaking and rotation changes in the counter-current direction. Here, toroidal torque flipping between ITG and TEM phases is related to reversal of the ion radial current or the J × B torque, indicating the coupling of particle and momentum transport channels in transient plasma responses.

In summary, significant progress has been made on modeling and understanding momentum transport in tokamak plasmas using full-f simulations. This progress is predicated upon the demonstration that toroidal angular momentum is conserved in such simulations. While full-f and delta-f models share common physics regarding the turbulent Reynolds stress, full-f models demonstrate the existence and importance of additional players, such as the neoclassical Reynolds stress, the polarization stress, and the J × B torque. Full-f simulations are particularly important to understand the intrinsic rotation generation in the plasma edge, where they are better suited to the task compared to delta-f models. Edge plasmas are characterized by large amplitude fluctuations and strongly distorted velocity distribution functions (e.g. via orbit losses due to magnetic ripple and proximity to the open field-lines), which are not captured in delta-f simulations. However, delta-f models may still be valid, provided the target plasma is well inside the separatrix. Full-f models have made good progress in modeling edge intrinsic rotation and rotation changes in transient plasma responses. However, modeling of the plasma boundary and including all required physics in the edge region remain significant challenges. Moreover, due to the complexity of the simulations and required computation resources, they are not well suited to fast, integrated modeling packages aimed at predicting behavior in future machines. Rather, they provide a complementary approach, aimed at identifying the most important physical processes and providing a guide and validation point for simpler model development.

In order to study momentum transport, a variety of experimental techniques have been developed. The best and cleanest way of separating diffusive and non-diffusive transport is by means of perturbation experiments, i.e. exploiting the additional information contained in the dynamic response of the plasma angular momentum to a time variation of the applied torque. Momentum perturbation experiments were carried out at several tokamaks using a variety of techniques to modulate the applied torque including NBI modulation, NBI blips, magnetic perturbations, and changes in plasma confinement. Exploiting the technique most appropriate to the machine, experiments were performed in JT-60U [507], JET [508–510], AUG [234, 511], DIII-D [512, 513], NSTX [514], MAST-U [515], KSTAR [516] and Alcator C-Mod [517, 518].

These early experiments were the first pioneering efforts to experimentally separate these fluxes. Due to the complexity of the problem, the analyses mostly employed simplified momentum transport models. The most common assumptions made included constant Prandtl numbers, temporally constant convection profiles, and two term momentum transport equations including only diffusive and convective components, i.e. neglecting any residual stress. Under this last assumption, an inability to match the experimental data provides indirect evidence of residual stress, for example, see [507, 519].

More recent works have improved the rigor of the momentum modulation analyses, pushing towards more quantitative theory validation through more complete models and demonstrating the importance of maintaining as much complexity as can be supported by the data. Several recent works have highlighted the importance of maintaining the time dependence of the transport coefficients as well as the impact of the ion heat channel on the momentum transport [511, 516, 520, 521].

Despite the relatively simple models employed, early torque modulation experiments were able to successfully demonstrate that the propagation of momentum in fusion plasmas cannot be modeled using diffusive terms alone, providing more quantitative evidence on the size and direction of non-diffusive fluxes. These works also demonstrated changes in the core momentum transport with changes in local plasma parameters, providing trends for comparison with theoretical predictions. As discussed in section 5.2.2, the Prandtl number is predicted to be of order one, but otherwise does not show any strong dependences that can be cleanly tested within the experiment. Nonetheless, this is, in itself, a prediction that can be examined and the experimental work exploring the behavior of the Prandtl number in tokamak plasmas is presented in section 5.3.2.

In contrast to the limited variation expected of the Pr number, theory provides clear predictions for how the Coriolis momentum pinch should scale, see section 5.2.3, together with a wide variety of predictions for the residual stress, see section 5.2.4. The experimental trends for non-diffusive fluxes observed in perturbation experiments are discussed in section 5.3.3. Few experiments have successfully separated the diffusive, convective, and residual stress components within a given experiment. Rather, experiments aiming to study the residual stress tend to focus on scenarios where the momentum and convective fluxes are expected to be negligibly small. These are presented in section 5.3.4.

Results from experiments at JET [522], which explored the momentum transport as a function of q, $R/L_{n_\mathrm e}$, and $\nu^*$, show good consistency with the theoretical predictions. They yielded Pr numbers around 2 and could not identify any correlation between Pr and the examined parameters. Dedicated gyrokinetic simulations of these plasmas using GS2 resulted in a similar lack of correlation with the scanned parameters and predicted Pr numbers lower than those determined experimentally. This under prediction may be due to the absence of residual stress in the analysis, or to other simplifications such as time-independent transport coefficients, and hence do not necessarily indicate a discrepancy between theory and experiment. This speculation is supported by recent results from AUG [511, 521]. In this work, NBI modulation experiments are analyzed retaining the diffusive, convective, and residual stress momentum fluxes as well as their time dependences and, thereby, demonstrates, for the first time, quantitative agreement between the measured and predicted momentum transport coefficients. In [511] the momentum transport in very similar hydrogen and deuterium plasmas is investigated, and found to be in good agreement, consistent with the theoretical expectation of only a very small isotope dependence of the transport coefficients in the plasma core. This can be seen in figure 23, where the D results from this comparison are shown together with the GKW predictions for the Prandtl and pinch numbers for both deuterium and hydrogen. The agreement between experiment and theory is within experimental uncertainties across the entire plasma minor radius, and the predicted difference between H and D is significantly smaller than the experimental error bars.

Zoom In Zoom Out Reset image size

The parameter dependences of the Prandtl number were also explored in JT-60U plasmas [519], where a variation of Pr between 1 and 3 was found, with Pr increasing with increasing β_N. Similar to the JET results and theoretical expectations, no dependences with $\rho_*$ or $\nu^*$ were observed. High Prandtl numbers have also been observed in spherical tokamaks, when the ion heat conductivity approaches neoclassical levels [515]. However, in the outer-half of NSTX Pr of order unity was obtained from modulation experiments [512, 514], similar to values found in conventional aspect ratio tokamaks. For ITER, and other large, low collisionality devices, it is expected that turbulent transport will dominate over collisional effects. As such both the ion heat and momentum transport are expected to remain anomalous in these devices.

The momentum diffusivity in DIII-D and NSTX was found to increase with increasing plasma collisionality [464], as did the convective velocity as will be discussed in the next section. The behavior of Pr was not reported, and it is worth noting that other plasma parameters were not held constant over these scans. The dependence of the momentum confinement time τ_φ on $\nu_*$ and $\rho_*$ was studied later in DIII-D. The momentum confinement time is not a measure of momentum diffusion, but rather a global measure of the net momentum transport, where higher confinement times can be the result of lower diffusion, stronger inward convection, intrinsic torques reinforcing the rotation direction, or combinations thereof. In these experiments τ_φ increased with decreasing $\rho_*$ and $\nu_*$, corresponding to a favorable scaling of the momentum confinement for ITER [523]. This work also reproduced the earlier result of lower absolute values of $\chi^\mathrm{inc}_\phi$ at lower $\nu_*$ and extended this to lower q. Both of these trends are favorable for ITER rotation [513].

The most recent in-house ITER toroidal rotation predictions are presented in [524]. The predictions indicate an ITER rotation of the order of 100 km s⁻¹. However, the calculation depends strongly on the assumption for the momentum confinement time τ_φ. The modeling approach in the ASTRA simulations uses accurate NBI torque calculation and assumes that the momentum transport can be described by a momentum confinement time τ_φ, with values as found in [504]. While this approach can, in principle, model the momentum diffusion appropriately, the predictions for the rotation may be incorrect in plasmas where the momentum pinch is significant. In particular, neglecting the momentum pinch may lead to the conclusion that the central rotation only has a weak dependence on the pedestal rotation. In contrast, earlier modeling studies for ITER plasmas that accounted for the momentum pinch concluded a strong dependence on the edge rotation values [522]. As a consequence, the ITER predictions presented there depend strongly on the edge intrinsic torque, which drives edge intrinsic rotation, and influences the core rotation via the pinch.

Finally, a multi-machine, momentum-transport-coefficient database was created within the ITPA TC group [525]. This database combined data from Alcator C-Mod, DIII-D, JET, NSTX and JT-60U and provides a means of examining the parameter dependences of the transport coefficients over a wider range than is accessible in any single device. This database demonstrated the strong experimental correlation between the ion heat conductivity and the momentum diffusivity, which results in limited variation in the Pr number, despite a two order of magnitude changes in the absolute value of the transport coefficients. By combining measurements from multiple machines, the database also notes a trend towards larger Pr numbers at lower collisionality and higher Pr values obtained in bigger machines (JET and JT-60U) compared to the smaller ones (C-Mod, NSTX, and DIII-D). Detailed gyrokinetic modeling of these datasets is required to understand these observations and remains pending.

When studying momentum transport with a two term (diffusion and convection) momentum transport model, it is tempting to compare the experimental convective term to predictions for the Coriolis pinch. However, the possibility of non-negligible residual stress should be retained, which re-casts such a two term description into a separation of diffusive and non-diffusive fluxes, opposed to diffusive and convective fluxes. Inherent to this statement is the assumption that the residual stress will tend to mimic pinch-like behavior in a perturbation analysis, but not alter the diffusive flux, which is not, a priori, proven. In fact, the residual stress is expected to scale with the momentum diffusivity as well as local plasma parameters, both of which may vary in time during momentum modulation experiments, depending on how the torque perturbations are engendered. This is particularly relevant for NBI modulation experiments at low frequency, which are often accompanied by clear modulation of the temperature profiles. Such changes can result in time dependent momentum transport coefficients including the residual stress, which would impact all aspects of the perturbation analysis. The impact of such perturbations on the analysis and derived transport coefficients are discussed, for example, in [511, 516, 520] together with methods to reduce the uncertainty in the deduced values.

However, for many experiments, the experimental steady state angular momentum profiles and the momentum propagation can be sufficiently well described by the two term (diffusive plus convective) equation, leaving little room to identify residual stress within the uncertainties. This ansatz, has been applied to modulation data in multiple fusion devices including JT-60U [507], JET [508–510, 526], AUG [234, 511, 520], DIII-D [464, 513], NSTX [512, 514], MAST-U [515] and KSTAR [516]. An important result from these experiments is that a finite inward convection was always required to best fit the experimental measurements of the momentum propagation, demonstrating the existence of convective transport and explaining why an effective Prandtl number is found that is significantly lower than 1, see section 5.3.1.

In multiple machines the inward pinch is observed to increase with increasing χ_φ [464, 519, 525], resulting in only moderate changes in the pinch number $-RV_{\phi}/\chi_{\phi}$. This correlation between the diagonal and off-diagonal components is expected from theory [437], making the pinch number a more relevant physics quantity for understanding momentum transport than the convective velocity itself. Moreover, this correlation has direct analogues with other transport channels, for example, impurity transport [160], which also demonstrates a clear correlation between the diagonal and off-diagonal elements 3.2.4. In both cases, when turbulent transport is strong, it is strong in both the diffusive and convective fluxes.

The parametric dependences of the momentum pinch were studied in JET. It was found that the magnitude of the inward momentum pinch depends strongly on the inverse density gradient length, with an experimental scaling for the pinch number of $-RV_{\phi}/\chi_{\phi} = 1.2R/L_n + 1.4$ [522]. This was also qualitatively reproduced over a much wider range in $R/L_{ne}$ in the multi-machine database [525], see figure 24. Here, two values of q₉₅ were isolated and for both data-sets the pinch number was observed to increase with increasing $R/L_{ne}$. For the JET dataset [522], gyro-kinetic simulations agreed qualitatively with experimental results, but the dependence on $R/L_{ne}$ was much weaker in the simulations and the absolute size of the pinch smaller. A similar result was obtained at AUG [520], where the measured pinch number was larger than predicted by gyrokinetic calculations, although with a similar profile shape. On DIII-D, the inferred inward pinch velocity from a rotation braking experiment [464] showed quantitative agreement to theoretical pinch predictions [437, 438] using simple fluid models that neglect the parallel dynamics and, as a result, maximize the predicted pinch velocity, see section 5.2.3. Gyrokinetic calculations that account for kinetic electrons would likely reduce the size of the predicted pinch, leading to the same under prediction of the pinch number seen in JET and AUG. In subsequent works from AUG [511, 521], the discrepancy in the pinch number between theory and experiment from the original publication [520] was shown to be the result of not including the time dependent residual stress self-consistently within the modeling. When all terms were included, quantitative agreement with the quasilinear gyrokinetic simulations for both the Prandtl and pinch numbers was achieved, as can be seen in figure 23.

Zoom In Zoom Out Reset image size

In the JET modulation experiments, no dependence of the pinch number on $\nu^*$ was observed, whereas the pinch depended weakly on the q-profile, with the pinch number decreasing with increasing q, qualitatively consistent with the predictions in Peeters et al [437]. The dependence of the pinch number on $\nu^*$ was also explored in DIII-D and NSTX plasmas where only a weak dependence was obtained from the data, consistent within uncertainties with the JET result [497, 523]. At JET, momentum transport was also studied via a database approach [506], which yielded clear evidence for non-diffusive transport. The parameter dependences of the non-diffusive portion were explored via regression techniques and compared to gyrokinetic predictions of the Coriolis pinch. For both the experimental data and the Coriolis predictions, the strongest dependences were on $R/L_{ne}$, q, and the inverse aspect ratio, ε, which is used as a proxy for the trapped particle fraction. Similar regression parameters were found for both data-sets. These results support the idea that the Coriolis pinch is responsible for much of the non-diffusive transport in those JET plasmas, but also left enough margin for a non-negligible residual stress momentum flux, particularly in the outer half of the plasma.

Since the theoretical identification of the Coriolis momentum pinch as a non-negligible momentum flux with the potential to provide significant rotation in future devices, a wide variety of experiments and analyses have been carried out with the aim of validating this prediction. These experiments successfully demonstrated the existence of an inward convective momentum flux of comparable magnitude to the diffusive flux and have demonstrated some basic parameter scalings that are qualitatively consistent with those predicted by theory. The measured non-diffusive fluxes are of similar magnitude to predictions of the pinch, but are often observed to be larger. Recent work suggests this discrepancy could be resolved by including time dependent residual stress in the analysis of the modulation data, indicating that the agreement between theory and experiment, at least in ITG plasmas, is potentially better than it already appears. That a strong inward convection of momentum is predicted by theory and seen in experiment has important consequences for modeling the rotation profile of future machines. In particular, predicting the value of the rotation at the plasma edge, or pedestal top, becomes critical to determining the direction and magnitude of the core rotation profile, as discussed in [522]. A recent review of experimental and theoretical work on the edge intrinsic rotation generation can be found in [418].

Evidence for the role of $\Pi_{\mathrm{RS}}$ in driving intrinsic rotation in fusion plasmas has been steadily building. Examination of the velocity profile evolution following changes from co- to counter-current NBI in JFT-2M plasmas demonstrated that there is a finite momentum flux even at zero rotation gradient and a finite offset velocity gradient even at zero momentum flux [406]. Further evidence came via a completely different approach, namely, observations of intrinsic rotation in C-Mod plasmas [407, 444], which demonstrated an edge localized intrinsic momentum source [527]. Subsequently, intrinsic rotation measurements were also made in many tokamak devices including TCV [528], DIII-D [529], C-Mod [530], JT-60U [531], JET [532], NSTX [533], AUG [38], MAST [471], KSTAR [534], and EAST [535]. The presence of co-$I_\mathrm p$ directed core-localized residual stress in NBI heated plasmas was demonstrated in novel torque balancing experiments at DIII-D [464]. The same technique was later applied in TEXTOR [536] and EAST [535] to demonstrate the presence of co-$I_\mathrm p$ intrinsic rotation in those devices and quantify their magnitudes. A different approach was taken in the J-TEXT tokamak [537], where a biased electrode was used to influence the edge rotation. With a negative bias, the naturally occurring co-$I_\mathrm p$ intrinsic rotation could be balanced, enabling the identification and quantification of a co-$I_\mathrm p$ intrinsic torque density in the edge of those plasmas. In contrast, in the presence of strong electron heating, counter-$I_\mathrm p$ directed, core localized, residual stress could be inferred from rotation measurements at DIII-D [538] and later at AUG [140].

These observations, together with observations of inward convection (section 5.3.3), demonstrated the existence of momentum transport mechanisms capable of competing with turbulent diffusion. As such, they became of interest to the ITPA TC topical group as potential mechanisms by which rotation would be generated in future devices. With this in mind, experiments were carried out to study the parameter dependences of the residual stress, identify dominant generation mechanisms, and explore their extrapolation to future machines.

5.3.4.1. Edge localized residual stress.

The generation of residual stress at the plasma edge is particularly interesting for future machines for two reasons: the $\rho_*$ scaling of the residual stress and the interaction with the Coriolis momentum pinch. Since residual stress terms are predicted to scale with $\rho_*$ (see section 5.2.4), they are expected to be smaller in future devices than in present day machines, see section 5.2.1. However, this scaling comes from an expansion of the gyrokinetic equation in $\rho_*$ $\sim \sqrt{T_\mathrm i}/B L_{\perp}$. In the core the plasma minor radius, a, is a good approximation for $L_{\perp}$. However, in the plasma edge this is not the case. Rather, when considering pedestal generated residual stress, it is the local perpendicular scale length that should be taken, i.e. $\rho_{*, \mathrm{edge}} \sim \sqrt{T_i}/\left(B L_{\perp,\mathrm{edge}}\right)$. However, without knowledge of the pedestal profiles in future machines, how this quantity will extrapolate is uncertain. If pedestal perpendicular scale-lengths in future machines are similar to those in present day devices, there is hope that pedestal generated residual stress at higher order in $\rho_{*, \mathrm{edge}}$ might be a relevant player. This could then provide sufficient edge rotation for the Coriolis pinch to convect inward and form peaked rotation profiles, even in the absence of significant applied torque.

The study of edge generated intrinsic rotation is complicated not just by the uncertainty on which parameters should be considered as important ($\rho_*, \rho_{*, \mathrm{edge}}, v_{E\times B}, \ldots$), but also by the presence of non-turbulence related rotation generation mechanisms. These include neutral drag [400], NTV torques [398, 399] (see also 6), ion orbit loss [397], and magnetic ripple [539, 540], which can strongly impact the plasma rotation via the loss of NBI fast ions as well as by increased ion orbit losses and ripple-associated NTV torques. Moreover, in the presence of RMPs, as foreseen for the baseline operating scenario in ITER, additional physics mechanisms can come into play that can drive toroidal plasma rotation [541, 542], (see also section 6).

Neutral drag is expected to be relevant in the SOL, potentially contributing to a boundary condition, but it is not expected to be important in the confinement plasma region. NTV torques and the impact of RMPs on plasma rotation are certainly important within the confined plasma and have been observed to cause strong changes to the plasma rotation, see for example, [541–543]. These effects are increasingly well understood, and in many cases, their predicted impact on the rotation has been compared to experiment and found to be in good agreement. Presently, a comprehensive model, including all of these sources of intrinsic rotation in the plasma edge, remains an outstanding challenge. Ultimately, simple models that reflect the physics of these effects must be built and incorporated in integrated modeling frameworks in order to provide a physics-based prediction of the rotation in the edge of future devices.

Most intrinsic rotation experiments are conducted under conditions where magnetic perturbations, and their associated effects, are not expected to be responsible for the observed changes. Ion orbit loss, on the other hand, is more difficult to separate from turbulence induced torques. Not only do they exist in roughly the same radial region and contribute on the same order of magnitude, they are also expected to depend on related quantities. Ion orbit loss depends on the velocity distribution of the edge ions, and hence on the square root of the ion temperature, and it is also sensitive to the radial profile of the edge ion temperature. Similarly, residual stress generation mechanisms are expected to be related to gradients of the edge kinetic profiles (or their curvature) and the shear in the radial electric field, which is driven by the main ITG.

The co-existence of multiple generation mechanisms necessarily muddles the interpretation of intrinsic rotation measurements in this region and their scaling with local parameters. Nonetheless, the scaling of the intrinsic rotation in the edge has been investigated in multiple devices with the joint aims of identifying driving mechanisms and providing empirical extrapolation of the intrinsic rotation to future machines. The first attempt to do this was an inter-machine comparison of intrinsic rotation [408] including data from six devices: JET, Alcator C-Mod, Tore Supra, DIII-D, JT-60U, and TCV. This work examined the change in rotation between L-mode and H-mode. Strictly speaking, this work did not focus on the edge rotation. Rather, the data included was measured in most cases near the plasma center, with the exception of the DIII-D data, which was in the edge on account of the counter-current core rotation present in the core due to the use of ECRH, see section 5.3.4.2. Therefore, the influence of momentum transport on the derived scalings cannot be ruled out. However, if the changes in the core rotation are dominated by changes in the edge intrinsic torque associated with the formation of an H-mode pedestal, empirical scalings of the core data would provide indirect information on the edge intrinsic torque. This work provided the first empirical scaling of the observed intrinsic rotation.

Several key observations resulted from this inter-machine comparison. First, no dependence of the intrinsic rotation is seen as a function of $\rho_*$  over a factor of three variation in $\rho_*$ within the data-set. Rather, each machine is defined by its respective value of $\rho_*$ and, at that value, is able to produce wide variations in intrinsic rotation values, which show significant overlap with those obtained in the other devices. This can be seen in figure 25. This demonstrates that the intrinsic rotation is strongly sensitive to other parameters and any hope of discerning a $\rho_*$ dependence has to simultaneously hold these unknown variables constant and employ a $\rho_*$ variation significantly larger than that in the matched parameters and their uncertainties. As the data from the different machines were not obtained under dimensionless and gradient-matched conditions, the lack of a clear scaling with $\rho_*$ is unsurprising and does not suggest that particular types of turbulence are or are not involved in the intrinsic rotation generation process.

Zoom In Zoom Out Reset image size

A second observation from this inter-machine comparison is that the change in the rotation in dimensional units scales almost linearly with the change in the plasma pressure (ΔP), the magnetic field ($B_\mathrm t$), and inversely with the plasma current ($I_\mathrm p$). This is similar to the observation made on several devices that intrinsic rotation scales with the plasma stored energy [544–546], or stored energy divided by the plasma current, $W_\mathrm{th}/I_\mathrm p$ [529, 547], an observation commonly referred to as the Rice scaling. In [408] significant rotation in ITER is predicted, of the order of 300 km s⁻¹. However, all scalings since then, which typically take into account additional physics considerations, result in significantly lower predicted rotation values. A work by deGrassie et al [548], finds that the intrinsic rotation data at the edge of DIII-D plasmas, but far inside the pedestal top, is well described by a dimensionless scaling of the form $M_\mathrm A \sim \beta_N \rho_*$, where $M_\mathrm A$ is the Alfven Mach number. This scaling also describes a set of C-Mod and JET data well and yields more modest predictions for the edge rotation in ITER of only 26 km s⁻¹. This scaling is also qualitatively consistent with a linear scaling found from nonlinear gyrokinetic studies, as discussed in section 5.2.5.

The scaling of the C-Mod core intrinsic rotation data was later revisited, connecting the change in rotation between L-mode and advanced confinement regimes to local, rather than global quantities [499]. In that work, the change in rotation is shown to correlate with changes in the edge temperature gradient and not with the edge density nor the pressure gradient, similar to NSTX [533], suggesting that the ITG in the pedestal, or a related quantity such as E × B shear, could be a dominant driver of pedestal generated residual stress. Due to the radial electric field well present in H-mode plasmas [549, 550], the sign of the E × B shear changes sign near the edge of the plasma (at the minimum of the well). This means the sign of the intrinsic rotation driven by this mechanism will also reverse at this location. To be consistent with the experimental observations, the sign of the intrinsic rotation in the inner-shear layer should be co-$I_\mathrm p$ directed. In the pedestal region, the contribution of the toroidal flow to the E × B shear (γ_E) is small and γ_E is effectively given by the diamagnetic pressure and poloidal flow, $\gamma_E\sim\gamma_E^*$, see equation (12). In the outer E_r shear layer, close to the separatrix, $\gamma_E^* \lt 0$ and the corresponding momentum flux drives counter-$I_\mathrm p$ directed rotation. In contrast, in the inner E_r shear layer, the E × B residual stress drives co-current rotation, consistent with the co-$I_\mathrm p$ directed intrinsic rotation seen in the experiment and the correlation with the edge ITG.

Other works have also identified correlations between intrinsic rotation generation and the ITG, for example, in internal transport barriers in LHD [498]. Similar observations were also made in [551] where the C-Mod data compare well against an E × B shear based residual stress model [451] and yield moderate core rotation (Mach numbers of $M_\mathrm i$ = 0.07, $V_{\phi}\sim 50$–75 km s⁻¹) predictions for ITER and much stronger rotation in ARC ($M_\mathrm i$ = 0.2, $V_{\phi}\sim 200$ km s⁻¹). These predictions are significantly lower than the first scalings based on global quantities in [408], but still twice as big as predicted by deGrassie et al [548]. It should also be noted that the residual stress model used in that analysis has been criticized for the use of inconsistent ordering, as discussed in [452]. In addition, while the model works well for C-Mod, it has not found success at ordering the edge intrinsic rotation data at DIII-D [552].

Empirical scalings provide estimates for the order of magnitude of rotation expected in a future machine and can capture the dominant dependences, which can help to guide theory. However, without a physics-based understanding of the mechanisms generating this intrinsic rotation, it is not possible to assess the reliability of such empirical scalings. Work in this direction was performed at DIII-D where the intrinsic torque profile was measured by zeroing the rotation profile using externally applied torque [464]. Using this technique a database of edge localized intrinsic torque was developed and was shown to correlate well with both the global parameter $W_\mathrm{th}/I_\mathrm p$ and the local edge pressure gradient $\nabla P$. The data are qualitatively consistent with the idea of pedestal residual stress generation by strong pedestal gradients. If this is the case, then the original scaling from [408] is optimistic, as the pedestal gradients in future machines will be limited by peeling-ballooning stability. The work on DIII-D was later expanded to cover a wider parameter space in [553] where it was observed that only considering a dependence on the local pressure gradient was insufficient to organize the experimental data well. However, when an ion orbit loss model [410] was taken into account, the intrinsic torque could be well described, demonstrating the importance of multi-mechanism interplay to accurately describe the rotation generation processes.

A third effect present in the pedestal is the impact of the magnetic configuration on the intrinsic rotation [446]. As discussed in section 5.2.4, there is a residual stress momentum flux associated with up–down asymmetry in the magnetic configuration. This was predicted theoretically in [446, 447] and demonstrated experimentally at TCV [449, 450], where it was shown that a poloidal tilt of the magnetic flux surfaces influences the local gradient of the intrinsic rotation. The effect is strongest at the edge, where the flux surfaces are most strongly asymmetric, but due to the momentum pinch, the entire rotation profile is altered by the change in the edge rotation gradient. The effect was shown to be quite large, and a change in the sign of the up–down asymmetry directly translated to a change in sign of the rotation gradient. Hence, the shape of the magnetic flux surfaces can be used as a tool to tailor the rotation profile, provided up–down asymmetric flux surfaces are present and compatible with future reactor operation. Moreover, this effect will increase the net pedestal generated intrinsic torque if the surfaces are tilted the right way and negate a portion of this torque if tilted the wrong way. This effect should be kept in mind when comparing the pedestal generated intrinsic torque between devices as well as between different scenarios in individual machines. In the presence of strong edge gradients, as present in H-mode plasmas, this contribution would be difficult to disentangle from gradient driven stresses.

This edge intrinsic torque model developed at DIII-D [553] was also applied to NSTX data and to DIII-D advanced inductive scenarios, which feature significantly higher intrinsic torque than observed in the original data set, and found to work well for both [554]. Interestingly, if this model is extrapolated to ITER conditions, the predicted edge intrinsic torque is only 2 Nm, which is smaller than the applied neutral beam torque in that region, and would be insufficient to drive significant rotation. This type of analysis requires the ability to zero the rotation profile through the application of counter-$I_\mathrm p$ directed external torque, which is only possible in a limited number of devices. Experiments of this type were recently performed at EAST [535], where they demonstrated the presence of co-$I_\mathrm p$ directed intrinsic torque, peaked at the plasma edge, consistent with the DIII-D result. More such experiments at DIII-D, EAST, and other larger devices with counter-$I_\mathrm p$ NBI would certainly be of interest for further exploration of the parameter dependences of the edge intrinsic torque, particularly with respect to local dimensionless parameters. In a recent work from AUG [511], the residual stress profile was extracted from NBI modulation experiments concomitantly with the diffusive and convective fluxes, and also demonstrates the presence of a co-$I_\mathrm p$ directed intrinsic torque, peaked at the plasma edge, consistent with the data from balanced beam experiments. The advantage of such a technique, is that all three components of the flux are obtained for the same plasmas, providing comprehensive theory validation possibilities. Moreover, this technique can also be applied to devices without balanced beams and over a wider range plasma parameters.

A similar experimental method was applied in a recent work at DIII-D that aimed to study the parameter dependences of the intrinsic torque as a function of $\nu_*$ and $\rho_*$ [523]. The impact of these parameters on the intrinsic torque is assessed through comparison of plasma discharges with otherwise matched dimensionless parameters, where an order of magnitude variation in $\nu_*$ and a factor of 1.4 in $\rho_*$ was achieved. This work shows a decrease in the intrinsic torque at lower $\nu_*$ and a surprising increase for lower $\rho_*$. This result holds for both the absolute torque measured and when normalized by the ion temperature, although it is not clear that $T_\mathrm i$ is an appropriate normalization when extrapolating to other machines. In general, normalization is required, as the extrapolated quantity should be dimensionless for a physically meaningful scaling. Moreover, a proper normalization should account for changes as a function of dimensional parameters that may also change when matching dimensionless parameters between machines.

In [523], the net extrapolation to ITER parameters taking into account the scaling of the momentum confinement time is positive, predicting 33 Nm of torque. This level of intrinsic torque is comparable to the projected NBI torque and may be enough to influence the ITER rotation profile, although estimations of the predicted ITER rotation including this torque yield quite low absolute rotation levels compared to present day devices [513]. The scaling of the intrinsic torque with $\rho_*$ was also explored in C-Mod plasmas and compared to the intrinsic torque observed in JET and DIII-D in [517]. In this work, the plasma discharges were selected by attempting to match the plasma shape, β, and q₉₅ and resulted in a similar level of total intrinsic torque for ITER (∼20 Nm) and lower levels predicted for SPARC and ARC (∼4 and 8 Nm), although proper normalization when extrapolating to future machines remains an open question.

There is an important caveat to these types of analyses. Many residual stress terms are predicted to scale with the normalized gradients of kinetic profiles or even their curvature. This leads to extreme sensitivities of the intrinsic torque to changes in the profiles that appear small when comparing absolute values. As such, despite displaying good matches in many dimensionless parameters, it is very difficult to assess if the first and second derivatives of the kinetic profiles are sufficiently matched to better than the variation achieved in $\rho_*$. A significantly larger variation in $\rho_*$ than can realistically be achieved in a single device is required to unambiguously investigate the $\rho_*$ dependence of the residual stress. In [517] almost a factor of 3 in $\rho_*$ was achieved, by inter-machine comparison, but with the inclusion of an additional complication. Namely, the choice of normalization is critical, this is discussed further in [517, 523]. In these multi-machine studies, matches in dimensionless parameters was required, but it is not clear what the normalization for the intrinsic torque should be. In [523] the ion temperature was used, while in [517], normalization of the intrinsic torque with R (rotation frequency) and R² (moment of inertia) were considered.

Residual stress in the edge pedestal is a potential mechanism by which sufficient toroidal rotation can be generated in a future machine. Many devices have demonstrated co-$I_\mathrm p$ directed intrinsic rotation associated with H-mode pedestal formation, and in particular, with the temperature pedestal. These studies are suggestive of the dominant role of E × B driven residual stress in this region or, alternatively, any residual stress model that is proportional to the edge temperature gradient, for example, the model presented in [555–557], which combines ion-orbit-loss physics with local turbulence gradients in the plasma edge, and has found success in describing TCV [557, 558] and AUG [559] data. The extrapolation of pedestal generated intrinsic rotation and intrinsic torque to future machines has been examined by a variety of experimental approaches. Most recent scalings agree that the intrinsic torque provided by such mechanisms will be modest, but could be sufficient to influence the pedestal top rotation, and hence, the core rotation via the convective flux. However, without a first principles understanding of the rotation generation mechanisms, significant uncertainties remain on such predictions. Further experimental work assessing the intrinsic torque generated in the edge pedestal as well as gyro-kinetic modeling of the residual stress in the pedestal could both help to reduce the uncertainties on such extrapolations. In addition, integrated models including not only turbulence generated residual stress, but also effects such as ion-orbit loss, ripple, and NTV would be required to provide a more comprehensive picture and assess the relative importance of the difference mechanisms.

5.3.4.2. Core localized residual stress.

This section of the Special Issue project documents the tremendous progress that has been made in understanding the physics of the influence of 3D fields on transport in tokamak and stellarator plasmas. The mechanism of magnetic field penetration into fusion plasmas and how it leads to either suppression of ELMs or even their mitigation is not discussed here, as it is very extensively reviewed in the Pedestal Chapter of this Special Issue [21]. Here, the focus is on the effects of the 3D fields on the various transport mechanisms in tokamak and stellarator plasmas. The scope of topics in this chapter covers experimental observations and theory/modeling/simulations in an integrated way for neoclassical and turbulent transport. Additional focus is provided for the LH- transition and the pedestal region, where there is an overlap between both transport channels and additional edge effects and MHD will influence the observations and underlying physics.

The chapter is organized as follows: after introducing the basics of various 3D topology concepts, the historical (mostly experimental) work that found the earliest hints of changes in transport due to external fields is presented. This is followed by section 6.2 on neoclassical transport, where, after the basics of the neoclassical picture, an influence on plasma torque and rotation is discussed. Section 6.3 discusses the effects of 3D fields on turbulent transport in both: tokamaks and stellarators. Additionally, in section 6.4 influence of the 3D field on the L-H transition in both tokamaks and stellarators is discussed. The final section (section 6.5) reviews some aspects of the transport changes at the pedestal due to magnetic perturbations.

There are many variations of how 3D fields are present in the magnetic fusion devices and how they influence the plasma transport. For instance, 3D fields can introduce a simple striation of the magnetic flux surface (like a ripple or the core of a stellarator) without breaking flux surfaces. This can have a mirror-like effect and alter the trapping/passing ratio of particles [575]. Altering the 2D magnetic field configuration [576–578] of a tokamak as well as the intrinsic 3D nature of stellarators [579] and RFPs [580] allows for multiple different magnetic geometries. Stronger 3D effects can create breaking of magnetic flux surfaces and form island structures [581], with a clear X- and O-point. This 3D magnetic geometry has a strong impact on measurements as the transport across the island is diametrically different as compared to the transport across the closed magnetic flux surfaces (see section 6.3.3.2). These islands or island chains form on rational surfaces and if their width is large enough and/or the rational surfaces are close to each other, the field lines in the volume encompassed by the island become stochastic [582]. It should be noted, however, that [582] discusses a case of fully developed ergodic volume, which is often not present in fusion plasmas due to, for example, screening of the external fields by plasmas [583]. A stochastic boundary brings an additional challenge as this breaks the simple picture in a 2D configuration [584], there is a clean split between field lines that do not interact with the vessel wall and field lines that do (Scrape-Off Layer). The creation of stochastic magnetic fields close to this boundary can lead to a region of field lines where so-called closed and open field lines intermingle.

Both tokamaks and stellarators have utilized 3D fields to alter the magnetic field geometry of the LCFS or separatrix. An interesting use of 3D fields was a concept of so-called ‘ergodic divertors’ developed at Tore Supra [576] and TEXTOR [578, 585]. The applied RMPs lead to such a strong corrugation of edge magnetic field lines that it formed an exhaust concept. Similarly, this exhaust concept was first employed on W7-AS [579] and more recently on W7-X [586]. The island chain can be altered using external 3D fields to alter the overall magnetic equilibrium, which directly affects the mixing of open and closed fields lines. Similar to the ‘ergodic divertor’ experiments on limiter tokamaks, this region of mixed field lines can be split into 2 zones [584, 587, 588], based on the length of the open field line with respect to the Kolmogorov length (which is a measure of the e-folding length of the exponentially growing separation of neighboring trajectories) [589]. Field lines shorter than the Kolomogorov length form the so-called laminar region, while those that are longer are in the ergodic region [577]. Both regions exhibit different transport and turbulence characteristics, which indicates that understanding the impact of 3D fields on TC requires a good understanding of the underlying magnetic field structure. This is often a challenge experimentally and errors in understanding the structure can lead to results that appear in contradiction.

While stellarators are designed to have 3D magnetic fields, ideal tokamak are axisymmetric and are often thought of that way. However, due to a finite number of coils that generate the toroidal B-field, a magnetic ripple is introduced, which can have a strong impact on confinement [539, 590]. Additionally, due to misalignment of the various coils, error-fields can be introduced, which can be corrected by using external and internal coils that produce a small radial magnetic field. These same coils can also be used to produce RMPs as well as Non-Resonant Magnetic Perturbations (NRMPs). The difference between a resonant and a NRMP is in its alignment to the pitch of the magnetic field lines. Resonant (contrary to non-resonant) magnetic perturbations are aligned with the pitch of the field lines on a certain resonant surface, i.e. $m/n = 3/1$ component of the RMP is aligned with the field lines forming the q = 3 flux surface. RMPs are designed to be predominantly resonant at the plasma edge and are regularly employed to suppress or mitigate ELMs. ITER will have several sources of 3D fields in addition to RMPs. Even with ferritic inserts, the toroidal ripple will be fairly large on ITER (∼0.5%), and test blanket modules containing ferritic steel will produce additional error fields [591].

The earliest reports on the effects of 3D fields on transport in tokamaks were published by Evans et al [592], who discussed the effect of an externally applied resonant magnetic field perturbation on TEXT plasmas. Multiple island resonances from a so-called resonant helical divertor were used to study the effects of a stochastic boundary and magnetic islands on local recycling, particle removal, and plasma density. Similar results were obtained at COMPASS, where it was shown that particle confinement could be either improved or degraded depending on the magnetic perturbation spectrum [593].

Early successful applications of the 3D fields lead to installation of the ‘ergodic divertor’ first at the Tore Supra [577] and later at the TEXTOR tokamak [585]. Both experiments reported density changes, in particular experiments at TEXTOR showed that particle confinement is attributed to a modification of the E × B shear rate near the resonance surface [594] due to open field lines in the stochastic boundary connecting deeper plasma regions with the plasma-facing components [595]. Measurements of density fluctuations with the reciprocating probes at TEXTOR showed that in the ergodic zone, the local turbulent particle flux reverses direction from radially outwards to inwards, which causes a steepening of the local pressure gradient [596]. The turbulence characteristics were strongly affected by penetrating 3D fields: energy redistribution occurred in frequency spectra and large scale turbulent eddies were suppressed, and the radial and poloidal correlation lengths decreased. Mordijck et al [597] analyzed L-mode DIII-D plasmas and found that the plasma transitions from an ITG to TEM regime at the plasma edge. Moreover, the application of RMPs to L-mode plasmas affects the power threshold of the L-H transition. These results indicated that the Rechester—Rosenbluth model [582] is too simplistic to cover the whole complexity of the transport changes at the plasma boundary due to external magnetic perturbations.

Although the first results on the suppression of Type I ELMs with RMPs at DIII-D indicated that there were not much changes in the H-mode pedestal parameters [598], the later analysis by Moyer et al [599] showed that the plasma electron and ion temperature and the plasma electron density at the pedestal top are reduced during the ELM suppression phases. An important finding of this work was that application of the RMPs leads to an increase in small amplitude, radially narrow events with increased magnetic and density fluctuations. Moreover, changes in density fluctuations resulting from the RMPs correlate with the density inverse scale length [600]. McKee et al measured with the Beam Emission Spectroscopy (BES) diagnostic a strong increase of the long-wavelength turbulence in the outer regions of DIII-D low-collisionality plasmas, when RMPs were used to suppress Type-I ELMs. Experiments with modulated currents of the perturbation coils suggested that turbulence is directly affected by the RMPs [601].

The profound influence of 3D fields on neoclassical transport has been well acknowledge in the stellarator community [602], which lead to optimization of neoclassical transport in Wendelstein 7-X [603]. In the case of tokamaks, the first observation of the NTV torque braking by magnetic perturbation was made experimentally [543] at DIII-D. The reduced NTV model [604] showed that E × B shear rates and poloidal and toroidal flows can be predicted by neoclassical theory.

The influence of 3D fields in both tokamaks and stellarators is a very complex phenomenon. In particular, in tokamak H-mode plasmas the modification of transport processes due penetration of the magnetic perturbation, results from the complex interplay between MHD physics, turbulence characteristics, and neoclassical transport. Both turbulent and neoclassical transport are strongly connected to radial electric field and the E × B shear and, thus, may influence each other.

A plasma with net toroidal, radially differential, rotation would be desirable for micro- and macro-instability suppression. As described in section 6.3, the rotation impacts turbulence through $E_\mathrm r$ shear. It is thus important to understand effects that can suppress any toroidal plasma rotation. Non-ambipolar neoclassical fluxes arise in the presence of 3D fields, yielding a radial current. This current drives a toroidal torque, known as NTV. NTV causes some damping of toroidal velocity if the velocity is larger in magnitude than the offset velocity. Under typical conditions, the ambipolar $E_\mathrm r$ is negative [570], called the ion root regime in the stellarator literature [607]. This results in an ion flow in the counter-current direction that is diamagnetic in size [608].

RMP will be applied using in-vessel coils in ITER to control ELMs [609]. While RMP fields are intended to resonate with the magnetic field lines, indicating that the poloidal, m, and toroidal, n mode numbers satisfy $m - n q = 0$ for the safety factor q on some rational surface, they also provide non-resonant perturbations. This will likely impact the rotation profile, since the NTV torque is largely driven by non-resonant perturbations. While the resonant electromagnetic torque remains localized near rational surfaces, see section 6.2.2, NTV torque tends to have a broader radial profile. In addition to RMPs, other sources of perturbations include toroidal ripple from the finite number of coils [539] and test blanket modules [610]. These 3D fields are further modified by the plasma response [611], as described further in section 6.2.2.

As the NTV arises as a kinetic effect, the scaling and form of the torque density depends on the scale of the effective collision frequency, ν, in comparison with the bounce frequency, $\Omega_{\mathrm{b}}$, E × B toroidal precession frequency, $\Omega_E$, and magnetic drift toroidal precession frequency, $\Omega_{\mathrm{m}}$. The early theory of Shaing [612, 613] describes the low-collisionality regimes when $\nu \lt \Omega_E$ and $\nu \gt \Omega_E$, where the viscosity scales as $\sqrt{\nu}-\nu$ and $1/\nu$, respectively. When $\Omega_E \gt \nu$, precession can effectively phase mix a perturbation along a particle trajectory, resulting in reduced transport. There remains a small fraction of particles satisfying a resonance condition between the bounce or transit frequency and the drift frequency,

$$\begin{equation} l \Omega_{\mathrm{b}} + n \left(\Omega_E + \Omega_{\mathrm{m}}\right) = 0,\end{equation} \tag{ 19 }$$

and experiencing elevated levels of transport. The resulting transport is said to be in the superbanana (SB) or superbanana plateau (SB-P) regime, depending on whether the resonance is broadened due to collisions or the orbit width, respectively. As described in section 6.2.3, these transport regimes arising due to 3D perturbations lead to enhanced particle transport in addition to rotation damping. Typically several of these transport regimes are present in different parts of velocity and physical space in a single discharge [614], and theory has been developed to connect the different regimes [615]. NTV theory has been well-examined with numerical modeling in different regimes [616–618].

Rotation braking due to NTV has been observed on many tokamaks including MAST [605], JET [619–621], DIII-D [604, 622, 623], NSTX [541, 624], KSTAR [625], TEXTOR [621], EAST [626], JT-60U [627], and Alcator C-Mod [628]. Rotation damping is also predicted for quasi-symmetric stellarators with the addition of symmetry-breaking perturbations [629, 630]. This effect has been observed during biased electrode experiments on the quasi-helical HSX stellarator [631, 632].

Rotation damping in the presence of 3D perturbation was first observed in DIII-D [633], where m = 1, n = 3 perturbations were applied to suppress NTMs and resulted in a reduction of the momentum confinement time by 83%. Experiments were then performed on JET with an m = 2, n = 1 perturbation [619], which demonstrated qualitative agreement with collisionless NTV theory in simplified cylindrical geometry. Experiments with n = 2 perturbations on NSTX have demonstrated a $T_i^{5/2}$ scaling of the torque [624], consistent with non-resonant, low collisionality transport [612]. Quantitative experimental validation of theory has improved with more sophisticated modeling efforts. Reasonable comparison between JET results and a bounce-averaged model including all collisionless regimes relies on accounting for both the electron and ion torques in addition to the plasma response. Modeling of a DIII-D QH mode discharge with an n = 3 perturbation [634] using a guiding center delta-f model [635] recovers the predicted total torque when the trapped and passing particles are accounted for in addition to the plasma response.

While NTV often damps rotation, it can also increase the rotation magnitude in low-torque regimes. Experiments on DIII-D with non-resonant n = 3 perturbations and varying NBI direction and magnitude verify an acceleration to the offset rotation in the counter-current direction [543, 622, 634]. The $T \sim V - V_{\mathrm{offset}}$ scaling of the torque and $1/\nu$ theory shows qualitative agreement with the rotation evolution [636]. This increase in rotation leads to an increase in energy confinement, see figure 26, likely due to an increase in the E × B shearing rate, see section 6.3.2, and enabling access to QH mode.

Zoom In Zoom Out Reset image size

The radial profile of the torque depends on the collisionality regime and $E_\mathrm r$ profile. The torque typically has a peak where $E_\mathrm r$ passes through zero, corresponding to enhanced $1/\nu$ transport. This can occur, for example, in the presence of a transport barrier (section 6.3.2), since $E_\mathrm r$ is typically positive in the core and negative in the transport barrier [637]. As the bounce-resonance condition depends on radius as well as rotation velocity, bounce-resonant transport also exhibits peaking near surfaces where the resonant condition (19) is satisfied, often near rational surfaces and in the core since the rotation must be sufficiently large. Studies on EAST with application of n = 1 perturbations indicate a double-peaking of torque density in the core and edge, consistent with bounce-resonant and $1/\nu$ transport, respectively [626]. The NTV torque density on JET with n = 1 perturbations is peaked in the core, which differs from NSTX and DIII-D experiments with higher-n perturbations, showing peaking at the edge [620].

It is desirable to maximize pedestal flow damping to allow maximal penetration in the edge, see section 6.2.2, while minimizing strong core flow damping, which can cause mode locking. While the edge and core rotation are often correlated, it may be possible to decouple them by reducing use of midplane RMP coils for ITER scenarios [638]. As RMP fields modify the ambipolar $E_\mathrm r$ in the pedestal, this can impact the E × B shearing rate, and thus microturbulence, see section 6.3.2. Gyrokinetic simulations of DIII-D plasmas indicate that the application of RMPs leads to a rapid change in $E_\mathrm r$ due to this neoclassical effect, leading to an increase in turbulent transport in the pedestal top [639].

There are opportunities to improve confinement through modification of the torque profile with carefully tuned 3D fields, analogous to the improvement of confinement through stellarator optimization. Rotation measured in the presence of an n = 3 EFC coil on NSTX allowed the identification of an optimal amplitude and phase to correct the intrinsic error field, which matches with the optimal parameters to reduce NTV through modeling [640]. Correction of this n = 3 error field has led to improved performance with sustained high rotation and an absence of core MHD activity [641]. NTV torque due to n = 3 perturbations and NBI have been used as actuators to control rotation on NSTX, successfully achieving the target rotation in much less than the momentum confinement time [642]. Non-axisymmetric perturbations designed to minimize the NTV torque density have recently been demonstrated to reduce flow damping and density pumpout on DIII-D and KSTAR [643].

As some ITER scenarios will operate above the no-wall stability limit, it is important to understand the sources and sinks for rotation for possible stabilization of RWMs [621]. Due to its relatively large moment of inertia and beam energy, the momentum provided will be relatively small in ITER. Thus the resulting rotation may be close to the neoclassical offset. A sufficiently large offset rotation may, therefore, be beneficial for confinement and stability [622].

The NTV torque from ripple is predicted to be smaller than the NBI and intrinsic torque drives on ITER [575]. However, the test blanket modules may further damp rotation, but it may be possible to correct with n = 1 fields [610]. The electron-driven NTV torque is expected to become important at the low collisionality expected for ITER, and consequently the offset rotation may be in the co-current direction [644]. As discussed in section 5.2.4, edge-generated intrinsic rotation, predicted to be in the co-current direction, is important for determining the core rotation due to the Coriolis momentum pinch. Even if the NTV torque is in the co-current direction, if it is significantly smaller in magnitude than other momentum sources, the resulting rotation may be reduced. Including NTV in integrated modeling is critical for determining the self-consistent rotation profile.

In a plasma described by ideal MHD, shielding currents arise on rational surfaces in response to RMPs in order to preserve the magnetic topology. With a finite amount of resistivity, reconnection can ensue, allowing magnetic islands to open at the rational surface. In a plasma, in which the rotation frequency is comparable to the resistive frequency, the change in topology is suppressed [645, 646]. Reduction of the resonant field by one to two orders of magnitude is observed when a resistive response is taken into account, for large enough plasma flows [647]. On the other hand, non-resonant components of the perturbation are not screened, can penetrate, and thus contribute to NTV, see section 6.2.1.

Shielding currents at rational surfaces induce a braking torque, known as an electromagnetic (EM) torque. This arises due to J caused by eddy currents near the rational surface and the magnetic field perturbation due to an RMP, providing a $\textbf{J} \times \textbf{B}$ torque [583]. The steady-state rotation frequency is determined by the balance between the viscous and EM braking forces. In contrast with the NTV torque, which typically scales as $T \sim V - V_{\mathrm{offset}}$, EM torque decreases with increasing rotation. Furthermore, EM torque is typically localized near rational surfaces, but can sometimes become more global due to radial diffusion [647].

The basic physical description of the interaction between 3D rotating magnetic structures and external magnetic perturbations (or resistive walls) was put forth by Fitzpatrick [645]. Here it was understood that rotating magnetic islands are stabilized by rotation, to some extent, if they are small enough. However, when the island width exceeds a critical value, mode locking occurs prior to complete stabilization. Because the EM torque achieves a local minimum for a given rotation frequency, below a critical rotation frequency, a bifurcation to very small rotation occurs [645, 646], known as mode locking. Slow rotation causes the penetration of the RMP and rapid growth of an island on a time scale associated with the resistive decay of current sheets. Due to the nature of EM forces, mode unlocking due to increasing of rotation occurs at significantly larger critical rotation value than that of mode locking, indicative of hysteresis.

The original picture has been further refined by improving on the plasma model by considering viscous-resistive effects, thus identifying three main parameters that characterize the response to a static, resonant, error-field of a rotating tokamak plasma in the vicinity of the rational surface. These are indeed the plasma viscosity, rotation, and resistivity. In total, eleven distinguishable response regimes have been identified [646]. An improved version of the original model of Fitzpatrick [646] was then proposed to include diamagnetic and semicollisional effects [648] to the layer response, that is finite equilibrium pressure gradients and electron thermal conduction. From the point of view of the layer response, electron diamagnetism can completely suppress a magnetic perturbation [648] for large enough values of the control parameter $\hat \beta = 0.5\beta_e L_s^2/L_n^2,$ via a screening of the perturbed current [649] (here β_e is the ratio of the kinetic to the magnetic plasma pressure, for the electrons, and L_s and L_n are some characteristic lengths for global shear and equilibrium density gradients, respectively).

The theory of resonant field screening and rotation is further modified in the presence of stochastic layers [650, 651]. A parallel current is established by electron streaming along field lines, while a neoclassical ion radial current arises to maintain ambipolarity. The associated self-consistent $E_\mathrm r$ modifies the screening of resonant fields, and the ion radial current leads to a modification to the rotation profile.

The inclusion of NTV (section 6.2.1) in the rotation balance can modify the locked mode physics. In a regime where NTV torque exceeds perpendicular momentum diffusion, the rotation is maintained at $V_{\mathrm{offset}}$ outside the resonant layer, and the threshold for locking typically increases with this effect [652–655]. Newer analytic models which include drift-MHD [656] and nonlinear [655] effects more accurately capture the linear scaling of the critical perturbation amplitude with density observed in JET, DIII-D, and Alcator C-Mod. Extrapolation of this theory to ITER indicates that avoidance of locked modes lies just within the capabilities of the error-field correction system [628].

Some aspects of the theory of the plasma response to magnetic perturbations at resonant locations have been investigated numerically [657], confirming the importance of semicollisional effects, but also pointing out the crucial role played by kinetic Alfvèn waves [649]. Other crucial phenomena, such as drift wave radiation [658], have been proposed to have an impact on the penetration of magnetic perturbations (either rotating or locked). A comprehensive theory [581] is yet not available, but high fidelity modeling [659] provides a general physical understanding in which screening effects increase for lower resistivity, stronger rotation, and smaller RMP amplitude. A further important effect is that of the nonlinear evolution of the radial electric field produced by RMPs. The E × B rotation due to such fields can be locally compensated by electron diamagnetic rotation. Then, RMPs can indeed penetrate and form magnetic islands. This is generally true for a large class of examples based on DIII-D, JET, and ITER type plasmas [659].

More direct experimental evidence of the interplay between magnetic perturbations, plasma rotation, and high fidelity layer modeling has been supported by extended MHD computations, which have been successful in reproducing the parametric space of ELM suppression via RMPs in the DIIID and KSTAR tokamaks [660]. Careful modeling and dedicated experiments brought the important conclusion that ELM suppression by RMPs can occur over wide q₉₅ ranges in DIII-D ITER similar shape (ISS) plasmas with n = 3 RMPs [661, 662]. Because the field penetration threshold decreases with density, the range of q₉₅ for ELM suppression increases with decreasing electron density. The access to this regime at substantially higher pedestal pressures is facilitated by applying higher toroidal mode number (n = 4) RMPs.

In conventional (far from quasisymmetric) stellarators, the NTV torque plays an outsized role in establishing the flow profile. The presence of $1/\nu$ transport in a boundary layer outside an island allows the EM torque to be overcome, and healing of the island [663], analogous to mode unlocking in a tokamak. The theory has been validated on TJ-II and LHD, which have observed the healing of magnetic islands with an increase in poloidal flow [664]. If the NTV is in the $1/\nu$ regime, it is observed to compete with viscosity and dominate the rotation profile for ITER parameters [636]. This leads to penetration in the edge and shielding in the core. The NTV torque has been computed via MHD-kinetic hybrid modeling of an ITER baseline scenario coupled with a drift-kinetic model. The results indicate that the EM torque dominates over the NTV torque for ITER parameters with n = 3 perturbations [665].

In practice, both the resonant and non-resonant torques play a role in determining the rotation profile in tokamaks with 3D perturbations. Plasmas respond to RMP fields by amplifying non-resonant harmonics that have the same sign of the pitch as resonant harmonics [647, 666] through a near-marginal kink-like perturbation. This leads to an effective increase in NTV that becomes more substantial near marginality. The effect has been observed in DIII-D [667], where the magnetic braking of an n = 1 perturbation was measured, showing, however, the importance of both resonant and non-resonant torques in determining the global plasma rotation profiles. Similar effects have been observed on KSTAR [668]. For non-resonant configurations, a “quiescent” behavior is observed in which no density pump out and loss of stored energy are detected. Resonant configurations feature stronger rotation damping along with confinement degradation and density pump out, in line with DIII-D results. The modeling of such phenomenon in ASDEX Upgrade (different n and phasing) [669] is consistent with the general picture: the resonant EM torque, the Reynolds stress, and NTV all contribute, to different extents, to the flow damping. The relative importance of each factor depends on the specifics of the configuration. Specifically, the NTV and the EM torque tend to be largest in the pedestal, for smaller Reynolds numbers.

The physics of impurity transport, neoclassical and turbulent, has been recently reviewed by Angioni [113] and is also discussed in section 3 of this chapter. One of the most important developments for neoclassical impurity transport studies is in the realization that poloidal asymmetries of the impurity density can develop [185, 199, 200], thus affecting their radial transport [188, 201, 202]. While this might be suggestive of specific features of nonaxisymemtric devices, where the electrostatic potential also shows variations on flux surfaces [670–672],it cannot be considered a truly 3D effect, but rather a collisional and finite poloidal Larmor radius effect [113]. Truly 3D effects are induced by the presence of (saturated) MHD modes, which can be mimicked with stellarator equilibrium codes [673]. However, contradictory results in JET hybrid scenarios [673] (where impurities peaking is observed and simulated) and ASDEX Upgrade ECRH centrally heated ones [253] (where an impurity hole is observed) are limiting our ability to make conclusive statements. Nevertheless, this area of research shows a great potential for synergistic tokamak-stellarator studies. Indeed, the observation of an impurity hole in a plasma with an ion internal transport barrier on the Large Helical Device [674] has been a key result, and a long standing open problem [675, 676] in the stellarator community. The analogy seems to apply also to the ITER pedestal, for which numerical studies have been performed to quantify the effects of non-axisymmetric fields on the ambipolar electric field and radial particle fluxes [677]. Here the authors studied both the toroidal field ripple and the externally applied fields for ELM control. We know that, in ITER, the toroidal field ripple effects on high-Z impurity transport are negligible [677]. However, ELM control fields cause a significant modification of the pedestal ambipolar electric field (becoming less negative), allowing for the appearance of multi-valued solutions for the electric field, analogous to core stellarator transport [607]. Thus, ELM control external fields significantly increase the outward neoclassical pedestal transport in ITER, for both the main plasma ions, and high-Z impurities.

Non-axisymmetry has an important effect on collisional particle transport for both impurities and the bulk species. For this reason, the formal tools used to address the physics of impurities and their 3D properties are the traditional ones of neoclassical theory [113, 678]. Since even in axisymmetry, magnetic field inhomogeneities due to the poloidal variation of the magnetic field strength cause the presence of trapped particles and, therefore, of neoclassical effects, it is not so surprising that the neoclassical theories of particle and impurity transport in tokamaks (whether 3D or not) and stellerators (whether modular or helical) are mutually interrelated. A ubiquitous particle transport effect induced by 3D fields is the large density pump-out phenomena routinely observed in ELM control experiments using RMP fields. It is now understood that this behavior is triggered by the same physics that produces the toroidal NTV torque, mediated by the drift resonance between the 3D field perturbation and the toroidal precession of trapped thermal particles [679]. Pump out becomes stronger at low collisionality and near the offset of the plasma rotation, which are relevant regimes for ITER. ELM suppression mechanisms reduce the E × B drift frequency at the pedestal, thus allowing for the resonance to occur. High-fidelity modeling of neoclassical transport in the presence of broken magnetic surfaces with the NEO and M3D-C1 codes has reproduced the density pump-out observed on RMP-suppressed DIII-D plasmas [680]. To date, the modeled pump-out level does not always match the experimental value. One of the reasons could be the effect of island stochasticity, which also affect the values of $E_\mathrm r$ [639]. Experimentally, a correlation between density pump-out, the increase of $E_\mathrm r$ with RMPs, and a spin-up in the co-current direction, possibly due to loss of electrons in the stochastic layer, is observed [681].

Several experiments have reported that after the application of RMPs, the changes that have been observed at the edge of both L- and H-mode plasmas can be reproduced using quasi-linear gyrokinetic modeling, without including any 3D deformations to the magnetic geometry. DIII-D H-mode plasmas showed that the density fluctuations in the steep part of the pedestal were directly correlated with the linear growth rates calculated by TGLF [600] and more recently similar results were obtained for the changes in electron temperature fluctuations [707]. Similarly, simulations with GTC showed that the changes in 3D field structure in similar H-mode plasmas on DIII-D were not sufficient to directly influence changes in turbulence characteristics [708].

In DIII-D L-mode plasmas, the applications of RMPs did lead to small changes in gradients at the plasma edge. Using a perturbative gas puff to probe particle transport, a reversal of the particle pinch from inward to outward was found, and this switch in direction coincided with a change from ITG to TEM driven turbulence according to TGLF simulations. This predicted switch from ITG to TEM was also observed in the reversal of the turbulent flow measured by the DBS [597]. None of these results can explain what caused the initial change in confinement, but indicate that once a new equilibrium is found that the 3D field geometry might not be the dominant component to the changes in turbulent transport in specific conditions. Similar observations on the role of local gradient are found when studying the impact of magnetic islands and 3D corrugation of magnetic flux surfaces, see sections 6.3.3.1 and 6.3.3.2.

The majority of the research on RMPs is linked to the suppression of ELMs by reducing the pressure gradient in the pedestal region. In the pedestal region, the E × B shear plays an important role in regulating and reducing turbulent transport [549]. In the pedestal region, the E × B shear is dominated by the pressure gradient and as a consequence any changes to the electron density and temperature profiles, will directly affect the E × B shear. Additionally, the changes in toroidal rotation through neoclassical effects, resonant and non-resonant braking (see section 6.2), cannot be neglected as a contribution to the E × B shear, especially in L-mode where the pressure gradient has a much smaller contribution.

The application of RMPs showed a flattening of the radial electric field in MAST L-mode consistent with changes in both perpendicular and parallel flows and a broadening of the ion saturation fluctuations [703]. These effects required the applied RMP field to exceed a minimum value and the effects were not toroidally localized. Similar, in AUG L-mode plasmas, using various RMP configurations also showed a threshold below which no changes to flow, $E_\mathrm r$ and density fluctuations were observed using the DBS [704]. Similar observations were observed in DIII-D L-mode plasmas, where below a threshold changes to the radial electric field and density fluctuations were within error bars [597]. Above the threshold, density fluctuations increased along a substantial change in the toroidal rotation and a change in the turbulent flow as measured by the DBS. Similar changes are observed before the L to H-mode transition in DIII-D, which will be discussed in section 6.4. Both the MAST and AUG experiments showed a differentiation between the turbulence response to RMPs inside the separatrix and in the SOL. This differentiation is attributed to difference in magnetic topology and the creation of a laminar and stochastic zone. This strong dependence of localized turbulence measurements on the underlying magnetic field structure highlights the need, when studying the impact of various RMP spectra in plasmas with different densities and q₉₅ values, to calculate the plasma response. A comparison between L- and H-mode plasmas revealed that in H-mode, even for very small levels of RMPs, a density pump-out was observed, counter to the need for a threshold in L-mode [709]. Aside from potential effects of the plasma response on the strength of the applied field, it was noted that H-mode plasmas are close to marginality and a small change in the E × B shearing rate (through for example neoclassical effects) could have a nonlinear effect on turbulence and transport, through a self-reinforcing cycle until a new transport equilibrium is achieved.

Further examination of how RMPs impact DIII-D H-mode plasmas shows that the E × B shearing rate at the top of the pedestal and inward $\Psi_N \sim 0.7$–0.9 decreases, along with an increase in density fluctuations measured both by the BES and DBS [600]. The reduction in the E × B shearing rate means that the linear growth rates of the dominant turbulent modes are now larger than the E × B shearing rate. This strongly limiting the E × B shearing rate ability to suppress turbulent transport. Gas puff modulations to measure the changes in particle transport indicated that the outward diffusion increases and the inward particle pinch decreases over the whole plasma radius. Identical results were obtained on EAST with n = 2 RMP perturbations and using the supersonic molecular beam to extract the perturbative particle transport coefficients [710]. EAST reproduces that the diffusive coefficient increases at the plasma edge, while the inward pinch decreases. A comparison with neoclassical simulations shows that the experimental obtained values are much larger. Similarly on DIII-D, a reduction of the E × B shearing rate was observed, along with an increase in density fluctuations. Further research on EAST found that there was a toroidal dependence of the E × B shear and fluctuations linked to the corrugation of 3D field when RMPs are rotated around the device [711]. KSTAR found an increase in perpendicular velocity fluctuations when the E × B shear is altered through neoclassical effects [712]. These results were confirmed using gyrokinetic modeling with the GTC code, where they found that changes in the E × B shear were necessary to explain the increase in density fluctuations with 3D fields [713].

Early research on limiter tokamaks such as TEXTOR [596, 714, 715] and TORE-SUPRA [705], and recent stellarator results [697, 716] showed that the underlying 3D geometry is important in how the turbulence characteristics are affected. In diverted tokamak configurations the rotation of RMPs and visible imaging of the X-point region showed that the separatrix is deformed and there is a corrugation [717, 718]. The 3D corrugation will affect long wavelength turbulence and turbulence flows directly. Moreover, the 3D homoclinic tangles [719]—topological structures created by splitting of the separatrix or any other resonant flux surface by resonant components of the magnetic perturbation—close to the X-point result in the creation of a region where closed and open field lines are intermixed. Experiments on TEXTOR studied the ergodic divertor concept, focused on how turbulence is affected by the mixing of chaotic field lines of different lengths. In the ergodic zone in TEXTOR, the turbulent particle flux reverses from outward to inward, which results in an improvement in particle confinement. Along with the reversal of the particle flux, the poloidal phase velocity of the fluctuations altered direction from the electron to the ion diamagnetic direction [596, 720]. The changes in turbulence in the laminar region are linked to a reduced rotational shear, and blobby transport is reduced in the SOL [715]. In the next sections, we will first address the impact of the simple 3D corrugation of magnetic flux surfaces, the impact of island structures, and finally the impact of chaotic field lines.

H-mode experiments in DIII-D in which the 3D-field perturbation was altered showed an immediate change in the density fluctuations as a function of the applied RMP field strength [601]. The increase in density fluctuations was observed over a wide span of q₉₅ values, from mid-radius outwards. This was independent of ELM suppression and the increase in density fluctuations was much larger than the local changes in density. The observed changes in density fluctuations were instantaneous and indicated that the 3D fields had a direct impact. GTC modeling using 3D fields obtained using plasma response model, were not able to reproduce these results [713]. To match the observations, the GTC model in addition to the 3D field required the altered E × B shear to be included. Regardless of the reliability of the physical model used to numerically describe turbulence at the L-H transition, a key aspect of a successful modeling is the study of the effect of the bootstrap current on the local value of the global shear, which has a great impact on the character of turbulence itself [721].

6.3.3.1. 3D corrugation of magnetic flux surfaces.

6.3.3.2. Magnetic islands.

Magnetic islands in the plasma edge are either spontaneously present, due to the MHD activities affecting magnetic flux surfaces with low rational safety factor [736], or result from external perturbation fields interacting with the plasma equilibrium [593]. As a magnetic island breaks the nested magnetic flux surfaces, it modifies plasma transport both locally and globally. The enhanced radial heat diffusivity ($\chi_\mathrm r$) due to the parallel transport along the field lines is increased over a region of about the island width w [737]. The effects are significant if the island width is larger than critical island width [737, 738]

$$\begin{equation*} w_\mathrm c = \left(\chi_{\perp}/\chi_{||}\right)^{1/4} \left(\epsilon_\mathrm s r_\mathrm s/L_{q0\mathrm s} n\right)\end{equation*}$$

where $\epsilon_\mathrm s = r_\mathrm s/R$ is aspect ratio of the resonant flux surface with radius $r_\mathrm s$ and $L_{q0\mathrm s} = q_0/q^{^{^{\prime}}}_{0}$. This section will provide a short summary of the impact of magnetic islands on turbulence and transport, and we refer the reader for a more extensive overview to Choi [739].

Earlier research had shown that heat pulses propagate at different speeds depending on whether they pass through the center of an island or the X-point in both stellarators and tokamaks [740, 741]. This geometric dependence on the heat flux resulted in a series of experiments on J-TEXT [742], DIII-D [743], HL-2A [744] and KSTAR [745] to better understand the origin and investigate the role of turbulence. All devices found a clear increase in density fluctuations using multiple diagnostics that measure density fluctuations at multiple wavelengths. Further research showed that the reduction in turbulence at the O-point is linked to the width of the island and the flattening of the temperature gradient [744, 746]. In general, most devices found good agreement between local gyro-kinetic models, based on the measured gradients, and the turbulence changes, which could include a change from ITG to TEM dominant, depending on the temperature ratio. As such, the magnetic O-point reduced gradients and thus transport, whereas no such reduction in gradients was observed at the X-point, nor a reductions in turbulence characteristics [742–753]. These results would suggest that flow shear and flow shear suppression are less relevant, but not all results fit within this framework. At the outer edge of the island, at the poloidal location of the O-point, turbulence is strongly suppressed, even with the reappearance of a temperature gradient. Experiments on KSTAR revealed a strong poloidal flow [754]. This agreed well with observations of radial electric field reversals in islands observed on LHD [755, 756] and TJ-II [757, 758]. Further experiments on KSTAR found that plasma turbulence becomes strongly inhomogeneous around a magnetic island due to the combined effect of the pressure gradient (driver) and flow shear (regulator) modifications by the island [739].

On LHD Jakubowski et al [759] showed that the natural $m/n = 1/1$ magnetic island enlarged with RMPs leads to an increased particle diffusion coefficient in the vicinity of the island. Temperature (L_T) and diffusion (L_n) scale lengths showed strong variation with the island width for the island width (w_t) larger than 5–6 cm. These findings are in line with gyrokinetic simulations of the static island by Bañón Navarro et al [690], showing a threshold in the island width significant for the transport changes. There are two distinct regimes depending of the island full-width, with a threshold occurring at ${W}_{t} = 33\rho_i$, where ρ_i is the ion gyroradius. Both turbulent radial particle and heat fluxes are reduced at the O-point while only the total heat transport is increased across the island region. The hysteresis relation between density fluctuation ($\tilde{n}/n$) and temperature modulation during the heat pulse propagation into a magnetic island at DIII-D provided experimental evidence of the turbulence spreading from the X- point to the O-point of the magnetic island [755].

Historically, a threshold in power has been observed to trigger the plasma to transition from L- to H-mode. The threshold is an important parameter which allows the identification of power requirements for future fusion reactors as well as acting as a simple metric to assess whether the L-H transition is affected by a change in plasma conditions and provide a scaling metric across multiple devices to aid in uncovering common physics.

The application of RMPs in L-mode, before the L-H transition occurs can result in an increase of the power threshold, which has been observed on DIII-D [760, 761], MAST [762, 763], AUG [764–766], NSTX [767] and, KSTAR [768]. All tokamak experiments indicated that the power threshold increases with increasing resonant magnetic field strength as calculated using plasma response models. This would mean that the power threshold is less affected if the applied field is mostly non-resonant [760–763, 766] and on K-STAR they did find that the power threshold in general, due to smaller error-fields was below those from DIII-D [768]. DIII-D also found that even with similar momentum injection, RMPs applied in ECH heated plasmas increased the power threshold more than RMPs applied in NBI heated plasmas [761]. Moreover, the applied RMPs further amplified the increased power threshold by a factor 2 at low collisionality, with little or no change at high collisionality. MAST and AUG were able to show that the power threshold is directly linked to the magnetic field as calculated using a plasma response model, instead of a vacuum calculation [763, 766]. Devices that employ a plasma response model to calculate the radial magnetic field strength all show a minimum requirement in the applied 3D magnetic field amplitude to increase the L-H transition power threshold. To better understand how the RMP spectrum affects the power threshold for the L-H transition the focus shifted to

• (i)  
  investigating the impacts of RMPs on L-mode plasma conditions before the transition and comparing those to physical parameters that are known to affect the L-H transition without the application of RMPs,
• (ii)  
  investigating how RMPs could alter the self-organization of turbulence through long-wave length turbulence flows which can facilitate the L-H transition.

Both observations and theoretical understanding will be discussed in the next sections with respect to both the impact of altered plasma conditions before the L-H transition as well as the changes during the transition.

6.4.1.1. RMP impact on pre L-H transition plasma conditions.

6.4.1.2. Plasma edge flows.

6.4.1.3. Plasma edge profiles.

The interplay between turbulence, main plasma flows and turbulent plasma flows, such as ZFs and GAMs are critical for the L-H transition dynamics. Two sets of experiments on DIII-D showed that the RMP induced changes in turbulence before the L-H transition depend on whether the RMPs were resonant or non-resonant [761, 770]. Both experiments observed that the application of RMPs result in a strong increase in density perturbations and a shift to lower frequencies, whereas the NRMPs slightly damp the density fluctuation measurements as measured by the BES at multiple radial locations at the plasma edge. These observations were found to be independent of the 3D corrugation of the separatrix. The changes were not consistent with the changes in turbulence drive as calculated with a quasilinear gyrokinetic code, thus indirectly pointing to a 3D effect on turbulence characteristics. Aside from changes in electro-static perturbations, an increase in magnetic fluctuations is observed on DIII-D with the application of both resonant and NRMPs [761].

6.4.2.1. Direct impact of RMPs on turbulent flows.

The observation of the improved-confinement regimes and the occurrence of ELMs in helical devices confirms the generic character of L-H-type transitions in toroidal confinement. The first observation of an L-H-type transition in a 3D device was reported in 1991 at W7-AS [772]. Further results followed by at CHS [773], Heliotron-J [774], LHD [775] and TJ-II [776].

While the improvement of the plasma confinement in the case of 3D devices is not as pronounced as in tokamaks with poloidal divertors, it occurs in a somewhat similar fashion. Indeed, density fluctuations can be associated with specific changes in radial electric field shear and perpendicular velocity at the plasma boundary, and a formation of a pressure pedestal (however, not as pronounced as in the case of tokamaks), which leads to an observable increase in plasma density and plasma energy content. The width of the pedestal (typically a few cms) seems to be independent of the machine size, which suggests that the barrier formation is defined by processes from the plasma boundary independent of the size of the core plasma [777]. The typical MHD activity associated with the L-H transition, manifesting as ELMs, has also been observed in W7-AS [778], LHD [779] and also TJ-II [780]. The details of the ELMy regimes in 3D devices, repetition frequency and associated particle and power loss to the SOL and ultimately to the divertor targets, depend sensitively on the configurations.

Obviously, if the transition to a high confinement regime is due to turbulence suppression, the plasma would ideally transition to the lowest possible level of transport, which in general is the neoclassic one (surely for tokamaks and optimized stellarators). Therefore, the most important difference between tokamaks and stellarator comes from the presence of a 3D magnetic equilibrium, which often impacts the transition to an H-mode due to an ambipolar radial electric field given by the neoclassical fluxes. This is due to the fact that neoclassical transport to large degree depends on the magnetic equilibrium in 3D configurations [603], see section 6.2.

6.4.3.1. Magnetic configuration.

6.4.3.2. Poloidal rotation and shear at the edge.

Another quantity that shows magnetic configuration dependence during the L-H transition is the shear of the negative radial electric field, which can increase in this phase. The shear of the negative radial electric field increases at the L–H transition by an amount that depends on the magnetic configuration and heating power. At TJ-II, fluctuation measurements show a reduction in the turbulence level that is most pronounced at the position of maximum $E_\mathrm r$ shear [782]. However, as opposed to the tokamak case, the poloidal rotation is non-negligible in the radial force balance equation. High temporal and spatial resolution measurements of the radial electric field and density fluctuations at TJ-II indicate that turbulence reduction precedes the increase in the mean shear flow, but is simultaneous with the increase in the low frequency oscillating shear flow. These observations are consistent with L–H transition models, developed in a tokamak context, predicting plasma bifurcations triggered by ZFs [783].

Similarly as in W7-AS [778], in TJ-II the classical H-mode occurs at high densities with NBI heating and changes in plasma radial electric field, plasma rotation and density fluctuations. An example is presented in figure 27. During an L-H transition at LHD in co- and counter-NBI heated discharges E_r is reduced by a factor of 2, which leads to significant increase in negative E_r-shear (figure 27(a)) and the perpendicular rotation velocity. The fluctuation level (figure 27(b)) shows a factor of five reduction at the location of the maximum velocity shear.

Zoom In Zoom Out Reset image size

The application of RMPs results in an increase in transport and a decrease in the electron density in tokamaks. This additional outflow of particles and the creation of a tangle structure interacting with the wall can alter the neutral density and thus fueling of the pedestal. To understand where the particles in this density pump-out wind up, a careful particle balance is performed. Unterberg et al [795] showed with a particle balance analysis that depending on the shape of the plasma in the divertor region, the particles can end up in the cryo-pump or in the wall and SOL and that the pumping efficiency can change [796]. While in low triangularity H-mode discharges in DIII-D the pumping decreases, in the ISS H-mode discharges in DIII-D, pumping increases and remains elevated throughout the RMP application. The changes in pumping were attributed to the distance of the traditional strike-point to the pump plenum entrance and the complex 3D geometry introduced by strike-point striations. Additional measurements show that at the midplane D_α increases by a factor of 2 when compared to the RMP-off phase of the experiment, which can affect the neutral fueling rate of the core plasma [797]. All these changes in particle balance suggest a change in plasma interaction with the walls. Additional modeling using SOLPS showed that while changes in recycling and pumping can lead to a strong reduction in pedestal density, a change in transport was necessary to explain the changes in the pedestal density gradient [798].

As discussed in earlier sections, RMPs can lead to enhanced turbulent and neoclassical transport reducing pedestal gradients both when ELMs are suppressed or only mitigated. To explain the difference between conditions in which ELMs are suppressed versus mitigated, based on the strong dependence on the toroidal rotation as well as the location of the resonant surfaces and the measured plasma response as discussed in section 6.2.2, observations indicate the formation of islands prevents the expansion of the pedestal to an unstable width thus suppressing ELMs [661, 799, 800]. The formation of magnetic islands in the pedestal foot can cause further density pump-out [786]. The induced changes to magnetic topology then impact turbulent transport as discussed in sections 6.3.3.2, 6.3.3 and 6.3.3.1. Nonlinear MHD modeling of KSTAR plasmas with RMPs, where the neoclassical transport code PENTRC has been coupled to the nonlinear 3D MHD code JOREK showed that a combination of kink-peeling, tearing response, and NTV (see section 6.2) can explain the pedestal degradation observed in experiments [801].

Experiments at DIII-D with low pedestal collisionalities ($\nu^* \lt 0.1$) and ITER-similar plasma shapes show clear variation of electron temperature and electron density at the plasma edge [784]. The coherent nature of the oscillations in this region is suggestive of a kink-like helical deformation of the plasma flux surfaces by the RMP. As described in section 6.2.2, the magnitude of the response is limited to certain q₉₅ resonant windows due to the density-dependence of the field penetration threshold. Those windows overlap with the ELM suppression q₉₅ windows, but especially for the pedestal density the window is wider than for the ELM suppression. In regions with high resistivity and and/or low electron fluid velocity, resistive MHD codes allow for penetration of the RMPs to the rational surfaces and thus formation of magnetic islands. EPED modeling [799] suggests that a magnetic island or stochastic layer caused by the RMPs may enhance transport near the top of the pedestal, which results in broader pedestal width. This model is consistent with the resonant nature of the changes in the pedestal densities and temperatures. Higher amplitude of RMPs were found at DIII-D [786] to cause the pedestal collapse in low confinement H-mode plasmas. It occurred when the applied magnetic perturbations exceeded the penetration threshold in the steep pressure gradient region of the ETB. Nonlinear MHD indicated that penetration in the steep gradient region of the pedestal can occur at a comparable RMP amplitude required for ELM suppression in model ITER equilibria.

Prior nonlinear MHD modeling showed an increase in radial convective E × B plasma transport to occur in the presence of the RMPs using nonlinear MHD modeling [802, 803]. Nonlinear MHD modeling [804] shows that the reduction in pedestal density observed in RMP experiments is consistent with the presence of E × B convection cells crossing the separatrix. Such a processes can operate when the resonant field amplitudes are reduced by plasma rotation, although the amplitude of E × B flows is reduced with the increasing plasma rotation velocity at the separatrix. Above some critical rotation, the enhanced density transport is completely suppressed, even though the E × B flow is only partially suppressed. This could mean that the rotation not only leads to screening of the penetration of the external fields, but also plays a secondary role in the reduction of the transport.

Application of RMPs leads very often to scenarios in which ELMs become more frequent and smaller in size. One of the assumptions on such scenarios was that the enhanced particle transport due to much higher frequency of the mitigated ELMs leads to the density drop at the pedestal [805] and would be detrimental for ITER operating at such scenarios.

In fact, as shown in figure 29, the changes in the density on AUG during experiments with n = 1 magnetic perturbation cannot be explained by more frequent mitigated ELMs. Absolute particle losses for natural and mitigated ELMs are identical (figure 29(a)) and depend only on the ELM frequency. At ASDEX Upgrade the particle loss rate ($\Delta n_\mathrm e\cdot f_{\mathrm{ELM}}$) is constant and does not depend on plasma density (figure 29(b)). Measured particle losses suggest that the higher the frequency of ELMs the smaller the particle losses per ELM are. In fact particle losses per second remain more or less constant independent of mitigated ELMs frequency. The application of the RMPs produces an increase in $f_{\mathrm{ELM}}$ and corresponding decrease in $\Delta n_{\mathrm{ELM}}$ consistent with $f_{\mathrm{ELM}}\cdot\Delta n_{\mathrm{ELM}} = \mathrm{const}$.

Zoom In Zoom Out Reset image size

The experiments with plasmas in connected double null or lower single null plasma shapes at MAST showed very similar results [574]. The frequency of the mitigated ELMs in these experiments was strongly related to the amplitude of the maximal resonant component of the applied field. For both shapes, the application of the RMPs lead to an increase in ELM frequency (up to ca. 600 Hz) and the corresponding decrease in $\Delta W_{\mathrm{ELM}}$ was consistent with $f_{\mathrm{ELM}} \cdot \Delta W_{\mathrm{ELM}} =$ 500 kW. The particle losses due to mitigated ELMs were smaller than for the natural (non-mitigated) ELMs. The relative loss of pedestal density due to mitigated ELMs ($\Delta n_\mathrm e / n_{\mathrm{ped}} \approx 0.02$) was two times smaller than for the natural ELMs ($\Delta n_\mathrm e / n_{\mathrm{ped}} \approx 0.04$). Similar effects of reduction of the ELM induced particle and energy losses due to mitigation of ELMs with n = 1 RMPs at LHD have been also reported [807]. Time resolved electron density profiles at the pedestal show much smaller modulation due to mitigated than due to non-mitigated ELMs. Based on these findings one may conclude that any changes in the pedestal height or width come from transport mechanisms not related to particle transport due to ELM losses.

As noted above, the ITPA database studied in [1] demonstrated a minimum in $P_\mathrm{th}$ as a function of density on several devices. The $n_\mathrm{th,min}$, was estimated on each machine in the ITPA database, and only L-H data having $n \gt n_\mathrm{th,min}$ were used to generate the scaling in equation (20). Inspection of $n_\mathrm{th.min}$ data in [1] suggested this critical density may have a positive dependence on $B_\mathrm T$ and a negative dependence on machine size R. Subsequently, detailed examination of the density dependence on multiple machines showed consistent minima in $P_\mathrm{th}$ and the prominence of a ‘low-density branch’, in which the threshold power decreases with increasing density [760, 811–813]. These results confirmed this leading order departure from the ITPA scaling. Specifically, C-Mod [812] and JET operated with its ITER-like wall [813] individually identified a positive, nearly linear, scaling of $n_\mathrm{th,min}$ with B_T. A tendency of $n_\mathrm{th,min}$ to increase with plasma current $I_\mathrm P$ has been observed on JET-ILW [814] and AUG [815].

Much of the leading order behavior of $n_\mathrm{th,min}$ can be captured by a model expression developed by Ryter et al [815] (and later generalized by Bilato et al [816]), inspired by evidence on AUG that a critical surface-integrated edge ion heat flux $Q_\mathrm{i,edge}$, rather than a total power, must cross the separatrix in order to facilitate the L-H transition. Based on the assumption that sufficient electron-ion coupling is needed to achieve adequate power in the ion channel, an expression is derived from a simultaneous solution of $P_\mathrm{th}$ in equation (20), the L-mode energy confinement time $\tau_\mathrm{E,L}$ [817], and an empirically determined critical value of $\tau_\mathrm{E,L}/\tau_\mathrm{ei}$, where $\tau_\mathrm{ei}$ is the volume-averaged electron-ion energy exchange time. The resulting scaling,

$$\begin{equation} n_\mathrm{th,min} \approx 0.7 I_\mathrm P^{0.34} B_\mathrm T^{0.62} a^{-0.95} \left(R/a\right)^{\left(0.4\right)}\end{equation} \tag{ 21 }$$

was found to agree relatively well with both AUG data and those from other devices reported in the literature [1, 812, 813]. Figure 31 illustrates the level of agreement with existing data and some indication of its potential predictive ability, in this case for ITER. Equation (21) has also been used to project the optimum density for accessing L-H transitions on SPARC [810].

Zoom In Zoom Out Reset image size

Other explanations have been proposed for the transition between low and high density branches. For example, numerical modeling provides indications that a critical ion pressure gradient at the plasma edge is needed to gain H-mode access, and that this is related to generation of ZFs and mean flow shear. Collisionality-dependent damping of ZFs may play a role in determining the shape of the $P_\mathrm{th}$ versus $n_\mathrm{th,min}$ curve. Power threshold modeling which accounts for these processes finds that the calculated $P_\mathrm{th}$ attains a minimum at low density and maximum ion heating, but increases at low density and increased electron heating [818] (see discussion in section 7.7. Other modeling work has identified changes in dominant turbulent modes in the edge in a collisionality scan, with a transition from ITG and TEMs toward resistive ballooning modes (RBMs) as density passes from the low to high density branch [819]. In the simulation this results in a changing critical edge temperature for turbulence stabilization and a local transition threshold. Yet another model, partially tested on JET, conjectures that L-H transitions occur when perpendicular transport time is comparable to the parallel Alfvénic time, and the transition between branches of $P_\mathrm{th}$ is explained by asserting the scrape-off layer transitions between sheath-limited and conduction-limited [820]. Models such as these offer tests that could be performed on multiple devices with the aim of deriving additional physics-based predictions of $n_\mathrm{th,min}$.

The sensitivity of the power threshold to main ion species is an important consideration for burning plasma devices. Most experience determining $P_\mathrm{th}$ today is based on DD plasmas, and multi-machine scalings all have this constraint built in (e.g. [1]). Naturally, burning plasma demonstration devices and pilot plants will instead target nearly balanced mixtures of D and T in order to optimize fusion gain. In addition, new burning plasma facilities can require initial operation in H or He plasmas, in order to gain experience without generating excessive levels of neutron fluence. This has been a particularly prominent feature of ITER, which prompted an ITPA Joint Experiment, TC-4, as well as experiments on individual facilities, to explore the main ion species effects on several tokamaks [808]. Examples of $P_\mathrm{th}$ observed in density scans with varied main ion species are shown in figure 32 for AUG and JET operating with metallic walls. Capturing the trends with density is found to be necessary due to the potential variability of $n_\mathrm{th,min}$ with the choice of main ion.

Zoom In Zoom Out Reset image size

First we consider the hydrogenic isotope scaling of $P_\mathrm{th}$. When applying $P_\mathrm{th}$ scalings to DT plasma scenarios that employ significant amounts of tritium, it is customary to apply a $1/A_\mathrm{eff}$ isotope scaling, where $A_\mathrm{eff}$ is the average atomic mass of the hydrogenic species. This implies a reduction in $P_\mathrm{th}$ of $\frac{2.0}{2.5} = 0.8$ when projecting DD $P_\mathrm{th}$ results onto a 50:50 DT scenario. This approach has been in common use dating back to 1990’s threshold experiments in H, D, T and DT plasmas on carbon-wall JET [823], and has found support through additional isotope studies on a variety of devices comparing H and D plasmas [765, 811, 824–827].

Ongoing comparisons on DIII-D, exploiting a carbon wall, continue to show substantially larger $P_\mathrm{th}$ in H plasmas than in D plasmas [760, 824, 825]. There are, however, indications that this gap in $P_\mathrm{th}$ closes at higher density values [827]. This density-dependence variation complicates the goal of scaling $P_\mathrm{th}$ with isotope mass, and is a widespread observation, as will be seen below.

A significant increase in $P_\mathrm{th}$ when going from D to H plasmas was confirmed on AUG with a tungsten wall, as illustrated in figure 32, top. AUG studies indicate that higher edge heat transport in hydrogen plasmas results in the higher threshold power [828–831], and that this is also reflected in a higher edge ion heat flux observed at the L-H transition in H plasmas, when compared with D plasmas [828, 831, 832]. $P_\mathrm{th}$ data from He plasmas are shown in the figure as well, and will be discussed later.

JET has made comparisons of $P_\mathrm{th}$ in H and D plasmas operated with W/Be PFCs [821, 833]. Figure 32, middle, shows trends for JET-ILW $P_\mathrm{th}$ in a specific configuration (horizontal strike point target, 1.8 T, 1.7 MA) comparing H, D and He plasmas [821]. A strong dependency of $n_\mathrm{th,min}$ on plasma species is seen. For densities below $n_\mathrm{th,min}$(H), $P_\mathrm{th}$(H) is significantly higher than $P_\mathrm{th}$(D), and $P_\mathrm{th}$(He) is much larger than $P_\mathrm{th}$(D) and $P_\mathrm{th}$(H) at medium densities (see the more detailed discussion on helium below). The difference in $P_\mathrm{th}$ between H and D plasmas in JET-ILW is much larger for NBI-heated plasmas, reaching up to a factor of 4. Investigations are ongoing to determine whether this is due to an effect of NBI torque on $P_\mathrm{th}$ (as observed on DIII-D [825]), or to the partition between ion and electron energy fluxes at the L-H transition (as observed for example on AUG [834]), or to radiation from Cu influx in NBI-heated plasmas.

Recent JET-ILW studies of $P_\mathrm{th}$ during the Tritium and DT campaign [822] illustrate the isotope dependencies of $P_\mathrm{th}$ in a higher field configuration (horizontal strike point target, 3.0 T, 2.5 MA), for which it is relatively easy to document the change in $n_\mathrm{th,min}$ and $P_\mathrm{th,min}$ with isotope. Shown in figure 32, bottom, is $P_\mathrm{th}$ as a function of density for D, DT and T plasmas. The value of $n_\mathrm{th,min}$ is lowest for T, followed by DT then D plasmas. A strong increase of $P_\mathrm{aux,min}$ with $A_\mathrm{eff}$ is observed, suggesting it may be easier to enter H-mode in T-rich plasmas in future DT fusion devices and experiments. In the high density branch thresholds increase monotonically with decreasing isotope mass, with $P_\mathrm{th}$(T) $\lt P_\mathrm{th}$(DT) $\lt P_\mathrm{th}$(D), yet a quantitative study of this dependency based on the new data has not been undertaken yet.

Altogether the accumulated results from various machines reveal no simple dependence on isotope mass, in scans of hydrogenic composition. The $1/A_\mathrm{eff}$ scaling for $P_\mathrm{th}$ should be considered a rule of thumb with significant uncertainty, with actual scalings departing quantitatively on each machine. TCV for example reports a H:D $P_\mathrm{th}$ ratio of 7:4 rather than 2:1 [826], and AUG has shown similar departures [828, 831]. Further, the density dependencies noted above complicate the ability to scale $P_\mathrm{th}$ with $A_\mathrm{eff}$, as does the variability with divertor geometry and input torque (discussed more below).

Initial mixed isotope campaigns on JET began to deliver data on $P_\mathrm{th}$ in tritium-containing plasmas [836–838]. In those experiments Ohmic T H-modes at low density agreed reasonably well with a $P_\mathrm{th}\sim 1/A_\mathrm{eff}$ scaling, but ICRH-heated plasmas at higher density did not. This was, in part, explained by a hardware problem in the RF system [837]. New measurements of $P_\mathrm{th}$ in plasmas of pure and mixed species, taken with good RF conditions, are reported in [838]. A comparison of various isotope mixtures in this study, including DT and HT, with pure isotope plasmas revealed that $P_\mathrm{th}$ can vary up to 65% for the same effective isotope mass $A_\mathrm{eff}$ at otherwise constant parameters [838]. It is found that different L-mode transport levels are present with different isotope compositions having the same $A_\mathrm{eff}$ (e.g. DD and HT), necessitating different power levels to reach the critical edge temperatures as needed for H-mode access.

Other isotope dependence studies have been performed in mixed hydrogen-deuterium plasmas. Transitions in 50:50 H-D mixtures in JET-ILW resulted in $P_\mathrm{th}$ and $n_\mathrm{th,min}$ values intermediate between the pure hydrogen and deuterium cases [833]. It was further concluded that, at least in the high density branch, $P_\mathrm{th}$ variation in HD mixtures was non-linear, with significant threshold variation at less than 20% or more than 80% H concentration, with little variation in between [839]. A different variability was observed in AUG HD mixtures: $P_\mathrm{th}$ and the edge ion heat flux are constant at the level of pure deuterium for $n_\mathrm H/(n_\mathrm H+n_\mathrm D)$ between 0 and 0.5 and monotonically increase to the level of pure hydrogen for $n_\mathrm H/(n_\mathrm H+n_\mathrm D)$ between 0.6 and 1 [831].

Lastly we consider helium plasmas, which are potential targets for early operation of burning plasma devices, as they may be useful in creating operational experience prior to introducing high neutron fluxes. Introducing He into tokamaks has been seen to influence $P_\mathrm{th}$ in a variety of ways. Often it is seen that He plasmas have a significantly larger $P_\mathrm{th}$ than in D, and prior to about 2008, indications from experiments were that $P_\mathrm{th}$ in He plasmas was simply about 40% higher than that in D [811, 840]. Subsequently, as part of the ITPA TC-4 Joint Experiment, the ratio of the thresholds $P_\mathrm{th,He}/P_\mathrm{th,D}$ was found to range between 1.0 and 1.8 on C-Mod, DIII-D, JET, MAST and NSTX, and was shown to depend on the choice of density on a number of these devices [808]. Similar findings were later reported on EAST [841]. Notably, the discrepancy between $P_\mathrm{th}$ in He and D is often seen to decrease at higher density values [808, 841, 842]. An illustration of this can be seen in figure 32, middle. On JET, for densities below the $n_\mathrm{th,min}$ for He, $P_\mathrm{th}$ is significantly higher for He than for D, but at higher density there is little difference [821]. Although less pronounced, a similar effect was observed in DIII-D (see [760], figure 1). The observation that in JET-ILW $n_\mathrm{th,min}$(He) is larger than $n_\mathrm{th,min}$(H) results in similar minimum $P_\mathrm{loss}$ for both species (partly due to enhanced W sputtering and higher radiation in He). Surprisingly, although AUG found $P_\mathrm{th}$ to be higher in He than in D with its carbon wall, following the transition to tungsten PFCs, the He power threshold was similar to the one in deuterium plasmas, even exhibiting the same non-monotonic density dependence [765, 811]. This is shown in figure 32, top.

Given the higher $P_\mathrm{th}$ commonly associated with H plasmas, and the potential desire to start up next step devices in hydrogen, experimental attempts have been made to lower the threshold with He doping. The results range from significant to negligible. In JET-ILW, helium-4 fueling into H plasmas was performed, resulting in a ∼25% reduction of the threshold with up to about 10% He concentration [839]. DIII-D shows a similar level of reduction in $P_\mathrm{th}$ but requiring larger concentrations of He (${\sim}30\%$) [843]. On AUG $P_\mathrm{th}$ stays constant on the level of pure hydrogen in hydrogen plasmas with helium doping up to $c_\mathrm{He} = 20\%$ [834, 844].

When considering extrapolation of $P_\mathrm{th}$ scalings to burning plasma devices, it is important to consider additional qualities which distinguish these future machines and which are not universally represented on present devices. In particular, data used to generate scalings such as equation (20) originate from multiple devices, many of which have features not expected for ITER, DEMO, SPARC, and fusion pilot plants. Examples include (1) significant torque driven by NBI and (2) plasma facing components made predominately of carbon.

In order to determine the influence of externally applied NBI torque on $P_\mathrm{th}$, the balance of co-current and counter-current beams was adjusted on power threshold experiments on DIII-D [825, 845]. Examples of $P_\mathrm{th}$ response to input torque are illustrated in figure 33 for several different main ion species plasmas. H-mode access was facilitated in all cases studied (favorable vs. unfavorable $B\times\nabla B$ direction, different main ion species) by having externally applied torque be more counter-current. At zero input torque, power threshold is modest and the L-H transition may be facilitated by an increased poloidal flow shear. These results suggest $P_\mathrm{th}$ should not be worse in a low-torque burning plasma device.

Zoom In Zoom Out Reset image size

Presently used scalings for $P_\mathrm{th}$ such as equation (20) were developed using data from mostly carbon wall machines, creating questions about their applicability to metal wall machines. Subsequently, two devices reported changes in $P_\mathrm{th}$ following changes from carbon PFCs to metal PFCs. AUG showed an approximately 25% reduction in $P_\mathrm{th}$ following changeover from carbon to tungsten [765]. This reduction was seen for both D and H independently, while the reduction for He discharges was much more significant. The AUG data in figure 32 were collected with W PFCs, and for D plasmas with $n \gt n_\mathrm{th,min}$ the $P_\mathrm{th}$ values are notably smaller than the ITPA scaling would predict. JET also compared its H-mode access with its ITER-like wall (W, Be) to prior results with C PFCs [813]. At plasma densities above $n_\mathrm{th,min}$, $P_\mathrm{th}$ in D plasmas was reduced by ∼30%, and by ∼40% when the radiation from the bulk plasma is subtracted, when going from JET-C to JET-ILW. One can see the prevalence of JET-ILW deuterium $P_\mathrm{th}$ values below the ITPA scaling curve in figure 32. Similarly, $P_\mathrm{th}$ data from Alcator C-Mod, which operated with molybdenum PFCs for its entire existence, showed optimal $P_\mathrm{th}$ significantly below the values given by equation (20), as seen in figure 34. Possible explanations for the reduction in $P_\mathrm{th}$ with C replacement by W are described in section 7.5.

Zoom In Zoom Out Reset image size

The significant differences between carbon and metal walls with respect to the H-mode power threshold have suggested the formulation of multi-machine scalings that include only devices with low carbon content. An ITPA activity (TC-26) has been initiated to this end, examining data sets from AUG, C-Mod and JET. The complete results of the database analysis are to be published, although the subset of data from JET-ILW has been used to produce an alternate scaling [836]. Beginning from the form of the expression in equation (20), a linear relationship between $P_\mathrm{th}$ and plasma surface area S was assumed and dependencies on field and density fitted to data for which $n \gt n_\mathrm{th,min}$. In aggregate a stronger density dependence is found than in the 2008 ITPA scaling, along with a striking separation in terms of divertor geometry. When JET-ILW is operated with a horizontal divertor target, $P_\mathrm{th}$ is systematically below the scaling prediction, as for the D points in figure 32. However when the outer strike point is placed on or near the vertical target plate, elevated $P_\mathrm{th}$ is observed. Direct comparisons of the JET-ILW threshold to the 2008 ITPA scaling database are shown in figure 34. Divertor effects such as this turn out to be not uncommon on tokamaks, as further discussed in section 7.3.4.

As tokamaks move toward demonstrating net fusion power, the emphasis on divertor power handling becomes crucial, and reactor concepts must develop improved divertor designs, sometimes incorporating advanced divertor concepts. In addition future concepts move increasingly farther from single null equilibria and toward double null configurations to facilitate exhaust management. As noted above, the specifics of the divertor geometry, including the placement of strike points and X-points, can have a sizable impact on the power threshold for H-mode on a given device. This was recognized at the time of the Update to the IPB [3], based in part on $P_\mathrm{th}$ comparisons in JET with and without a divertor septum and with scans of X-point height [847]. There a decrease in $P_\mathrm{th}$ was observed as the X-point moved closer to the divertor floor, while strike points were on mostly horizontal targets. Similar trends have also been reported on DIII-D with the outer strike point on a horizontal target, although equilibrium shaping was also varied [760]. MAST also performed experiments demonstrating a substantial reduction of $P_\mathrm{th}$ as the X-point was lowered [848]. Subsequently, L-H transitions on TCV under Ohmic heating showed a clear trend of reduced $P_\mathrm{th}$ with X-point closer to the vessel floor [849].

While the above experiments indicated that longer outer divertor leg length correlated with higher H-mode threshold power, an inverse trend was reported on Alcator C-Mod [850]. However, the substantial decreases in $P_\mathrm{th}$ on C-Mod corresponded also to the movement of the outer strike point from a vertical to a horizontal target, and also to a configuration with greater neutral baffling [846] (see figure 34). These results in addition to the JET-ILW results with vertical vs. horizontal targets [836], may indicate a greater importance for divertor shaping in setting $P_\mathrm{th}$. Despite the notable influence of divertor geometry and strike point placement on many devices in at least one configuration on TCV, little effect was seen on $P_\mathrm{th}$ from variation in both the major radius of and the poloidal flux expansion at the strike point [851].

Systematic studies of the effects of divertor configuration and placement of X-points and strike points have been attempted under Joint Experiments originating in the ITPA PEP Topical Group. The first was PEP-28, ‘Physics of H-mode access with different X-point height’, inspired mainly by the JET, DIII-D and MAST results obtained prior to 2011. Increasingly, as findings emerged on these and other devices suggesting complex dependencies on divertor geometry and plasma shaping, PEP-28 transitioned to PEP-39, ‘Impact of the SOL and divertor conditions on the L-H threshold’. Work continues in this area on multiple devices, with a focus on single null discharges with the $\nabla B$ drift in the favorable direction for H-mode access, i.e. with $B\times\nabla B$ directed toward the active X-point.

Tokamaks commonly observe that power threshold in SN discharges is significantly elevated when the direction of $B\times\nabla B$ is reversed to point away from the active X-point. Typical experience suggests a factor of ${\sim} 2$ increase [852–854]. An interesting exception was noted in JET, in which $P_\mathrm{th}$ for accessing H-mode was seen to agree relatively closely across a $B_\mathrm T$ direction switch [855]. Notably, operating JET in the unfavorable direction did, in some cases, increase the power requirements to reach an H-mode with an improved density pedestal and large ELMs. Some other devices have associated similar phenomena with I-mode (see section 7.4), a confinement regime that is distinct from L-mode and H-mode, and which tends to be observed at intermediate power levels when operating with unfavorably directed $B\times\nabla B$. In scans of the magnetic balance on various devices, increases in $P_\mathrm{th}$ can be observed in SN plasmas that are near the DN point, and transitioning from the favorable to unfavorable case [852, 856, 857]. However, there is also evidence of improved H-mode access in nearly exact DN configurations, as first reported on MAST [858] then in a joint study among MAST, NSTX-U and AUG [859], and in experiments on EAST [860]. Power reductions in connected DN configurations on the low aspect ratio devices were found to be a factor ${\sim}2\times$ while on AUG the $P_\mathrm{th}$ decreased by only ≈30% in DN.

Details of equilibrium shaping are not captured by typical $P_\mathrm{th}$ scalings, other than indirectly through the surface area dependence (see equation (20)). Direct tests of the influence of plasma triangularity δ on $P_\mathrm{th}$ have been attempted, although design of such is complicated by the correlation of plasma equilibrium shape changes with X-point movement and strike point placement within divertors, aspects which as seen above have a large impact on $P_\mathrm{th}$. For example, indications from JET that H-mode access was facilitated at increased δ could not readily be decoupled from transitions from vertical to horizontal divertor targets [861]. Later experiments with the JET-ILW suggested that, for a given δ associated with the lower chamber divertor, increasing the upper δ away from the divertor actually raised the power flux needed to trigger the L-H transition [813]. Increasing power requirements with higher triangularity was also shown in a scan of δ from 0.2 to 0.8 on DIII-D [853].

Whereas most tokamak experience consists of plasmas with non-negative values of triangularity, a growing body of research into negative triangularity (NT) has come about in recent years [862]. NT plasmas have demonstrated improved energy confinement over their positive triangularity analogs without the need to form an H-mode pedestal, as demonstrated on TCV [863–865] and DIII-D [866, 867]. Indeed the L-H transition can be suppressed up to much higher levels of loss power than in positive triangularity [867, 868]. Single null discharges with negative, but nearly zero, average triangularity on AUG exhibit elevated $P_\mathrm{th}$ relative to positive triangularity cases, when $B\times\nabla B$ is in the unfavorable direction, but show $P_\mathrm{th}$ similar to PT cases when the drift direction is favorable [869]. The mechanism behind the suppression of H-mode access for sufficiently negative values of δ is actively being investigated. One promising explanation for the suppression is the tendency of NT plasmas to be limited to low edge gradients by high-n ballooning instabilities, which prevents the formation of significant $\textbf{E}\times\textbf{B}$ shear and pedestal turbulence suppression [870]. This model has enabled prediction of H-mode suppression for a wide range of possible NT reactor parameters [871]. Additional discussion of NT is found in a companion article [21].

Departures from axisymmetry exist in real toroidal devices, and these are known to have effects on H-mode access. One case of special interest for tokamaks is the influence of toroidal field ripple. The finite number of toroidal field coils in a tokamak results in a ripple in $B_\mathrm T$, which is known to influence the degree of ion orbit loss and which could impact flows and radial electric field in the edge of a plasma as it approaches an L-H transition. In order to assess the impact of ripple on H-mode physics, JET utilized a feature of its TF configuration to scan its ripple field from 0.08% to 1.10%, a range intended to span the anticipated ripple field in ITER at its nominal operating field [861]. Although significant changes in toroidal velocity at the edge were observed in the scan, there was no impact on $P_\mathrm{th}$. Subsequently, DIII-D performed experiments with a set of in-vessel coils intended to simulate the level of TF ripple commensurate with the presence of ferromagnetic test blanket module inserts on ITER [760]. A local ripple value of 3% was achieved and $P_\mathrm{th}$ for H-mode was evaluated for both NBI and ECH heated discharges, at low values of input torque, consistent with ITER and other burning plasma targets. The presence of the applied ripple field had no appreciable effect on access to H-mode.

In addition, a wide variety of 3D magnetic perturbations may be applied in tokamak plasmas with the intent to affect features such as core rotation, particle inventory and the prevalence of ELMs in the H-mode pedestal. The influence of applied 3D perturbations has been studied on several devices, including DIII-D [760, 761, 770], AUG [764–766], MAST [762, 763], NSTX [767] and KSTAR [768]. Tokamak experiments generally observe increases in $P_\mathrm{th}$ with increasing strength of RMPs, with indications that non-resonant MPs have less impact on power thresholds [760–763, 766]. Experiments on AUG determined that while the increase on $P_\mathrm{th}$ was maximized for the perturbation alignment planned for ELM suppression on ITER, the radial field perturbation level needed to realize this elevated threshold power actually exceeds the level required for ELM suppression [766]. If the AUG results extrapolate to ITER, access to H-modes with RMP becomes promising. An overview of the effects of RMPs is found in section 6, where discussion is included of the effects of these perturbations on L-mode plasma conditions before the L-H transition, and on how RMPs could alter the edge turbulence and flow characteristics in the vicinity of an L-H transition.

L-H transitions are also observed in toroidal devices with significant intrinsic 3D fields. Details of observations in stellarators are discussed in section 6.

In addition to conducting database studies of engineering parameters to produce global scalings, it is instructive to examine data sets of local edge parameters. Mapping out the operational space of the edge in a given device, and potentially on combinations of devices, can inform the boundaries between confinement regimes in local terms and facilitate comparison with transition models.

The UIPB reported on efforts to compare C-Mod and DIII-D edge data to theoretical predictions for critical parameters at the L-H transition, based upon the finite beta stabilization of drift-RBMs [3]. Extended work on C-Mod examined the operational space in terms of ballooning-normalized pressure gradient $\alpha_\mathrm {MHD}$ and a diamagnetic term $\alpha_\mathrm {dia}$ and found a separation between L-mode and H-mode points in qualitative agreement with predictions from theory, for both Ohmic [885] and RF-heated plasmas [812]. On AUG a clear separation between L-mode and H-mode can be shown directly in the operational space of separatrix electron temperature and density. In this picture, the dominant turbulence is more drift wave like at low densities and more ballooning like at high densities. The transition into H-mode occurs in both cases if the energy transfer from the turbulence to the mean flow exceeds the turbulence drive [886]. Figure 36 illustrates the different operational limits of the AUG L-mode and H-mode domains for $B_\mathrm t = 2.5$ T, $I_\mathrm p = 0.83$ MA, plotted in terms of the separatrix (LCFS) density and electron temperature. Here, in addition to the MHD resistive ballooning limit and the L-mode density limit, the experimentally observed L-mode and H-mode operational domains are delineated well by the calculated criterion for H-mode access, when the energy transfer from the turbulence spectrum exceeds the turbulence drive.

Zoom In Zoom Out Reset image size

There are numerous experimental indications that gradients in edge kinetic profiles influence access to H-mode. Several of these in turn are consistent with the hypothesis that a critical $\textbf{E}\times\textbf{B}$ shear is necessary to trigger L-H transitions, as discussed in section 7.5.2. This is due to the leading dependence of radial electric field $E_\mathrm r$ on the diamagnetic term proportional to the ion pressure gradient. On AUG, at density values high enough to ensure electrons and ions are coupled, the H–L back-transition occurs at the same electron pressure gradients as the L–H transition [887, 888]. The threshold power for the H–L back-transition exhibits the same increase at low densities as the L–H transition, as described in section 7.3.1. At modest to high densities on Alcator C-Mod, a critical edge $T_\mathrm e$ is also observed at the L-H transition [812, 854].

Observations in JET-ILW indicate that, at least in the high density branch, the favorable isotope effect on $P_\mathrm{th}$ is a consequence of the favorable isotope effect on transport at the L-mode edge. At a given line density, $n_\mathrm e$ and $T_\mathrm e$ profiles at the transition for H, D, H and T mixtures, T, and DT plasmas match, despite very different power threshold values and relation to $n_\mathrm{th,min}$ in each case [833, 838]. This picture is also supported by experimental work on AUG [831] and theoretical investigations [829, 830].

The impact of PFC material on $P_\mathrm{th}$, described in section 7.3.3, was examined on AUG through comparison of kinetic profiles. A reduction in $P_\mathrm{th}$ by about 25% was observed after the plasma facing components have been changed from carbon to tungsten [765]. A possible explanation for this reduction in power threshold could lie in the changed density profiles, which are steeper in plasmas with W PFCs [889]. The higher particle and energy reflection coefficient of W when compared to C has been shown to steepen the density profile [890]. Another possible contributor to the reduction in $P_\mathrm{th}$ when going from C to W PFCs is the reduction in low Z impurities and reduced edge $Z_\mathrm{eff}$. Gyrokinetic modeling of pre-L-H conditions on JET-ILW suggests a prevalence of instabilities with growth rates that decrease as $Z_\mathrm{eff}$ decreases, which can be linked to experimentally reduced $P_\mathrm{th}$ at sufficiently high collisionality [891]. In contrast, at low collisionality (higher edge temperature), where ITG/TEM modes are expected to dominate, recent DIII-D experiments and corresponding gyro-fluid/kinetic simulations indicate an upshift of the ITG critical gradient, and hence reduced instability growth rates and thermal fluxes with increasing $Z_\mathrm{eff}$, resulting in reduced $P_\mathrm{th}$ [892]. This could be an important consideration for discharges with impurity seeding for dissipative divertor operation, also.

The impact of sheared $\textbf{E}\times\textbf{B}$ flow on stabilizing plasma turbulence and improving confinement in magnetized plasmas is widely accepted, and is thought to play a key role in the origin of pedestal formation in tokamaks [893]. Typically, the radial electric field well at the edge of the confined plasma shows a distinct minimum, where the minimum corresponds to the strongest $\textbf{E}\times\textbf{B}$ velocity (see for example [549, 894, 895]). Examples of fine scale variation in $E_\mathrm r$ have been observed in some circumstances, as will be discussed in section 7.6. More generally, however, smooth profiles are reported, and observed differences in low-power plasmas are often invoked to explain changes in observed $P_\mathrm{th}$, for example in the case of a reversal in the $\nabla B$ drift direction [896, 897]. Under the assumption of a critical shear in the poloidal rotation from $\textbf{E}\times\textbf{B}$, it is possible to predict scalings for L-H power threshold that are linear in density and field (see for example [898]). Numerous experiments have sought to determine whether a critical value of sheared $\textbf{E}\times\textbf{B}$ velocity can indeed explain thresholds on a given device.

On AUG a threshold in the $\textbf{E}\times\textbf{B}$ velocity minimum, serving as a proxy for the $\textbf{E}\times\textbf{B}$ velocity shear in the pedestal region, has been found at the L-H transition in favorable drift direction plasmas with a single X-point geometry [832, 899]. Experimental indications are that, for discharges of a sufficient density (i.e. above the low density branch), the plasma transitions into H-mode if this critical velocity is exceeded. Under these conditions the evolution of the $E_\mathrm r$ near the L-H transition is dominated by the $\nabla p_i/n$ term [889, 897, 899], and emphasis is placed on the accurate identification of the ion pressure gradient in L-mode under various conditions. The $\textbf{E}\times\textbf{B}$ velocity threshold at AUG was found to be highly uniform across a range of toroidal magnetic field strengths, densities and for both D and H isotopes, with $P_\mathrm{th}$ varying by a factor of 3 for this database. Using these findings, a simple derivation of the $P_\mathrm{th}$ scaling is proposed giving a physics interpretation of the $B_\mathrm T$, density, isotope and surface area dependence of $P_\mathrm{th}$ [832].

As discussed in section 7.3.1, a linear dependence on density of the ion heat flux at the L-H transition can explain the non-monotonic density dependence of the L–H threshold power [815]. On AUG, results are persistent in H plasmas, helping to explain the isotope effect on L-H power [828]. Based on this key role of the ion heat channel and the L-mode scaling of the energy confinement time, a phenomenological model for the L-H power threshold was developed that explains the dependence of $P_\mathrm{th}$ on plasma current $I_\mathrm p$ in the low density branch on AUG [816]. In the context of this model it is shown that the assumption of a critical velocity shear (shear suppression paradigm) leads to the prediction of the ion heat flux being linearly dependent on density (key role of ion heat flux).

On the other hand, JET-ILW observations contradict the hypothesis of a critical $\textbf{E}\times\textbf{B}$ edge shear flow at the L-H transition across the range of parameters explored (including variations in density, ion mass, heating method and divertor configuration) [895]. No significant change in the edge $\textbf{E}\times\textbf{B}$ shear flow is measured along the power ramp preceding the L-H transition in D plasmas [836, 900]. In D plasmas, but not in He plasmas, there are indications that the inner shear layer, farther from the separatrix, may play a role in setting the threshold, rather than the outer shear layer [895]. The measured $E_\mathrm r$ profile in the SOL shows only a modest dependence on divertor configurations and cannot explain the different power thresholds with horizontal and vertical targets described in section 7.3.4 (although a significant impact on $E_\mathrm r$ inside the confined region was seen) [900]. These experimental results contrast with earlier findings from modeling in which lower $P_\mathrm{th}$ is correlated with $E_\mathrm r$ spikes in the near SOL [901].

The aforementioned work on the critical ion heat flux for L-H transitions arose from studies of the so-called ‘low density branch’ (section 7.3.1), in which electrons and ions are not well coupled by collisions and in which the heat transport of electrons and ions can be separately determined. For example, this is done in AUG experiments by solely heating the electrons via electron wave heating (ECRH). It could be shown that the ion heat flux at the plasma edge exhibits a linear dependence on the line averaged density, whereas the electron heat flux does not [765, 815]. Examples for two values of plasma currents are shown in figure 37. The linear dependence of the ion heat flux at the L-H transition with density can explain the non-monotonic density dependence of the L–H threshold power: the L-H transition mechanism possibly requires a critical edge ion heat flux and due to insufficient coupling of electrons and ions at low density, a much higher total input power is needed to reach the required edge ion heat flux. With pure ECRH heating the threshold power is also higher at low densities than with additional NBI heating, confirming the role of the ion heat flux across the plasma edge [815].

Zoom In Zoom Out Reset image size

Extended tests of the critical heat flux hypothesis were made by examining Alcator C-Mod discharges heated with ion cyclotron power (ICRH), and calculating the contributions of power flow into the ion and electron channels [902]. Results were found to be consistent with a linear scaling of $Q_\mathrm{i,L-H}$ with density. A fit to a combined data set of AUG and C-Mod data resulted in a power law expression for the critical edge ion heat flux at the L-H transition

$$\begin{equation} q_\mathrm{i,L-H} = 0.0021 \times \bar{n}_\mathrm e^{1.07\pm 0.09} B_\mathrm T^{0.76\pm 0.20}\end{equation} \tag{ 22 }$$

where $q_\mathrm{i,L-H}$ is in units of MW m⁻², $\bar{n}_\mathrm e$ is the line averaged electron density in units of $10^{19}\,\mbox{m}^{-3}$, and $B_\mathrm T$ has units of T. Multiplying both sides of this equation by surface area S yields a power law expression for the surface-integrated ion heat flux $Q_\mathrm i$ with dependencies not unlike ‘nBS’ scalings for $P_\mathrm{th}$ such as equation (20). Figure 38 illustrates the quality of the fit for AUG ECH-heated plasmas and C-Mod ICRH-heated plasmas. Some NBI-heated plasmas on AUG are found to have higher $q_\mathrm{i,L-H}$ than predicted by equation (22). This is also, on average, the case for NBI-heated deuterium discharges on JET [903], also shown in figure 38. The role of torque, discussed in section 7.3.3, provides a possible explanation for this discrepancy. A role for toroidal velocity in determining $Q_\mathrm{i,L-H}$ was identified by Ryter et al [815]. When significant co-current NBI torque is applied, a reduction in the magnitude of $E_\mathrm r$ naturally occurs from the increased $\textbf{v}\times\textbf{B}$ in the radial force balance. This can result in a higher ion heat flow being needed to reach a sufficient $\nabla p_\mathrm i /n$ and critical velocity shear, enabling the L-H transition. However, as noted in section 7.5.2, JET evidence points against a critical velocity shear value defining access to the L-H transition.

Zoom In Zoom Out Reset image size

Analysis of I-mode threshold data on AUG also agreed with the hypothesis of a critical $Q_\mathrm i(n)$ [835], as figure 37 illustrates. In this analysis, the $Q_\mathrm i$ is strictly linear with line density with no offset, while an offset density is needed to fit the I-mode $P_\mathrm{th}$ trend [880]. Independently performed gyrofluid simulations targeted at AUG I-mode edge conditions additionally suggest requirements on local edge conditions, namely low collisionality near the separatrix and an ion-to-electron temperature ratio well above unity [904].

On AUG limit cycle oscillations (LCO, also called ‘I-phase’ (distinct from ‘I-mode’ discussed in section 7.4) are observed when close to the L-H power threshold. At very low density, Doppler reflectometry data are interpreted in a way that GAM shear flow can be dominant in the first phase of the L-H transition giving rise to a turbulence-flow interaction, which generates the LCOs [922]. LCOs were also observed in DIII-D, associated with low frequency (${\sim}2$ kHz) poloidal flow oscillations consistent with the measured Reynolds stress drive, before edge pressure gradient oscillations become dominant leading up to the transition from LCOs to H-mode [905, 913]. By way of example, figure 39 shows contour plots of the normalized density fluctuation level $\tilde{n}/n$ and the $\textbf{E}\times\textbf{B}$ velocity, as measured by a multi-channel Doppler Backscattering system in DIII-D, in comparison to the D_α recycling light evolution [769]. The data demonstrate clearly the periodic formation of the well structure in the $\textbf{E}\times\textbf{B}$ velocity during LCO, starting at the LCFS, and the concomitant periodic reduction in normalized edge turbulence level $\tilde{n}/n$, as well as the radial inward expansion of the edge flow layer. At the transition from L-mode to LCOs, a 90^∘ phase shift between $\textbf{E}\times\textbf{B}$ velocity and density fluctuation level has been observed in several devices, including TJ-II [783] and DIII-D [905]. At later times during LCO, and in a wider density range, plasma flow and pressure gradient are in phase, indicating that the diamagnetic component of the $\textbf{E}\times\textbf{B}$ flow dominates under those conditions. Near the LCO-H-mode transition, LCOs show a behavior more similar to type-III ELMs [926, 927], with magnetic precursors on Mirnov signals sometimes present.

Zoom In Zoom Out Reset image size

In contrast to the DIII-D measurements, the turbulence-driven flow contribution is measured to be small in AUG [928] and in JFT-2M [917, 918], demonstrating the dominant neoclassical contribution to the background flow in the vicinity of the L-H transition. Another interpretation of the LCOs in I-phase stems from measurements of up–down asymmetric parallel current fluctuations occurring during the I-phase in ASDEX Upgrade. They are interpreted as a response to strongly ballooned plasma transport coupled with the Stringer spin-up mechanism [929].

DBS measurements on JET have revealed radial structure in E_r profiles, within the well, in L-mode near the high density H-mode transition. The structures were shown to be quite stationary and were found consistent with ZFs. A conjecture is that stationary ZFs play a stronger role in mediating the L-H transition on JET than oscillatory GAMs or LCO mechanisms [930]. In most JET L-H transitions achieved with slow power ramps, a low frequency n = 0, m = 1, up–down, magnetic oscillation is observed in the H-phase, named the M-mode at JET [931]. A multi-diagnostic comparison of JET and AUG I-phase shows that the M-mode and the I-phase are the same [932] While M-mode frequency has been shown on JET to scale with the poloidal Alfvén velocity, I-phase frequency scaling is more complex [926] and no universally accepted scaling has been established.

The interaction of edge turbulence with edge flows is also significant for H–L back-transitions [906, 909]. Substantial power hysteresis has been observed in DIII-D between the L–H transition threshold power (separatrix loss power) and the required ‘holding’ power below which an H-L back-transition is initiated. These observations contrast recent results in AUG, where a significant L-H/H-L power hysteresis was only obtained if the density difference between the pre-L-H transition L-mode plasma and the pre-H-L transition H-mode plasma was taken into account [897]. In regimes where LCOs are present, the H-L back-transition is typically characterized by transients resembling type I or type III ELMs, although these transients can be triggered at edge pressures well below peeling-ballooning stability limits. Following a significant reduction of the pedestal pressure gradient, LCOs are often observed, with the cycle frequency increasing before edge turbulence levels rise and an LCO-L-mode transition occurs. The time evolution observed in the L-LCO-H-mode progression is therefore reversed in the H-LCO-L back-transition sequence, and a reciprocal evolution of the LCO frequency is observed during the back-transition [926].

Direct measurements of the Reynolds stress gradient across the L-H transition and the energy transfer from edge turbulence into directional poloidal flows have been carried out for example in DIII-D [910], and C-Mod [914]. On DIII-D in edge plasmas nearing the L-H transition, rapid ($\lt$1 ms) increases are measured in turbulence amplitude, Reynolds stress gradient and poloidal flow. This burst is seen to immediately precede the transition, and his highly localized to the region 1–2 cm inside the plasma separatrix. Figure 40 shows examples of this profile evolution. The connection of this to the L-H transition is further reinforced by the infered non-linear energy transfer from the turbulence into ZFs. Prior to the L-H transition, the normalized energy transfer rate was shown to exceed both the effective turbulence growth rate $\gamma_\mathrm{eff}$ and the turbulence decorrelation rate just inside the LCFS. At smaller minor radii, the energy transfer rate remained below $\gamma_\mathrm{eff}$, indicating that the L-H transition is triggered near the LCFS in this case [910]. Measurements on C-Mod yielded similar findings, with the ratio of Reynolds work (which generates low-frequency flow) to the turbulence drive exceeding unity shortly before an L-H transition, and a subsequent drop in the turbulence kinetic energy and an increase in poloidal flow, leading finally to a rapid buildup of an edge pedestal [914].

Zoom In Zoom Out Reset image size