Experimental determination of the three components of toroidal momentum transport in the core of a tokamak plasma

A new approach to infer the momentum transport in tokamak core plasmas via perturbation experiments is presented. For the first time, the analysis self-consistently includes all momentum transport components and their time dependencies, which are essential to separate the momentum fluxes and closely match the experiment. The quantitative agreement between the experimentally inferred transport coefficients and the gyrokinetic predictions provides an unprecedented validation. This work shows that the new methodology and gyrokinetic predictions can now be utilized on the route to physics-based prediction of momentum transport in future reactor plasmas.

Whereas particle and heat transport arise from diffusive and convective mechanisms, momentum transport has in addition a non-Fickian component (proportional to neither the rotation nor its gradient), which is usually called residual stress [19,21,22].The residual stress acts effectively as an intrinsic source of torque that can accelerate the plasma and thus substantially influence the resulting plasma rotation profile.The effective net toroidal rotation caused by this effect is referred to as intrinsic rotation [23][24][25].This is a remarkable characteristic of a self-organizing system, where small-scale turbulence dynamics can convert into large-scale flows.
Experimental validation of momentum transport theory requires concomitant measurements of diffusive, convective, and residual stress fluxes over a wide range of experimental conditions, which has proven challenging.Strong simplifications have often been employed in experimental analyses.One approach is to measure the intrinsic rotation or a related intrinsic torque over a range of plasma conditions and then extrapolate to future tokamak plasma [24,[26][27][28], based on the dependencies observed in present-day devices.However, such an approach does not provide direct information on the transport processes and neglects the impact that the individual components of the transport (diffusion, convection, and intrinsic torque) can have upon each other in the extrapolation to future devices, as a consequence of their different dependencies on fundamental plasma parameters.
A second approach is to design scenarios where the diffusion and convection can be neglected by using co-and countercurrent neutral beam injection (NBI) to zero the rotation profile, which allows to cleanly estimate the intrinsic torque.This is, among others, a technique used in [29,30].Based on these works, the authors discussed an extrapolation to a future fusion reactor [31,32].These kind of experiments can, however, only be performed in devices with balanced beams and no validation of the diffusive and convective fluxes is obtained.
The method of choice to separate all three fluxes is torque perturbation, where the external momentum source is varied temporally and the plasma reaction analyzed.A Fourier decomposition of the perturbation propagation is used together with the steady-state profiles and can separate all three components [33], although in many previous analyses an intrinsic torque was not consistently included [25,[33][34][35][36][37][38].Using high power NBI modulation as a time-varying momentum source, the applied torque can be estimated from beam dissipation codes, but the transport coefficients are, themselves, also affected by the absorbed heat and injected particles that cause the plasma kinetic profiles to vary in time.The inclusion of this time dependence becomes critical to the analysis, as will be shown herein, but was neglected in previous analyses [25,35,38].This resulted in discrepancies between the experimentally inferred and gyrokinetically predicted momentum transport coefficients [25,33,35,38] rendering extrapolation to future devices problematic.In this work, none of the mechanisms is fixed to theory predictions, the intrinsic torque is included and full time dependencies are retained in all of the transport coefficients and kinetic profiles, in contrast to the previous work [33] of the authors.This results in a stable, unique, and physics-based description of the momentum transport in the core of ASDEX Upgrade (AUG) H-mode plasmas.
The momentum transport modeling in this work employs the toroidal momentum conservation equation as formulated in [39].Here, the radial flux of toroidal momentum is written as with m the main ion mass of the plasma, n the main ion density, χ φ the momentum diffusivity, v φ the toroidal velocity (positive sign denotes co-current rotation), R the local major radius of the tokamak device, ∂/∂r the derivative with respect to the minor radius r, V c the convective velocity, and Π int an additional, intrinsic flux, which gives rise to an intrinsic torque τ int = −∂V/∂r Π int with V the enclosed plasma volume.Here, this flux is not specifically referred to as turbulent residual stress, even though that is expected to be the dominant contribution.As this torque is found to be edge localized, it may include contributions from non-turbulent mechanisms such as, e.g.ion orbit losses [40].The convective term in this model also includes possible contributions from particle flux, but from an analysis of the particle transport they are expected to be small in the experiments studied herein.
In the forward model, the number of free parameters is minimized by physically sensible choices for the radial dependencies of the three terms.The momentum diffusivity χ φ is taken to scale with the experimental ion heat diffusivity χ i via a linear function of the normalized toroidal flux radius ρ φ .Together, they define the Prandtl number Pr = χ φ /χ i , which is demanded to be positive for the entire radius, as a negative diffusivity is unphysical.The ion heat diffusivity is calculated via a power balance and is assumed to be representative of the effective transport present.Effects of stiff or non-linear heat transport are not explicitly corrected.The methodology assumes that the χ φ calculated in this manner is comparable to the incremental momentum diffusivity obtained from partial derivatives, which would guarantee that the diffusivity obtained from the analysis of the modulation is also representative for a steady-state description.This assumption is consistent with the observed temporal responses of the rotation profile to the applied torque modulations for the analyzed discharges.The convection (via the pinch number −R V c /χ φ ) and the intrinsic flux (via Π int /χ φ ) profiles are prescribed as radial cubic polynomials, with the condition to vanish at ρ φ = 0 to guarantee continuity there.These functions were selected to provide sufficient flexibility to closely reproduce the experimental dynamics whilst employing as few parameters as possible.Together, this model contains 8 free parameters.
It is important to note that through Pr, χ φ also retains a time dependence on the experimentally estimated ion heat diffusivity χ i = −Q i /(n ∂T i /∂r) (Q i the ion heat flux, T i the ion temperature), which enters explicitly into the formula for both the pinch number and the intrinsic torque.Including a scaling with χ i is an essential proxy for a time-dependent scaling with the turbulence intensity, which not only scales the momentum diffusion, but also dictates the residual stress amplitude [39,41] and compensates for the cross-channel coupling induced from NBI modulation.
The momentum conservation equation is solved using the transport code ASTRA [42,43], where the heat fluxes and neutral beam torque are calculated using the Monte Carlo code NUBEAM [44], which is part of the TRANSP code suite [45].The NUBEAM results were observed to be stable under small changes to the input kinetic profiles.ASTRA evolves the toroidal rotation profile, based upon a set of transport coefficients and an experimental rotation boundary condition at the edge of the plasma core at ρ φ = 0.8.The transport coefficients are varied to best match the experimental steady-state, amplitude, and phase profiles of the toroidal plasma rotation.There, the amplitude and phase profiles of the first harmonic of the modulation are used; no signature of higher harmonics was found in the analysis of the experimental data.The steady-state profile corresponds to the time-averaged toroidal rotation profile of the considered time-span.Among these three profiles, no hierarchy is applied in optimizing the experimental response to the model parameters unlike that, for example, in [35,38], as this can severely affect the inferred transport coefficients.The analysis framework includes a statistical error analysis, relying upon an error-weighted χ 2 red cost function.Error bars on the fit reflect one standard deviation from the best solution.The variation is mapped out by randomly varying the fitting parameters, tracking changes of the cost function up to a variation of 1.5 in χ 2 red , which is equivalent to one standard deviation.By doing so, also a wide parameter space was scanned to ensure that the final transport coefficients represent a unique, global solution.
The fitting procedure employed in this study is based on the fundamental assumption that the perturbed Fourier profiles contain additional information about the steady-state profiles of diffusion and convection, and that the perturbation does not alter the background turbulence.Unlike heat diffusivity, which is influenced by temperature gradients through turbulence, the momentum transport coefficients are not expected to strongly depend on the rotation or its gradient within the range they are varied in this study.This represents a fundamental assumption of the analysis methodology.
This assumption was previously examined theoretically for the Prandtl number [46].In that study, the predicted dependence of the Prandtl number on u ′ , the normalized rotation gradient, was assessed.In the experimental data set studied herein, the effect of the modulation is only between 10% and 15% of the absolute value of u ′ ≈ 0.2-1.5 for all cases and radial positions.When comparing these values to the findings in [46], it becomes evident that the changes in u ′ during the modulation are too minor to significantly affect the Prandtl number.The authors assume this also applies to the ratio of convective velocity to momentum diffusivity, the pinch number, and the intrinsic torque.This is supported by the fact that the gyrokinetic calculations, which will be discussed later in this letter, computed growth rates an order of magnitude larger than the E × B shearing rates for the cases under investigation.Therefore, a modulation of 10%-15% of the absolute value of u ′ , directly translating into the shearing rate, is assumed to have minimal to no impact on the turbulence state.Additionally, it is worth noting that the observed absolute phase relationship between the rotation gradient and the ion heat diffusivity is nearly in phase across most of the radial domain.
If momentum transport coefficients were to depend on the steady-state rotation or its gradient, it could create an underconstrained fitting problem, leading to significant uncertainties in the coefficients.However, as later demonstrated, this is not the case.Furthermore, modeling steady-state and timedependent quantities is feasible using the same transport coefficients.This suggests the validity of the fundamental assumption that additional information about the steady-state transport coefficients can be extracted from the Fourier profiles.
This analysis framework was used to analyze a deuterium H-mode discharge in the AUG tokamak (#40076 from 2.0 to 4.2 s) with a magnetic field of 2.5 T, a plasma current of I p = 0.8 MA, and a line-averaged density of 6.4 × 10 19 m −3 .The plasma was heated with 0.6 MW of electron cyclotron resonance heating and 4.8 MW of steady, on-axis, NBI heating (extraction voltage U ex = 53 kV).A torque perturbation was induced by a reduced power off-axis oriented NBI with a modulated power of 0.7 MW (U ex = 50 kV).A 5 Hz, symmetric duty cycle beam modulation was applied and its effect on the plasma rotation was tracked using charge exchange recombination spectroscopy [47,48] of the intrinsic boron impurity on a 10 ms integration time base, which is sufficient to resolve the modulation.The modulation is slow compared to the J × B and collisional slowing down time inferred from the NUBEAM results.Herein, the impurity rotation is assumed equal to the main ion rotation, i.e. neoclassical corrections are not applied.Previous works on a similar discharge showed that possible offset corrections can, partly, be compensated for via the boundary condition given [33].Neoclassical corrections in gradients were observed to be smaller than experimental uncertainties and have no effect on the assessed Pr and −RV c /χ φ .For discharges with strong differences in profiles gradients, this assumption should be reassessed.
The measured toroidal rotation was modulated by ≈10%.Other key quantities were less perturbed, i.e. the plasma-stored energy (≈5%), and the ion temperature and heat diffusivity (≈3%).With a modulation of the density profiles at mid radius of ≈0.5% and stable density gradients, no modulation of the particle transport is expected.As a result, the changes in plasma momentum are dominated by the changes in rotation velocity.This justifies a model minimization against the easy measurable rotation velocity, as opposed to the angular momentum, which is the physically conserved quantity.The region of interest was set to 0.2 < ρ φ < 0.8, to minimize the effect of sawtooth and edge-localized mode plasma instabilities.
Experimental data and modeled results are shown in figure 1 for the steady-state, amplitude, and phase profiles (Panels (a)-(c)) of the plasma toroidal rotation.The excellent agreement is a considerable improvement compared to previous results, e.g. in [25,35,36,38].Corresponding momentum transport coefficients are shown in Panels (d)-( f ).The extracted Pr (Panel (d)) is of order 1, as expected [17], and increases with minor radius.The pinch number (Panel (e)) increases from mid-radius towards the outer core.In Panel ( f ), the assessed intrinsic torque is shown that is found to be near zero in the core.However, a strong, edge-localized, co-currentdirected intrinsic torque is required to reproduce the experimental data, that is consistent with previous observations from balanced beam experiments [29,30], but in contrast to those works, here, the intrinsic torque is obtained together with its interaction with diffusion and convection.Moreover, within the analysis for this work, it was observed that an intrinsic torque in the outer core, through inward convection and diffusion, leads to an effectively co-current intrinsic rotation in the inner core.This underlines the importance of diffusion and convection on the formation of rotation shearing and in setting the profile shape, also by their interaction with the intrinsic torque.This interaction highlights the need for disentangling the transport mechanisms instead of concentrating solely on the scaling of intrinsic rotation, which might hide dependencies that significantly impact, for example, extrapolations.
Figure 2(a) shows the assessed momentum fluxes.The diffusive flux is mostly balanced by the torque from NBI and the intrinsic momentum flux.The sum of these fluxes (in black symbols) is zero, leading to stationary profiles when time averaged.This illustrates the consistency of these results.In figure 2(b), the modulation amplitudes of the same fluxes are plotted.The time-dependent effect of the perturbation is strongest towards the edge of the plasma, but remains non-negligible in the plasma core.One can see that the NBI torque flux has, by far, the largest modulation relative to its time-averaged values.In figure 2(c), the phase profiles of the fluxes are shown.Their values were shifted to zero at the simulation boundary to facilitate comparison.One can see that most fluxes have very similar phase profiles, with the exception of the diffusive flux, which picks up some additional delay, due to its dependence on the gradient of the rotation, which has a larger phase delay than the rotation itself.This difference is likely caused by variations in the time scales of the mechanisms involved.Firstly, switching on the modulating beam applies an almost instantaneous J × B torque to the plasmas, which increases the overall value of the rotation profiles and the associated convective momentum flux.Subsequently, collisional torque coupled with transport leads to a peaking of the rotation profiles on a delayed time scale.This underlines the importance of the diffusive process and the interaction of the mechanisms for the modeling of the phase profiles.
Next, other fitting approaches are tested to vindicate the inclusion of time dependencies in all terms.First, a model with constant diffusion and convection and no intrinsic torque is applied (not shown in figure 1 for clarity).A second variant only removed the intrinsic torque, but retained time-dependent diffusion and convection, shown by purple dashed lines in figure 1.In the third test, a temporally constant intrinsic torque was employed with the results shown in dash-dotted pink lines.All resulting fits were statistically significantly poorer with extracted coefficients drastically different from the complete fit especially in the outer core where the effect of time dependencies is strong, see figure 2(b).Although additional ingredients improved the experimental agreement, only the complete fit showed agreement consistent with experimental uncertainties.It is salient to repeat that adding these time dependencies does not introduce additional fit variables, so, the improvement is not a result of over-fitting.The authors speculate that neglecting the intrinsic torque or time dependencies in previous works caused a systematic, and often unphysical, deviation in the inference of the transport coefficients.Here, the free inclusion of an intrinsic torque that scales with χ φ reproduces the experimental data well, whilst remaining consistent with theory [39,41].
Local, quasi-linear, gyrokinetic calculations were performed with the GKW code [49] for comparison of the experiment with theoretical predictions.The fluxes were averaged over a spectrum of five binormal wavenumbers between 0.2 < k y ρ i < 0.9.This generates values for Pr and pinch numbers, but not for the intrinsic torque, which would require global, non-linear, calculations that are beyond the scope of this work.Also, one cannot directly compare the unnormalized fluxes, e.g.χ i or χ φ independently, between experiment and gyrokinetic prediction, as the gyrokinetic fluxes quasi-linear scale with the growing potential fluctuation in the calculation.All cases reported herein converged to unstable modes where the ion temperature gradient (ITG) mode dominated.The calculated Pr, see black symbols in figure 1(d), agree with the experimental data within error bars.Also, the increase in Pr with increasing major radius is reproduced and the linear trend agrees with earlier predictions [17].The predicted pinch numbers also show good agreement, see figure 1(e).As very few physical constraints were imposed upon the profiles of these quantities in the modeling, these observations should be taken as a validation of the gyrokinetic predictions for this discharge.
To generalize this validation, further ITG-dominated H-mode deuterium and hydrogen plasmas were studied.Main plasma parameters are shown in table 1.These were chosen to represent a variation in plasma current and consequently safety factor q and shear s.They differ in heating mixture and cover a range of line-averaged core densities ne and normalized density gradients (R/L n ≈ 0-4).Both shear and density gradients are predicted to be key quantities for the convection of momentum [15,20], so these discharges provide a strong test for assessing the sensitivity of this method.
The momentum transport components were extracted using the complete analysis framework.The need to include time dependencies was again confirmed by these additional cases, completely analogous to the results presented in figure 1.The pinch number profiles (middle column of figure 3) are all peaked, with the exception of #39015 which has the lowest density gradients in the plasma inner core and contains no ECRH in contrast to the other discharges.The derived intrinsic torque profiles (right column of figure 3) increase towards the outer core, although for the discharge #29216 in the bottom row, it remains close to zero over the entire minor radius.This is attributed to lower pedestal pressure and pedestal gradients of this discharge.The intrinsic torque at ρ φ = 0.7 shows a clear correlation with the pedestal pressure gradients, see figure 4(a), similar to observations in [30,50].
Asymmetries in the shown error bars can be caused by the conditions set for Pr and −RV c /χ φ to be positive, as a negative diffusion or outward pinch would be unphysical.This condition holds for the entire radial range of the ASTRA simulation, not just the range over which the analysis was done.This influences in particular the Pr with its linear shape.For example, for #41551, slightly lower Prandtl number values could have been acceptable solutions according to the fitting criteria described in the shown radial range, but this is not reflected in the displayed error bars.GKW predictions match experiment for all analyzed discharges, see figure 3.They recover the factor of 3 variation in the experimental Pr and reproduce well the variations in pinch profile shapes and magnitudes seen in the experiment.Figure 4(b) shows the predicted and experimental pinch numbers against the density gradient at ρ φ = 0.35, where the largest variation of the logarithmic density gradient is obtained in the available experimental dataset.The predicted pinch number scaling with the density gradient [15,20] is mirrored by the experiment.
In summary, a new methodology was presented that is able to extract momentum diffusion, convection, and intrinsic torque over a wide range of tokamak plasma discharges.Maintaining time dependencies in all transport coefficients is shown to be crucial to reliably decouple the non-linear effects of the induced perturbations.Only then can the experimental data be matched by the transport model, accurately separating the momentum fluxes.The determined Prandtl and pinch numbers show quantitative agreement with gyrokinetic predictions, for the first time in such an analysis.This represents a milestone in the study of turbulent momentum transport, as previously reported discrepancies between gyrokinetic predictions and experimental transport coefficients are now resolved, adding credibility in exploiting gyrokinetic theory to predict momentum transport in future devices.As the intrinsic torque has been included consistently in the modeling, this work opens a new experimental possibility for the validation of intrinsic torque predictions and integrated modeling approaches.Next, essential, steps are to revisit previously assessed parameter scalings of the transport coefficients and to probe more extreme plasma conditions, e.g. with larger fraction of electron heating, closer to conditions in a reactor plasma.This should now lead to the first, consistent, physicsbased and validated predictions for a future reactor scenario.A first detailed application of this methodology is demonstrated in [51].

Figure 1 .
Figure 1.(a) Steady-state, (b) modulation amplitude, (c) and phase profile of the rotation induced in discharge #40076 (2.0-4.2 s).Experiment is shown in brown solid lines and complete modeling in green dashed lines.Alternative modeling discounting the intrinsic torque is shown as dashed purple, while modeling including intrinsic torque, but neglecting its time dependence, is shown as pink dash-dotted lines.Panels (d)-( f ) show the fitted Pr, pinch number, and intrinsic torques.The gyrokinetic predictions for Pr and pinch numbers are shown as black circles.

Figure 2 .
Figure 2. (a) Modeled contributions to the time-averaged momentum flux for the discharge #40076 (2.0-4.2 s).A negative flux is co-current/inward directed.The sum (crosses) shows that all components balance, demonstrating the consistency of the assessed coefficients.(b) The modulation amplitude of the same fluxes increases towards the outer core.(c) The modulation phase shifts are similar for the fluxes with the exception of the diffusive flux, which incorporates strong effects of its scaling with the rotation gradient.

Figure 3 .
Figure 3.Comparison of gyrokinetic predictions (circles) and modeling (dotted lines) of Pr and pinch numbers (left and central columns) for further tested discharges.The modeled intrinsic torque is shown in the right column.Error bars for #39015 are very large in the outer core, which is why they were not displayed in favor of clarity.

Figure 4 .
Figure 4. (a) Scaling of the intrinsic torque at ρφ = 0.7 plotted versus the edge pressure gradient in the steepest gradient region.#39015 has not been included as it features very large error bars in the outer core.Linear fit to guide the eye.(b) Scaling of the pinch number with the density gradient for ρφ = 0.35, where variation of the density gradient is largest.Experimental data is shown as green crosses, gyrokinetic prediction as open black symbols, matching within error bars.The linear fit of the experimental data (dashed green line) reproduces the linear fit of the gyrokinetic predictions (solid black line).

Table 1 .
Comparison of the main plasma parameters of the investigated discharges.