effects of α particles on thermal transport, enhancing under-
standing towards the extrapolation to burning plasmas.

4.2.3. *Theoretical understanding of ion stiffness mitigation.*

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 ($^3$He) – D ICRH [332], and recent experiments on AUG succeeded in producing very peaked $T_i$ profiles optimizing ICRH in (H) – D minority [285]. Global GENE simulations in the inner half of the plasma show a narrow region $\rho_{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 $\beta$, 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_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 ($^3$He) – 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 $\beta_e$ is found when including fast ions with respect to the case of pure thermal $\beta_e$. Two phases for the stabilization process are recognized, as shown in figure 14(b) for the case with highest $\beta_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].

Fast ions are not a necessary ingredient for nonlinear e.m.
stabilization, although they greatly enhance its effect. For the

<figure><img src="image_16.png" /><figcaption>**Figure 13.** GENE nonlinear e.m. simulations of the JET high NBI case of figure 11(b) at $\rho_{tor} = 0.33$, showing the stabilizing effect of rotation and fast ions [316]. Reproduced from [272]. © IOP Publishing Ltd All rights reserved.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 14.** (a) Nonlinear $q_i^{gB}$ from GENE nonlinear e.m. simulations for a case based on parameters of JET high NBI shot with increasing $\beta_e$ values, with and without including NBI fast ions. (b) Time trace of $q_i^{gB}$ and of zonal shearing rate from GENE nonlinear e.m. simulations with fast ions for a case based on parameters of JET high NBI shot, and $\beta_e = 0.012$. The vertical black line marks the time of the transition between phase I and II. Reproduced from [337]. © EURATOM 2019. All rights reserved.</figcaption></figure>

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 $\beta$. 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 $\beta$ stabilization, as discussed in [298]. Also in JET and AUG Hybrid scenarios, the dominant stabilization is by the thermal $\beta$, 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-$\beta$, 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-$\beta$ case. Analytically, it
<figure><img src="image_18.png" /><figcaption>**Figure 15.** (a) $q_i^{gB}$ vs $R/L_{T_i}$ from CGYRO nonlinear e.m. simulations for the DIII-D experiment with high co-injected NBI power at $\rho_{tor} = 0.4$. The e.s. cases ($n_{field} = 1$) are compared against the e.m. cases ($n_{field} = 2$), with and without rotation shear. The experimental values are indicated by the grey bars. Reproduced from [298]. © EURATOM 2019. All rights reserved. (b) Ion heat flux vs $a/L_{T_i}$ from GENE nonlinear e.m. simulations for an AUG high $\beta$ advanced tokamak scenario, comparing the e.s. case to the e.m. cases with and without fast ions. Reproduced courtesy of IAEA. Reproduced from [323]. © 2017 Max-Planck-Institut fur Plasmaphysik. All rights reserved.</figcaption></figure>

was possible to relate the coupling efficiency to a triplet correlation time $\tau$ [345, 346]. While an exact calculation of $\tau$, including stable modes, is complex, a proxy for $\tau$ based on couplings between two unstable ITG modes and a ZF was suggested, which increases with $\beta$. This could be used to account for these effects in quasilinear estimates. In addition, also for the thermal $\beta$ 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].

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_{pFD}$, 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_c}$ and $R/L_{pFD}$, 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.

### 4.2.4. Status of quasilinear modeling of ion heat transport.
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_e/T_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_{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 \times 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, $\alpha$ 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 $\beta$ 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_{T_i}$ 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.

### 4.3. Electron-scale heat transport and multi-scale interactions

Besides the obvious importance of studying and understanding ion heat transport, since ions are the particles involved in fusion reactions, research on electron heat transport is in fact equally important. Reactors will feature dominant electron heating through $\alpha$ heating and, as a consequence, the ions will be mainly heated by collisional heat exchange with electrons. Therefore any electron heat transport mechanism preventing the achievement of peaked $T_e$ profiles will also strongly penalize the $T_i$ peaking and fusion power production. The electron heat transport in present tokamak cores is mainly ascribed to ion-scale TEMs/ITG [354] micro-turbulence and electron-scale ETGs (electron temperature gradient modes) [355]. Both TEM and ETG modes are destabilized with increasing $R/L_{T_e}$, with approximate linear thresholds given by

$$R/L_{T_e}^{\text{crit,TEM}} = (0.357 + 0.271/\sqrt{\epsilon}) [4.90 - 1.31R/L_n + 2.68s + \ln(1 + 20\nu_{\text{eff}})], \quad (7)$$

from [356], where $\epsilon = a/R$ is the inverse aspect ratio, $\nu_{\text{eff}} \simeq n_e [10^{19} \text{ m}^{-3}] Z_{\text{eff}}/T_e$ [keV] the effective collisionality and $Z_{\text{eff}} = \sum_i Z_i^2 n_i/n_e$ (sum on ion species) the effective charge, and

$$R/L_{T_e}^{\text{crit,ETG}} = \max \{(1 + \tau)(1.33 + 1.91s/q)(1 - 1.5\epsilon) \times (1 + 0.3\epsilon d\kappa/d\epsilon), 0.8R/L_n\} \quad (8)$$

from [301], where $\tau = Z_{\text{eff}}T_e/T_i$ and $\kappa$ is the elongation. The corresponding turbulent transport, due to the nonlinear saturation of the micro-instabilities, can in principle set a limit on the $T_e$ peaking. Early studies of electron heat transport with techniques like heat flux scans or $T_e$ modulation were facilitated by the availability of ECH, which guarantees direct electron heat deposition with narrow profiles. As presented e.g. in [291, 357–361], these early studies were made in ECH only plasmas, which were found TEM dominant and with a mild level of electron stiffness. More recently, the same type of experiments were performed in plasmas with NBI+ECH heating, with dominant ITG modes, and here the results were found to be very different, with high electron stiffness, difficult to explain with TEMs alone [296, 361–367]. In this review we will focus on these latest studies, which are more relevant for the reactor situation, always featuring ITG dominance. For cases where the interplay of electron and ion heating results in large values of $R/L_{T_e}$ in conjunction with sufficiently low values of $\tau$, the ETGs could be linearly unstable. However, the nonlinear physics is much more complicated than the linear one and, moreover, the ion-scale ITG/TEM modes can nonlinearly interact with high-k ETGs, as well as with ZFs and ETG streamers (see next paragraph). This could result in cross-scale effects that can manifest in a multi-scale character of the electron heat transport, where a scale separation between low-k and high-k modes cannot be invoked at all. Therefore, the evaluation of the relative impact of ITG/TEM- and ETG-driven turbulence on the electron heat transport is nontrivial,
and has to be investigated in order to be able to predict the
behavior of the electron temperature profiles in present and
future machines.

**4.3.1. Impact of ETGs on the electron heat flux.** 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_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_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_i \lesssim 1$, where $k_\theta$ is the poloidal wave number and $\rho_i$ is the ion Larmor radius), testing if the ion-scale contributions to $q_e$ are sufficient to explain the experiment. If there is a gap between the ion-scale gyrokinetic $q_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_e \lesssim 1$, where $\rho_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 con-
tribute to the experimental $q_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_e$ is underpredicted by ion-scale GYRO
gyrokinetic runs for those cases, while the experimental $q_i$ is
explained, the ETGs are invoked as possibly responsible for
the missing $q_e$. In [361], JET pulses are considered. In figure
14 of [361], the comparison of gyrokinetic GKW nonlinear
ion-scale $q_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} > 9$. These indicates a
possible role for ETGs in setting the $T_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 par-
ticular, 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_\theta$
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 elec-
tron scales and electron scales on ion scales are present. In
figure 16(B) the reduction of the efficiency of ZF shear gen-
eration 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_e$ (green line) with decreasing $a/L_{Ti}$,
with a maximum effect that is obtained when the ITGs are mar-
ginally 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_e$ (yellow line
in figure 16(B)) with respect to the result of ion-scale sim-
ulations (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 exper-
imental 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 previ-
ously observed in [375], where a slightly different mechan-
ism was deduced for ETG quenching by ion scales, with ETG
streamers suppressed by ITG eddies instead of ZFs.

An additional strong dependence of these cross-scale couplings on the $E \times 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 $\gamma_E$ that are smaller than the experimental one, at experimental $\gamma_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 $\gamma_E$ (with increasing $\gamma_E$ for this case) the electron-scale→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 $\gamma_E$ value, the experimental fluxes can be matched with a similar proportion of $q_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 $\tau$ 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
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Ion (A) and electron (B) heat fluxes obtained from scans of the long wavelength turbulence drive, $a/L_{Ti}$, at constant $a/L_{Te}$, using ion-scale (red diamonds) and multi-scale (blue squares) GYRO simulations. Reproduced from [362]. © 2016 IAEA, Vienna. All rights reserved.</figcaption></figure>

using a critical gradient model [376]) on $\tau$ is found consistent with the ETG threshold (equation (8)). The results of the non-linear analysis indicate that the ETGs are needed to explain the experimental $q_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_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_e \sim T_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 \times B$ shearing has only been considered for the TCV cases with NBI (negligible rotation at the radius of analysis for JET and much smaller $\gamma_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_{FI}/n_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_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_e$ is slightly underestimated, while for TCV $q_e$ is strongly underestimated because of a synergic effect of fast ions and $E \times 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_{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_{eff} = 1$, are shown in figure 18. A stabilizing ~20% effect of impurities is found on $q_e$, approximately equally divided between ion-scales ($k_y < 1$) stabilization and electron-scales ($k_y > 1$) stabilization. These results, however, have to be taken with caution, since a relatively small statistics ($\Delta t_{cs}/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.

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
<figure><img src="image_20.png" /><figcaption>**Figure 17.** TGLF stand-alone fluxes, compared with exp. points and GENE ion-scale and multi-scale fluxes, for AUG (a) and JET (b). Except for the results of the new GENE nonlinear ion-scale (red diamond) and multi-scale (red star) simulation taking into accounts impurities. Reproduced from [367]. © EURATOM 2021. All rights reserved.</figcaption></figure>

value of $R/L_{Ti}$ that could allow to match the experimental $q_i$ for smaller $R/L_{Ti}$ than the one needed to match $q_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_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
<figure><img src="image_21.png" /><figcaption>**Figure 18.** Electron heat flux spectra in gB units for the JET GENE nonlinear ion-scale (dotted) and multi-scale (solid) simulations, with $R/L_{Te} = 8$, $Z_{eff} = 1$ (blue), $R/L_{Te} = 11$, $Z_{eff} = 1$ (black) and $R/L_{Te} = 11$, $Z_{eff} = 1.4$ (experimental value: red). The results of the runs with $Z_{eff} = 1$ are adapted from [367].</figcaption></figure>

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_e$ when $T_e \sim T_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_e/T_i \sim 3$. 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, microtearing modes (MTM) can be destabilized with finite $\beta$ 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 \times 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.

### 4.3.2. Simplified criteria to predict the nonlinear ETG relevance based on linear results and reduced models.
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_{ETG}^{\max}/\gamma_{ITG}^{\max}$ of the maximum linear growth rate in the electron-scale and ion-scale k-ranges is larger than $\sqrt{m_i/m_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_e$ when $(\gamma/k_y)^{\max}_{ETG}/(\gamma/k_y)^{\max}_{ITG} > 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 \times B$ shear. The latter criterion reduces to the former in the approximation $k_{ETG}^{\max} \sim \rho_{Te}^{-1}$, $k_{ITG}^{\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_{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_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_{eff} = 1$ assumption, and the impossibility to perform a fine $R/L_{Ti}$ scan down to ITG marginal stability to spot possible electron-scale → 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.

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)^{\max}_{ETG}/(\gamma/k_y)^{\max}_{ITG} > 1$ criterion (ETG are turned on when the criterion is approximately met). TGLF, on the other
<figure><img src="image_22.png" /><figcaption>**Figure 19.** Top: linear benchmark between GKV and GENE for the multiscale parameters for growth rates. Both ion and electron scales are shown with a logarithmic x-axis. The growth rates and frequencies are normalized by the wavenumber $k_y$. Bottom: GENE and GKV electron heat flux spectra in a logx plot. The y-axis is $Q_e \cdot k_y$, such that the areas under the curves in the logx plot correspond to the total heat flux. Less than 10% of the electron heat flux is from electron scales. The heat fluxes are in gyroBohm units with $T_{\text{ref}} = T_e$ and $L_{\text{ref}} = R_{\text{maj}}$. Reproduced from [384]. © 2022 IAEA, Vienna. All rights reserved.</figcaption></figure>

hand, employs a more refined zonal-flow mixing model [33]. However, the inverse electron scale → ion scale effects like the enhancement of ITGs due to ETGs [362, 375](or quenching, depending on $\gamma_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_e$, $q_i$, the percentage electron temperature fluctuation $\delta T_e/T_e$ and the perturbative electron heat diffusivity $\chi_e^{\text{pert}} = -(1/n_e)dq_e/d(\nabla T_e \cdot \hat{r})$, varying $R/L_{Te}$, $R/L_{Ti}$, $R/L_n$ and $Z_{\text{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 $\exp[-((\gamma/k_y)_{\text{TEM}}^{\text{max}}/(\gamma/k_y)_{\text{ITG}}^{\text{max}} - 1)^2/0.09] > 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_e$ and $q_i$, it is possible to explain the experiment by fortuitous agreement, since $\delta T_e/T_e$ and/or $\chi_e^{\text{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_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_e > q_i$ and TEMs that fill a spectral gap between ion and electron scales, while in [299] a regime with $q_e < q_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.

### 4.4. Conclusions and future work

Core ion and electron heat turbulent transport continues to be a focus of research within ITPA, with significant progress in understanding and predicting, due to new experimental results and theoretical developments. In the last decade, two specific topics have been the objects of an intense collaborative work: the effects of electromagnetic stabilization and fast ions on ion heat transport and the role of electron scale turbulence and multi-scale interactions on electron heat transport. The state of the art on both topics is discussed in this paper. Regarding the first topic, a wealth of consistent experimental results in JET, AUG, DIII-D, is now available concerning thermal $\beta$ stabilization and stabilization due to low energy fast ions ($<1$ MeV), so no further experiments are deemed necessary. The understanding of these results seems complete in terms of local gyrokinetic simulations, which for such low energies appear adequate, as has been shown by comparison with global gyrokinetic simulations for selected cases [387]. On the contrary, very little has been done for the effects of high energy fast ions and $\alpha$-particles, which deserve further investigation. Experimentally,
the D-T campaign on JET has produced many results which still require detailed modeling work. More experiments with high energy fast ions, such as those produced by ICRH in the three-ion schemes, would also be valuable. On the theory side, the validity of the local limit in presence of high energy ions is questionable and highly demanding global flux-driven simulations are necessary to assess whether stabilization effects will remain relevant in burning plasmas. Another open point is the upgrade of existing quasilinear models to include these nonlinear EM stabilization effects. This is a work item of highest priority, since present models, used commonly for integrated modeling and projections to future machines, do not include the necessary physics and basically underestimate performance for all scenarios where these effect are important, such as Hybrid and Advanced Scenarios. Regarding the role of electron scale turbulence and multi-scale interactions on electron heat transport, the focus on experiments has been on plasmas with mixed ion and electron heating, since the two conditions for a significant contribution of ETG flux are a ratio $T_e/T_i \le 1$ and high $R/L_{Te}$. Hints of the existence of high electron stiffness and an ETG 'wall' are found in C-Mod, JET, AUG, TCV. Here the measured normalized $q_e$ vs $R/L_{Te}$ cannot be ascribed to ion-scale turbulence only. Multi-scale flux-tube simulations yield a non-negligible electron heat flux due to ETGs when neglecting impurities, but the picture is less clear when impurities are included, particularly on JET, also due to the difficulty of having a large set of heavy multi-scale simulations with sensitivity studies. This remains a point on which further work is needed. Another topic for future work is the improvement and validation of quasilinear models including multi-scale physics. In fact, simplified linear criteria for ETG relevance are not always found in agreement with nonlinear multi-scale simulations, and our confidence in ETG predictions with quasilinear models needs to be increased by further development and validation against GK simulations and experiments. In particular, while the effect of ion scales on electron scales is included in existing models, albeit with simplified treatment, the inverse effect of electron scales on ion scales is presently not included and may be important for integrated modeling. Finally, more studies about the impact of electromagnetic effects on ETGs should be performed, and experimentally it would be important to complement the transport analysis with measurements of density and possibly temperature fluctuations at electron scales, to be compared with synthetic diagnostics applied to multi-scale GK runs.

## Acknowledgments

This work was conducted under the auspices of the ITPA Topical Group on Transport & Confinement. The Authors would like to thank all of the colleagues who have contributed to the ion and electron critical gradient and profile stiffness working group, for the collaborative spirit that led to many interesting experimental and theoretical results and for the valuable discussions during the ITPA T & C meetings.

# 5. Momentum transport

R.M. McDermott, C. Angioni, Y. Camenen, B. Duval, Y. Idomura, J.E. Rice and T. Tala

## 5.1. Introduction

The plasma rotation and its shear are important parameters in fusion devices as they influence the plasma turbulence [388, 389], L-H power threshold [390], both high- and low-Z impurity transport [144, 183, 189, 193, 245], and plasma confinement and transport barriers [391, 392]. They also provide protection against locked modes [393-395] and other deleterious MHD instabilities [396]. The goal of the rotation community is to understand the sources, sinks, and transport of momentum in fusion plasmas in order to provide a predictive capability, as well as to seek ways of optimizing the rotation profile. While some sources and sinks of momentum, such as NBI, are well understood, others, such as ion orbit loss [397], NTV [398, 399] and neutral drag [400] remain active areas of research. In addition, aspects of our theoretical understanding of turbulence driven momentum transport still require experimental validation and a complete theoretical model providing the full turbulent momentum flux within a practical, transport-modeling approach is not yet available. As such, it is not yet possible to regularly include rotation predictions in the modeling and simulations of future devices.

This chapter will focus on the progress made in understanding toroidal, or parallel, momentum transport, as this has been the main focus of research performed within the ITPA TC group. However, NTV and the impact of resonant magnetic perturbations (RMPs) on the rotation are partially covered in section 6. In addition, although significant progress has also been made in the study of poloidal rotation and understanding the differences between main ion and impurity ion rotation, these topics are considered outside the scope of this contribution.

The radial flux of the main ion toroidal angular momentum, $\langle n_i m_i R v_\phi \rangle$, can be written as

$$ \Pi_\phi = n_i m_i R_0 \left( -\chi_\phi R_0 \frac{\partial \Omega_\phi}{\partial r} + V_\phi R_0 \Omega_\phi \right) + \Pi_{RS}. \quad (9) $$

Here $n_i$ and $m_i$ are the density and mass of the main ions, $R_0$ is a reference major radius, $\chi_\phi$ is the diagonal viscosity (diffusion), $\Omega_\phi$ is the toroidal angular velocity, $\Omega_\phi = \langle v_\phi/R \rangle$ with $v_\phi$ the toroidal velocity, $R$ the local major radius, and $\langle \cdot \rangle$ the flux surface average, $V_\phi$ is the convective velocity, and $\Pi_{RS}$ represents the residual stress. Note that, strictly speaking, although the conserved quantity in a tokamak is the toroidal angular momentum summed over all plasma species, see e.g. [401], in practice, considering main ions only is a good approximation provided the impurity concentration is sufficiently low. It should also be noted here that, while the terms in equation (9) are expected to represent the dominant momentum transport mechanisms present in tokamak plasmas, there are other turbulence driven sources of intrinsic rotation that are not included here. For example, turbulent acceleration [402, 403], which
cannot be written in terms of the divergence of a flux, is neglected in these considerations. The first two terms on the right hand side of equation (9) represent diffusive and convective contributions. The final term on the RHS, the residual stress, corresponds to the sum of all the terms that are not proportional to either the angular velocity nor its gradient and can generate a finite flux of momentum in the absence of both rotation and its shear, i.e. the plasma can spin up from rest resulting in what is termed intrinsic rotation.

The following sections present key highlights from the work performed within the rotation community focusing on the ITPA TC activities in the last fifteen years. These works aimed to improve our understanding of momentum transport in fusion plasmas and include advances in theoretical understanding, numerous international joint experiments, a wide variety of new experimental methodologies, and several multi-machine databases that have facilitated cross-machine comparison. In 2007, the ‘Plasma confinement and transport’ section of the IPB [3] featured a collection of experimental observations and concluded that our understanding of momentum transport is insufficient to the task of predicting the rotation profile. In comparison, today, we have a considerably more robust and well developed theoretical basis that includes the physical understanding of the origins of the turbulent convection as well as the identification of key residual stress generation mechanisms. Moreover, we have global simulations capable of reproducing intrinsic rotation observations and advanced experimental techniques capable of validating theoretical predictions. These advances are presented in the sections that follow.

This chapter is organized as follows. In section 5.2, an overview of the development of momentum transport theory is given. This section is divided into a discussion of progress made in understanding the diffusive 5.2.2, convective 5.2.3, and residual stress 5.2.4 momentum fluxes, together with a section dedicated to advances in simulating these fluxes 5.2.5. Section 5.3 discusses experimental works that measure some of the components of the momentum flux, with a focus on those that have been used to validate theoretical predictions. This section presents experiments aimed at identifying diffusive versus non-diffusive transport in sections 5.3.1–5.3.3 and those aimed at obtaining information on the residual stress in section 5.3.4. Finally, a summary of the progress together with a discussion of open areas of research is presented in section 5.4.

## 5.2. Development of theoretical understanding

From early work on turbulent momentum transport, [404] ion heat and momentum transport were expected to be strongly coupled. As a result, in the absence of a strong external momentum drive, such as NBI, the magnitude of the core rotation would be negligibly small due to strong turbulent momentum diffusion. This presented a potential problem for future, large machines, such as ITER, in which the applied beam torque will be small compared to the plasma inertia, such that any remaining rotation might be insufficient to provide the required stability. However, experimental measurements

began to show evidence for an inward momentum convection [405] and significant rotation in the absence of external momentum drive [406, 407]. These measurements demonstrated an unexpected complexity within the momentum transport channel and provided hope that alternative transport mechanisms could be strong enough to compete with turbulent diffusion, and result in significant rotation, even in large, non-NBI dominated machines. Intrinsic rotation was, in fact, widely observed and inter-machine scaling studies suggested significant self-generated flow was possible without externally applied torque [408]. Recent reviews of experimental rotation observations include: [409–414].

Experimental observations motivated a re-evaluation of the theoretical picture of momentum transport that resulted in the identification of the Coriolis momentum pinch as well as a variety of residual stress mechanisms. This section provides an overview of the progress that has been made on the theoretical understanding of turbulent momentum transport. As such, it does not intend to be comprehensive, but rather informative. For a more complete review, there are several papers that can be consulted [27, 415–418]. In addition, there are several works that consider the electromagnetic contribution to momentum transport and its role in driving intrinsic rotation. These are not considered specifically in this contribution, but could be important. More information on these can be found in [419–421].

### 5.2.1. Ordering.
The local delta-f approximation of the gyrokinetic equation retains only the lowest order terms in a $\rho_*$ expansion ($\rho_* = m_i V_{th}/L_{\perp} e B_t$), where $V_{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_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 \times 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:

$$ \frac{\Pi_{\phi}}{n_{i}m_{i}R_{0}} = -\chi_{\parallel}v_{\parallel}^{\prime} - s_{j}\chi_{E}\gamma_{E} + R_{0}V_{\phi}\Omega_{\Phi} + C_{FS} + \frac{\Pi_{*}}{n_{i}m_{i}R_{0}} \quad (10) $$

where $v_{\parallel}^{\prime}$ is the flux-surface-averaged parallel flow shear, $\gamma_E = \frac{1}{B} \frac{\partial^2 \Phi}{\partial r^2}$ is the $E \times 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 \times B$ shearing rate are coupled and proportional to the toroidal flow shear:
$v'_\parallel \sim \frac{RB_i}{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_i}{B} - \chi_E \frac{B_p}{B}$. The next two contributions are the Coriolis pinch $R_0 V_\Phi$ and the residual stress driven by any up- down asymmetry in the magnetic equilibrium, $C_{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.

**5.2.2. Diffusive flux.** The first work that showed that the ion heat diffusivity, $\chi_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, $\text{Pr} = \chi_\Phi/\chi_i \sim \chi_\parallel/\chi_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_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_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 gyrokinetic simulations of ITG turbulence with kinetic electrons to show that $\chi_\parallel/\chi_i$ varies between only 0.7 and 1.2 over a wide range of plasma parameters. In addition, $\chi_\parallel/\chi_i$ was predicted to depend only weakly upon plasma $\beta_e$, magnetic shear, and $R/L_{T_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_i > 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 \times 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 > 0$, making the total momentum diffusivity, $\chi_\Phi = \chi_\parallel \frac{B_i}{B} - \chi_E \frac{B_p}{B}$, lower and thereby reducing the Prandtl number. For small toroidal rotation shear, typical of intrinsic rotation, a Prandtl number as low as 0.5 can be obtained when including the $\chi_E$ contribution [433]. For negative magnetic shear, the situation reverses with $\chi_E < 0$ [433]. In contrast to $\chi_\parallel/\chi_i$, which is mostly independent of $v'_\parallel$, $\chi_E/\chi_i$ decreases with increasing $\gamma_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.

**5.2.3. Convective flux.** 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 com-
ponent perpendicular to the toroidal direction and, therefore,
feels a Coriolis force, $\vec{F}_C = 2m_i v_{\parallel} \mathbf{b} \times \mathbf{\Omega}_{\phi}$ with **b** the unit vec-
tor 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} > 0$ (or reduces it for
$v_{\parallel} < 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} > 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} < 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 coup-
ling 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_{ne}$. 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 elec-
trons, 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$, as well as dependences on the magnetic shear, safety
factor, and ITG [437] as shown in figure 20.

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 pre-
dicted under most conditions in the presence of TEM turbu-
lence. For both ITG and TEM turbulence, in standard aspect
ratio tokamaks at moderate plasma $\beta$, the resulting convection
is inward, meaning it acts to strengthen the rotation. However,
for high-$\beta$ and low aspect ratio, as found in spherical tokamak
plasmas, sign reversals of the convection have been obtained
within gyrokinetic simulations [442, 443].

<figure><img src="image_23.png" /><figcaption>**Figure 20.** The parameter dependences of the Coriolis pinch number were explored via linear gyrokinetic simulations in [437]. The dependence as a function of $R/L_{ne}$ is shown in red crosses, the shear 3s as blue stars, the safety factor $q$ as green circles, $R/L_{Ti}$-6 as magenta diamonds, and $8k_{\theta}\rho_i$ as black squares. Reprinted figure with permission from [437], Copyright (2007) by the American Physical Society.</figcaption></figure>

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.

**5.2.4. Residual stress.** 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 compar-
ison 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_{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 $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_{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 \times B$ flow [453, 454]. The relationship between the toroidal flow shear, $\frac{\partial \Omega_{\phi}}{\partial r}$, and the $E \times B$ flow shear $\gamma_E$ is given by [119, 455]:

$$ \gamma_E = -s_j \frac{1}{|\nabla r|} \frac{RB_p}{B} \frac{\partial \Omega_{\phi}}{\partial r} + \gamma_E^* \quad (11) $$

with

$$ \gamma_E^* = s_j \frac{1}{|\nabla r|} \frac{RB_p}{B} \frac{\partial}{\partial r} \left\langle \dot{u}_{\theta} \frac{B_t}{R} + \frac{1}{Z_s en_s} \frac{\partial r}{\partial \Psi} \frac{\partial p_s}{\partial r} \right\rangle. \quad (12) $$

Here, $\dot{u}_{\theta} = (v \cdot \nabla\theta)/(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 $\dot{u}_{\theta}$, one obtains:

$$ \gamma_E^* \sim -\rho_* v_{\text{thi}} \frac{\partial}{\partial r} (k_i R / L_{Ti}) + \rho_* v_{\text{thi}} \frac{\partial}{\partial r} R / L_{pi} \quad (13) $$

where $v_{\text{thi}} = \sqrt{T_i/m_i}$ is the ion thermal velocity, $k_i$ the neoclassical flow coefficient ($k_i \sim 1$ in the banana regime), $R/L_{Ti}$ the normalized ITG and $R/L_{pi}$ 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 $\gamma_E$ in the presence of steep profiles such as those encountered in the edge pedestal. They enter the parallel and $E \times B$ flow shear contributions to the momentum flux as:

$$ \left( \frac{\Pi_{\phi}}{n_i m_i R_0} \right)_{\parallel+E \times B} = - \left( \chi_{\parallel} \frac{B_t}{B} - \chi_E \frac{B_p}{B} \right) \frac{\partial \Omega_{\phi}}{\partial r} - s_j \chi_E \gamma_E^*. \quad (14) $$

As discussed in section 5.2.2, for positive magnetic shear, the diffusive flux driven by $E \times B$ flow shear is positive, $\chi_E > 0$. In this case, and with $\gamma_E^* > 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 \times B$ residual stress often refers to the momentum flux related to $\gamma_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 \times B$ momentum flux, related to $\gamma_E$ and the $E \times 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 \times B$ shear is closely linked to the tilting of the turbulent eddies away from the low-field side mid-plane ($\theta_0 = 0$). The resulting stress is proportional to the tilt angle, $\theta_0$, which is determined by the amplitude and sign of $\gamma_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 \times B$ shear was developed and proved to accurately capture the $E \times B$ momentum flux in quasilinear models (together with the turbulence stabilization by $E \times B$ shear).

Together with the $E \times 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 $\theta_0$. 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 \times B$
residual stress and the resulting momentum flux can be formally written as

$$ \left( \frac{\Pi_{\phi}}{n_i m_i R_0} \right)_{2^{\text{nd}}\text{der.}} = -s_j \chi_n \gamma_n - s_j \chi_T \gamma_T \quad (15) $$

with $\gamma_n \sim \rho_* v_{\text{thi}} \frac{\partial}{\partial r} R/L_n$ and $\gamma_T \sim \rho_* v_{\text{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 > 0$ and $\chi_T \sim -\chi_E$ is observed, meaning that positive $\gamma_n$, as typically encountered in the edge region, will drive an inward flux of co-current momentum whereas positive $\gamma_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 > 0$ was still observed, but $\chi_T \sim 0$. Similarly to the $E \times B$ residual stress, there is a saturation of the driven flux due to $\chi_n$ and $\chi_T$ decreasing towards zero at high values of $\gamma_n$ and $\gamma_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.

### 5.2.5. Momentum transport modeling.
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 gyrokinetics 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 \times 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 \times 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

$$ u' = - \frac{\Gamma_{\text{RS}}}{\chi_{\phi}} - \frac{R_0 V_{\phi} u}{\chi_{\phi}}, \quad (16) $$

where a normalization to the ion thermal velocity, $v_{\text{th,i}} = \sqrt{2T_i/m_i}$, has been introduced and the equation written in terms of the normalized rotation gradient, $u' = (-R_0^2/v_{\text{th,i}}) d\Omega_{\phi}/dr$, normalized rotation $u = \Omega_{\phi}R_0/V_{\text{th,i}}$, and normalized residual stress, $\Gamma_{\text{RS}} = \Pi_{\text{RS}}/(n_i m_i v_{\text{th,i}})$. 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, $\theta_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 $\theta_0$ is possible as

$$ u'_{p} = - \frac{\Gamma_{\text{RS}} \theta_0}{\chi_{\phi} \theta_i} - \frac{R_0 V_{\phi} u}{\chi_{\phi}}. \quad (17) $$
Hence, when comparing quasilinear results to experiment, the value of $\theta_0$ 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 \times B$ flow to explain hollow rotation profiles, and a mechanism dependent on the sign of the mode frequency, to explain co-$I_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 \times B$ flow as the dominant driver of hollow rotation profiles, while the co-$I_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

<figure><img src="image_24.png" /><figcaption>**Figure 21.** Global GTS simulations succeed in quantitatively reproducing measured intrinsic rotation gradients observed in DIII-D. Reprinted figure with permission from [472], Copyright (2017) by the American Physical Society.</figcaption></figure>

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_{th,i}$ were obtained.

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.

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
<figure><img src="image_25.png" /><figcaption>**Figure 22.** A strong correlation between the normalized flow gradient, $u'$ and the second derivative of the density profile was found within global simulations of a database of L-mode intrinsic rotation measurements at AUG. Reproduced from [466]. © 2018 IAEA, Vienna. All rights reserved.</figcaption></figure>

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 \times 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 \times 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 \times 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 \times 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 gyrokinetic 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, $\gamma_n$ and $\gamma_T$, at low values of shearing and then a saturation when the level of shearing increases. Since $\gamma_n$ and $\gamma_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 \mathcal{L}_\phi + \partial_\psi (\Pi_{\phi,turb} + \Pi_{\phi,neo}) + \partial_\psi T_\phi = \mathcal{J}$, where $\mathcal{L}_\phi$ is the gyrocenter angular momentum, $\Pi_{\phi,turb}$ and $\Pi_{\phi,neo}$ are turbulent and neoclassical components of the Reynolds stress, $T_\phi$ is the field term stress, and $\mathcal{J}$ is the $J \times 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 $\mathcal{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 $\mathcal{O}(\rho_*^3)$ under the low flow ordering, $v_\parallel/V_{th} \sim \mathcal{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 $\mathcal{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_* \le 220$ showed that the profile shear stress scales as $\Pi_*/(n_i m_i R_0 \chi_\phi) \sim v_\parallel^3 \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_i/L_B \sim \rho_i/R$ and the third order gyro-center transform in $\epsilon \sim \rho_i/L_n \sim \rho_i/a$, where $\rho_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 $\mathcal{O}(\epsilon_B^2) \sim \mathcal{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 \times B$ torque is negligible in the plasma core. However, in the presence of transient particle transport, e.g. in modulation heating experiments, the $J \times 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 \times 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 $\theta_0$,
which depends on the profile shear and the $E \times B$ flow [489]. According to the higher order ballooning theory [492], $\theta_0$ is determined by the sum of the profile shear and the $E \times B$ shear. In addition, global simulations show that $\theta_0$ is affected not only by the $E \times B$ shear, but also by the magnitude of $E \times B$ flow itself [489, 493]. Although the profile shear stress theory [457] was proposed without taking account of the $E \times B$ flow, it enables a unified treatment of these effects through the tilt angle $\theta_0$.

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_{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 \times B$ shear. By comparing simulations with and without an adaptive momentum source, it was shown that the intrinsic rotation enhances the $E \times 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 \times 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 \times 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 \times 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.

### 5.3. Experimental work and theory validation

The starting point for the momentum transport studies within the ITPA TC group was the observation that the experimentally determined effective Prandtl number, $Pr_{eff} = \chi_{\phi,eff}/\chi_{i,eff}$, was well below one in JET, typically $Pr_{eff} = 0.2-0.3$ [504], in seeming contradiction to theoretical predictions [356], see section 5.2.2. However, $Pr_{eff}$, by definition, assumes the momentum flux can be modeled as a purely diffusive term, (see equation (18)),

$$ \Pi_{\phi} = n_{i} m_{i} R \left( - \chi_{\phi,eff} R \frac{\partial \Omega_{\phi}}{\partial r} \right). \qquad (18) $$

It, therefore, also includes all non-diagonal contributions such as the Coriolis momentum pinch and residual stress. It is in the absence of these terms, that one expects $Pr_{eff} \sim 1$. Measurements showing $Pr > Pr_{eff}$ were, therefore, indirect experimental evidence of a non-negligible, net-inward-directed, off-diagonal, momentum flux comparable in magnitude to the diffusive flux. Similar conclusions were reached
by an analysis of a large database of JET discharges [505, 506]. Consequently, the motivation to search for this flux in experiments was strong and became realizable in several devices equipped with improved diagnostic systems, using specific experimental techniques and new modeling capabilities. Section 5.3.1 provides an overview of the various ITPA TC joint experiments carried out to separately identify diffusive and non-diffusive momentum fluxes in multiple tokamak devices.

**5.3.1. Experimental separation of diffusive and non-diffusive fluxes.** 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 Prandtl 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.

**5.3.2. The diffusive momentum flux.** Results from experiments at JET [522], which explored the momentum transport as a function of $q, R/L_{ne}$, and $v^*$, 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.

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 $\beta_N$. Similar to the JET results and theoretical expectations, no dependences with $\rho_*$ or $v^*$ 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,
<figure><img src="image_26.png" /><figcaption>**Figure 23.** Pinch and Prandtl number derived from an NBI modulation experiment of a D plasma in ASDEX Upgrade (green lines) compared to theoretical predictions from the GKW code (symbols) for D (brown crosses) and substituting in the H mass ratio in the calculations (purple circles). Reproduced from [521]. © 2023 The Author(s). Published by IOP Publishing Ltd on behalf of the IAEA.CC BY 4.0.</figcaption></figure>

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 $\tau_{\phi}$ 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 $\tau_{\phi}$ 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_{\phi}^{\text{inc}}$ 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 $\tau_{\phi}$. 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 $\tau_{\phi}$, 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.

5.3.3. The convective momentum flux. 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 $\chi_{\phi}$ [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_{95}$ 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

<figure><img src="image_27.png" /><figcaption>**Figure 24.** A multi-machine database of core momentum transport coefficients was assembled in [525] and demonstrated an $R/L_{ne}$ dependence of the pinch number at fixed $q_{95}$ using data from multiple devices. Reproduced from [525]. © 2012 IAEA, Vienna. All rights reserved.</figcaption></figure>

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.

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, $\epsilon$, 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 experi-
ments 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 agree-
ment 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, predict-
ing the value of the rotation at the plasma edge, or pedestal
top, becomes critical to determining the direction and mag-
nitude 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].

**5.3.4. Residual stress.** Evidence for the role of Π<sub>RS</sub> 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<sub>p</sub> 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<sub>p</sub> 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<sub>p</sub> intrinsic rotation could be balanced, enabling the identification and quantification of a co-I<sub>p</sub> intrinsic torque density in the edge of those plasmas. In contrast, in the presence of strong electron heating, counter-I<sub>p</sub> 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 gener-
ation 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 ρ* scaling of the residual stress and the interaction with the Coriolis momentum pinch. Since residual stress terms are predicted to scale with ρ* (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 ρ* ~ √Tᵢ/BL⊥. In the core the plasma minor radius, a, is a good approximation for L⊥. 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. ρ*<sub>,edge</sub> ~ √Tᵢ/(BL⊥<sub>,edge</sub>). 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 ρ*<sub>,edge</sub> 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_{*,\text{edge}}, v_{E \times B}, \dots$), 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, poten-
tially contributing to a boundary condition, but it is not expec-
ted to be important in the confinement plasma region. NTV
torques and the impact of RMPs on plasma rotation are cer-
tainly 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 incor-
porated 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.

<figure><img src="image_28.png" /><figcaption>**Figure 25.** The Alfven Mach number as a function of $\rho_*$ for a database of intrinsic rotation data collected from six tokamak devices. Reproduced from [408]. © 2007 IAEA, Vienna. All rights reserved.</figcaption></figure>

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 ($\Delta P$), the magnetic field ($B_t$), and inversely with the plasma current ($I_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_{th}/I_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$^{-1}$. 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_A \sim \beta_N \rho_*$, where $M_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$^{-1}$. 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 \times 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 \times 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 innershear layer should be $\text{co-}I_p$ directed. In the pedestal region, the contribution of the toroidal flow to the $E \times B$ shear ($\gamma_E$)
is small and $\gamma_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^* < 0$ and the corresponding momentum flux drives counter-$I_p$ directed rotation. In contrast, in the inner $E_r$ shear layer, the $E \times B$ residual stress drives co-current rotation, consistent with the co-$I_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 \times B$ shear based residual stress model [451] and yield moderate core rotation (Mach numbers of $M_1 = 0.07, V_\phi \sim 50-75 \text{ km s}^{-1}$) predictions for ITER and much stronger rotation in ARC ($M_1 = 0.2, V_\phi \sim 200 \text{ km s}^{-1}$). 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_{th}/I_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_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_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_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_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_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, $\beta$, and $q_{95}$ and resulted in a similar level of total intrinsic torque for ITER ($\sim$20 Nm) and lower levels predicted for SPARC and ARC ($\sim$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^2$ (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-Ip 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 \times 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.
The previous section discussed scalings of the intrinsic rotation with a focus on the plasma edge, which is expected to be more relevant in future machines than rotation in the core. However, the physics of the intrinsic generation in the plasma core is also relevant for the edge and, in many ways, provides a cleaner backdrop for the study of residual stress. Whereas in the plasma edge multiple mechanisms are at play, in the plasma core, changes in intrinsic rotation can often be directly attributed to residual stress alone.

A wide range of observations of core localized residual stress are available. Measurements by Solomon et al [497] demonstrated the existence of a core-localized, co-Ip directed intrinsic torque in NBI heated plasmas using balanced beam experiments as well as a strong counter-current torque in an ECRH heated plasma. In those experiments, no scaling of the intrinsic torque with the local pressure gradient was observed in the core, in contrast to the observations at the plasma edge, but qualitatively consistent with the scaling of the $E \times B$ shear across the two regions. The observation of strong counter-current residual stress torque with ECRH is consistent with older DIII-D observations of rotation changes in the counter-current direction, when the plasma was heated with ECRH power [538]. Similar observations were made in other machines including JT-60U [519], KSTAR [560] and AUG [140, 234]. The magnitude of the intrinsic torque seen in the core of tokamak plasmas with direct electron heating is comparable to that typically seen in the edge pedestal, but, as they are often in the opposing directions, if present and non-negligible in future machines, they could act to reduce the core rotation. As such it is of interest to understand the residual stress drive and quantify its impact on the rotation in future machines. It is worthwhile to note, however, that while rotation is generally observed to decrease when RF power is applied to NBI heated discharges, this decrease is not necessarily due to an intrinsic torque. Rather, adding RF heating generally also increases the turbulent viscosity, leading to flattening of both the ion temperature and toroidal rotation profiles [561].

The drive mechanisms generating the residual stress can be probed by studying the parameter dependences of the observed intrinsic rotation and by comparison of these to gyrokinetic simulations. Intrinsic rotation experiments in the absence of
ion heating, and hence strong gradients in $T_i$ at the edge, are particularly interesting for this purpose. ASDEX Upgrade Ohmic L-modes, ECRH heated L- and H-modes, and ICRF heated L- and H-mode plasmas are dominantly electron heated and feature only very weak ion temperature pedestals. Under these conditions the edge rotation is only a few km s$^{-1}$ in the co-$I_p$ direction, resulting in a momentum pinch that remains approximately constant and, almost, negligibly small. Hence, the core is effectively decoupled from the edge and momentum transport becomes a balance between diffusion and the residual stress. In steady state, the normalized rotation gradient, $u'$, then gives an indirect measurement of the residual stress, normalized by the momentum diffusivity, see equation (16).

In [465] a database of intrinsic $u'$ measurements, from a variety of confinement regimes and plasma conditions, was examined and the simple observation was made that the core rotation magnitude is directly correlated to the rotation gradient at mid-radius. Put another way, the plasma turbulence at mid-radius builds a positive, or negative, rotation gradient setting the profile shape, and this results in a higher or lower central rotation. Thus, the absolute sign and magnitude of the central rotation can be useful as a diagnostic of the rotation gradient at larger major radii, but do not inherently contain any information on the underlying transport driving that rotation. Rather, it is the rotation gradient at the location of the maximum residual stress, i.e. the location of maximum rotation gradient, that is important. This was also reported in Camenen et al [450], which examined the impact of the magnetic configuration on the edge intrinsic rotation gradi-ent. That work showed that the absolute value of the rotation at the edge can be matched across different magnetic configurations, but the residual stress acts upon the gradient, which then determines both the shape of the resultant profile and the magnitude and sign of the core rotation. The AUG database also demonstrates that there are no fundamental differences between the core localized intrinsic rotation in L-mode and H-mode, once decoupled from the plasma edge [465]. In all cases, the intrinsic rotation profiles and their gradients were observed to behave identically as a function of local plasma parameters, irrespective of heating method or confinement regime, with both peaked and hollow intrinsic rotation profiles being observed in both the L- and H-mode plasmas.

Previous observations proposing more fundamental differences between L- and H-mode plasmas [408, 475] examined cases in which edge effects were more dominant. In particular, in L-mode plasmas, magnetic geometry can strongly impact the edge rotation and rotation gradient, which then influences the rotation at smaller radii as already discussed [450]. In contrast, in H-mode plasmas with strong pedestals the edge intrinsic rotation is typically non-negligible. Together with the pinch, this strongly impacts the core profile. Thus, the core rotation behavior in L- and H-mode can appear fundamentally different, due to the influence of different edge effects, while the physics of the core-localized residual stress in both confinement regimes, in fact, remains the same.

The AUG database contained both positive and negative $u'$ values and quasilinear gyrokinetic simulations were performed over the entire data-set at mid-radius to determine

the frequency of the dominant turbulent mode, $\omega_r$. The database did not show a clear correlation between $u'$ and $\omega_r$. [140, 465], inconsistent with a residual stress mechanism that changes sign with $\omega_r$ being the (only) dominant mechanism. Instead, the database showed a linear correlation with $R/L_{ne}$. Peaked rotation profiles are correlated with flat electron density profiles, and hollow rotation profiles with peaked ones. However, significant scatter exists. While $R/L_{ne}$ was clearly identified as the strongest correlating parameter, strong correlations with $R/L_{T_i}$ and collisionality were also observed [140], demonstrating the complex multi-dimensional aspects of the phenomenon. A strong correlation between $R/L_{ne}$ and $u'$ was also recently identified limited TCV plasmas [562] and was examined within the C-Mod intrinsic rotation database [475, 476, 563, 564]. There, reversals were identified with no changes in the electron density profile, also demonstrating the presence of other equally important parameters contributing to the changes in residual stress. This is suggestive of a bifurcation, which is supported by the temporal evolution of the profile in C-Mod and TCV at the reversal, see discussion in [414].

These dependences were further explored in the context of the 'rotation reversals' that are observed in Linear to Saturated Ohmic Confinement (LOC-SOC) transitions in tokamak plasmas. This phenomenology is described in detail in many works, see [565] and references therein, and is also discussed in section 9 of this chapter. For diverted plasmas in the LOC regime, at low plasma density, or more correctly low collisionality as collisionality has been shown to better unify the experimental observations [38, 566], $u'$ is generally positive around mid-radius leading to peaked or flat core rotation profiles and, therefore, to co-$I_p$ rotation. As the density and collisionality are increased, the plasma transitions toward the SOC regime, $u'$ becomes negative leading to hollow rotation profiles and often to counter-$I_p$ rotation in the plasma center. It should be noted that this reversal is observed in the opposite direction in limited plasmas [414, 562], a phenomenon which is not yet understood. Historically, the term 'rotation reversal' referred to the sign change of the central rotation value that is often observed under these conditions. However, the name is not apt, as it is possible to have a change in the sign of the rotation gradient, without a change in sign of the core rotation, as was demonstrated, for example, in KSTAR [534]. In addition, it is possible to have a change in the sign of the rotation without a change in the sign of the gradient, as discussed in [562]. Hence, the identification of a reversal should always be made with respect to the sign of the gradient. For diverted plasmas, at even higher collisionalities a second reversal from negative $u'$ back to positive has been observed in AUG and KSTAR [38, 465, 534].

A transition from TEM turbulence to ITG driven turbulence is often predicted in the same parameter range in which the LOC-SOC and the first (positive to negative $u'$ with increasing collisionality) rotation reversal takes place. Some devices have presented evidence for such a change in turbulence coincident with the change in rotation [567, 568], while other devices observe a change in the turbulence propagation direction only after the rotation gradient has switched sign [38, 562, 569]. As
the tilt angle of the turbulent eddies changes sign from TEM to ITG, residual stress mechanisms that are proportional to the angle have been considered to explain the observed rotation reversal. However, this picture is oversimplified for a number of reasons.

The change of dominant turbulence regime in these plasmas is not a switch, but rather an evolving process characterized by strong, co-existing, sub-dominant modes that spreads in both space and time [467, 562, 568]. Moreover, while such a mechanism is potentially able to qualitatively explain the first reversal, which takes place near a change in dominant turbulent mode, it cannot explain any reversal that takes place in ITG dominated conditions. Rotation reversals observed in ICRF heated plasmas in ITG dominated conditions [563] can also not be explained by a change in the underlying turbulence type. However, these observations do not necessarily rule out profile shearing or intensity gradient as relevant mechanisms, first because the background turbulence is a mix of evolving ITGs and TEMs as the plasma changes and second, even if the dominant turbulence type were clear, the sign dependences of these residual stress terms is more complicated than just the relationship to the turbulent propagation direction, see section 5.2.4.

As such, global gyrokinetic modeling is required to study the total residual stress produced for different experimental conditions and the mechanisms by which they are produced. It would be extremely interesting to apply this analysis to rotation reversals seen in ITG dominated plasmas such as the negative to positive $u'$ switch observed in AUG at in even higher collisionality L-modes or the reversals seen in ICRF heated plasmas in Alcator C-Mod [563]. Similarly, the AUG database was modeled via the global gyrokinetic simulations with GKW [466]. This work demonstrated that global simulations could reproduce the hollow experimental observations, similarly to the DIII-D analyses [472, 477]. Additional modeling of the AUG data-set, particularly the co-$I_p$ directed rotation profiles at low and high collisionality could provide interesting comparison points to the DIII-D work and insight into dominant momentum transport mechanisms in the core of ITG dominated plasmas. Moreover, rotation reversals with extremely similar background profiles [475, 476, 563, 564] make for interesting targets for future global, modeling efforts. To aid such efforts, an ITPA profile database of intrinsic rotation measurements from AUG, C-Mod, TCV, and JET was formed and is freely available to the modeling community.

It should be noted that plasma turbulence is not necessarily responsible for the whole variety of intrinsic rotation and toroidal rotation reversals observations. For example, NTV has been proposed to explain some experimental features of the reversals [570–572]. The NTV mechanism requires the presence of 3D field structure in the plasma core, such as internal kink modes, which results in a core-localized NTV torque, the direction of which depends on the plasma collisionality. At low collisionality, the ‘electron root’ solution to the ambipolarity equation results in co-current directed intrinsic rotation, while at higher collisionality there is a bifurcation to the ‘ion root’ solution, which causes a switch to the counter-$I_p$

direction. This change of torque direction could contribute to co-$I_p$ to counter-$I_p$ reversals that occur with increasing electron density, as happens at low-density in LOC to SOC transitions. The second transition back to the co-$I_p$ direction, however, at even higher density is not compatible with an NTV explanation.

### 5.4. Summary and future work

Excellent progress in understanding momentum transport has been made in recent decades and much of this work was performed within the auspices of the TC ITPA group. Our theoretical understanding has expanded from a purely diffusive picture to the identification of a sizable, inward convection capable of competing with turbulent diffusion to provide significant peaking of the core rotation profile in tokamak plasmas. In addition, the presence of large, residual stress momentum fluxes were demonstrated conclusively in experiments, motivating the community to develop the accompanying theoretical description and modeling capability.

The residual stress mechanisms that are believed to be dominant in present day experiments include the stress arising from the up-down asymmetry in the magnetic equilibrium, which is the only residual stress term to lowest order in $\rho_*$, the residual stress arising from first order contributions in $\rho_*$ to the flow shear ($E \times B$) and, finally, residual stress arising from radial derivatives and second derivatives of the background plasma profiles. The implementation of these turbulent momentum transport mechanisms in gyrokinetic codes have greatly improved the predictive capability of the observed rotation behavior in present day experiments and laid the foundation for theory based prediction for future devices.

In particular, the gyrokinetic models for the Prandtl and pinch numbers are believed to be quite robust. Experimentally, Pr has been shown to exhibit only limited variation over a wide range of plasma parameters, similar to theoretical expectations. The pinch was shown to exist and to be inward under most conditions and was also observed to increase with increasing $R/L_{ne}$, consistent with theoretical predictions. Although the number of works that provide quantitative comparison to gyrokinetic predictions is limited, those that have shown agreement, generally, to a factor of two. Moreover, recent works that include the residual stress and time dependence within the analysis are starting to provide even more stringent validation of the theoretical predictions, achieving agreement between theory and experiment within the experimental uncertainties. These results were obtained in ITG dominated plasmas. Validation of the theoretical predictions for TEM dominant conditions or in mixed-mode conditions, as is expected to be the case in future machines, is not yet available. In addition, the inter-machine database revealed a correlation between Pr and machine size, which has not been examined through corresponding gyrokinetic simulations. Understanding this trend is important, as larger Pr numbers in bigger machines would lead to significantly higher momentum diffusivity, and flatter rotation profiles, than presently expected.
Great advances have also been made in modeling the residual stress driven intrinsic rotation profiles. It was shown that physics contained within local models is insufficient to reproduce the experimental *u'* values and, instead, a global approach is required. Experimental intrinsic rotation profiles have been quantitatively reproduced via a variety of methods including Eulerian Vlasov codes like GKW and by integrating the residual stress momentum flux produced by global particle in cell codes like GTC. While these results are an important step forward in our predictive modeling capability, it is not clear how applicable such methods are to future machines as these codes demonstrate strong sensitives to small changes in background plasma conditions. Moreover, the experimental validation of such sensitivities present an extreme challenge to experiment, as it requires measurements of not only first, but second derivatives of the background kinetic profiles. The global simulations also shed light on the dominant mechanisms responsible for the intrinsic rotation seen in the experiment. However, to extract this information diagnostics must be implemented within the codes and different approaches have been employed within the literature, making it difficult to compare results. However, no strong inconsistencies presently exist and zonal $E \times B$ flows seem a clear contender to generate hollow rotation profiles. The implementation of the same diagnostics in the different codes would allow more stringent, and conclusive, assessments of the relative roles of the different mechanisms.

Altogether, the modeling of momentum transport in the plasma core using quasilinear GK models that are fast and inexpensive to run is quite good, so long as the residual stress is small compared to the other momentum fluxes. Under this condition, the resultant predictions for the core rotation profile are reasonable provided an appropriate edge boundary condition is given. Looking towards modeling of the rotation profile in future devices, it remains unclear if the core residual stress will be large enough to significantly impact the rotation profile. This requires additional work on the scaling of residual stress with $\rho_*$ likely via a multi-machine approach. However, for the same background core gradients and curvatures, it is reasonable to expect that the core produced residual stress will be smaller than in present day devices, making predictions that rely only on diffusion and convection in the core appropriate. The edge boundary condition remains crucial and is, by far, the largest unknown.

The pedestal generated residual stress may still be significant in future devices on account of the strong radial gradients in this region. Unfortunately, the modeling capability for residual stress in the edge pedestal is still lacking as are modeling frameworks that integrate all effects known to be important in this region. In particular, the impact of magnetic perturbations on the edge rotation may be particularly important to include, given that RMP ELM-suppression is envisaged as the baseline operating scenario in ITER. Lacking this kind of integrated first principle modeling, estimates of the rotation profile in future machines still rely on empirical scalings of the pedestal generated intrinsic torque. A significant amount of effort has been invested in this topic via a variety of techniques and have resulted in predictions for ITER ranging from 2 to

33 N·m. While the former is negligible compared to the NBI torque, the latter would be sufficient to have a positive influence on the edge rotation and positively influence the core rotation and rotation gradient. The uncertainty on these scalings, however, emphasizes the need for continued work. Pedestal intrinsic torque measurements from larger machines, such as JT60-SA would be particularly interesting for constraining empirical scalings. On the theoretical side, what is missing are reduced models and integrated modeling frameworks to include all effects believed to be important in the edge and to use these to perform quantitative comparisons against experiment. Another element that is presently missing from our modeling capability is the connection of the core flows with those in the SOL. This is an important element for developing reliable models for the rotation at the boundary, which can impact the $E \times B$ shear in this region with implications for pedestal transport and stability as well as the L-H transition.

## Acknowledgments

This work was conducted under the auspices of the ITPA Topical Group on Transport & Confinement. The Authors would like to thank all of the colleagues who contributed to improving our understanding of momentum transport over the years. In addition, the authors are very grateful for all of the contributions to this work. In particular, the authors appreciate the discussions with and support from W. Solomon, C. Chyrstal, and B. Grierson.

# 6. Impact of 3D fields on transport

M. Jakubowski, S. Mordijck, E.J. Paul and A. Zocco

## 6.1. Introduction

In the last two decades there has been an acceleration in the study of the application of small 3D modifications to the magnetic fields in tokamak plasmas. The initial seed for the flurry of research was the experimental observation that a small relaxation of the axisymmetric magnetic field into a non-axisymmetric, three-dimensional (3D) configuration can lead to mitigation or suppression of the ELMs in H-mode on the DIII-D tokamak [101]. ITER aims to mitigate ELMs by employing small external magnetic perturbations. However, since these perturbations can also impact plasma confinement, it is crucial to comprehend the influence of 3D fields on plasma transport. The objective is to identify a regime where Type I ELMs can be effectively mitigated without significant reductions in plasma performance as unintended consequences. The interpretation of experimental observations [227, 573, 574] complicates the understanding of the underlying physics. In certain plasma conditions, there is a noticeable depletion of plasma density (so-called 'pump-out'), potentially leading to reduced fusion efficiency in ITER, while the electron heat transport in the pedestal region remains low.
Non-axisymmetric plasmas have been thoroughly invest-
gigated in other (non tokamak) magnetic fusion devices like
stellarators and Reversed Field Pinches (RFPs). The effects
of 3D fields on the plasma are governed by a complex inter-
play between the plasma's response to the perturbation and
its impact on transport processes in both the plasma core
and the boundary. Therefore, it is crucial to gain a compre-
hensive understanding of how different transport channels are
affected by the external fields. Consequently, the investigation
of 2D topology breakdown in tokamaks has expanded beyond
the study of pedestal gradient relaxation at the plasma edge
to a broad spectrum of topics. It now encompasses a wide
range of magnetic field configurations that modify confine-
ment and affect neoclassical and turbulent transport proper-
ties. This review presents a comprehensive overview of the
literature, which deals with the influence of 3D fields on trans-
port in toroidal confinement devices. Specifically, we focus on
the recent advances in understanding how non-axisymmetric
field (external or intrinsic) affect neoclassical and turbu-
lent transport in tokamaks and stellarators. Additionally the
effects of external fields on pedestal and L-H transition are
presented.

**6.1.1. Structure of the paper.** 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.

**6.1.2. Magnetic field configurations.** There are many vari-
ations 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 mag-
netic 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 mag-
netic 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, how-
ever, 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 cre-
ation 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 trajector-
ies) [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 differ-
ent transport and turbulence characteristics, which indicates
that understanding the impact of 3D fields on TC requires a
good understanding of the underlying magnetic field struc-
ture. 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].

### 6.1.3. Experimental observations of confinement changes with 3D.
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 \times 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 \times 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 \times B$ shear and, thus, may influence each other.

## 6.2. Neoclassical transport

When non-axisymmetric magnetic fields are considered, either introduced externally or generated by saturated large-scale instabilities [605], neoclassical transport is no longer intrinsically ambipolar [606]. This means that the radial plasma current only vanishes for specific values of the radial electric field, $E_r$, and a corresponding rotation velocity called the offset velocity, $V_{offset}$. In this way, neoclassical physics strongly influences the rotation profile as described in section 6.2.1. The penetration of a resonant perturbation, and thus its impact on neoclassical transport, are tied to the ideal plasma response. The screening of RMPs on rational surfaces provides an additional rotation damping mechanism as described in section 6.2.2. Because of the role of the $E \times B$ shear rate in regulating turbulence, neoclassical and turbulent transport are intricately linked. Furthermore, the enhanced neoclassical transport due to 3D perturbations has a direct impact on particle transport, leading to the density pump-out effect and modification to impurity transport as described in section 6.2.3. Neoclassical transport in non-axisymmetric magnetic fields has been thoroughly reviewed by Shaing et al [399]. Because of their inherent 3D nature, neoclassical transport can have a substantial effect on the profile evolution in stellarators. For this reason, much of the theory for tokamaks with 3D perturbations was developed alongside the stellarator theory. Connections
between the two configurations will be highlighted in the following section.

### 6.2.1. Influence of 3D fields on NTV.
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_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_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-nq=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, $\nu$, in comparison with the bounce frequency, $\Omega_b$, $E \times B$ toroidal precession frequency, $\Omega_E$, and magnetic drift toroidal precession frequency, $\Omega_m$. The early theory of Shaing [612, 613] describes the low-collisionality regimes when $\nu < \Omega_E$ and $\nu > \Omega_E$, where the viscosity scales as $\sqrt{\nu-\nu}$ and $1/\nu$, respectively. When $\Omega_E > \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,

$$l\Omega_b + n(\Omega_E + \Omega_m) = 0, \quad (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_{\text{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 \times B$ shearing rate, see section 6.3.2, and enabling access to QH mode.

The radial profile of the torque depends on the collisionality regime and $E_r$ profile. The torque typically has a peak where $E_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_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
<figure><img src="image_29.png" /><figcaption>**Figure 26.** NBI is ramped toward zero at constant inject power. Without $n=3$ fields, rotation drops quickly, and QH mode is lost. With $n=3$ fields, larger rotation is maintained for same NBI torque and QH mode is maintained. Reproduced from [634]. © 2011 IAEA, Vienna. All rights reserved.</figcaption></figure>

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_r$ in the pedestal, this can impact the $E \times 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_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.

### 6.2.2. Resonant braking and screening effects at the boundary.
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