Comparison of momentum transport in matched hydrogen and deuterium H-mode plasmas in ASDEX Upgrade

Advanced momentum transport analysis is used to study matched hydrogen (H) and deuterium (D) plasmas in the core of ASDEX Upgrade. The aim is to validate gyrokinetic theory and assess a possible isotope dependence. The methodology extracts momentum diffusion, convection, and intrinsic torque as a function of time from experiments employing neutral beam injection (NBI) modulation. H and D plasma scenarios with comparable ion heat fluxes, NBI torque, electron densities, and several dimensionless parameters were designed to highlight any mass dependency. Linear gyrokinetic simulations predict that, for similar background gradients, the Prandtl and pinch numbers should be similar for H and D. This was confirmed by the experimental momentum transport analyses. The assessed intrinsic torques were found to be similar between H and D, co-current directed and located near the outermost region of the plasma core. The strength of the intrinsic torque is correlated with the amplitude of the plasma pressure gradient in the pedestal. Finally, a robust error analysis demonstrates the uniqueness of the parameters obtained together with their uncertainties. Neglecting the intrinsic torque, or its time dependence, systematically distorts the assessed momentum diffusion and convection. This is the first method to separate all three transport mechanisms from experimental data by retaining their time dependencies, that is found to match, quantitatively, the gyrokinetic predictions for Prandtl and pinch numbers, within experimental uncertainties.


Introduction
Understanding momentum transport is crucial to reliably predict the plasma rotation profiles in fusion devices.Rotation is known to influence impurity transport [1][2][3][4][5], contribute to the avoidance of MHD instabilities [6][7][8][9][10][11], and affect turbulence through E × B shearing [12][13][14].Understanding the sources, sinks, and transport of momentum is necessary to reliably predict the rotation profiles in future machines and its impact upon the stability and confinement of those plasmas.
While there have been several studies investigating mass and isotope effects upon momentum confinement [15,16], the majority of experiments focused on the dynamics of momentum transport in D plasmas [17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32].However, future reactor operation is foreseen with D and tritium (T) isotope mixtures and both helium (He) and H are being considered as working gases for initial experimental campaigns in ITER to avoid the complications of nuclear operation in the commissioning phase.Successful validation of our theoretical understanding of transport for different main ion species provides a more robust predictive capability.
This paper is presented as follows.In section 2, an advanced momentum transport analysis methodology is presented that can extract momentum diffusion, convection, and intrinsic torque by retaining the full time dependencies of all three mechanisms.Section 3 describes the performed NBI modulation experiments for H and D. Section 4 presents the assessed momentum transport coefficients and corresponding uncertainties.Alternative simplified transport models are benchmarked in section 5.The results are compared in section 6 with gyrokinetic simulations and, lastly, in section 7, a summary is given.

Methodology
The fundamentals of the methodology used in this work to study momentum transport have been described in an earlier work [33].This paper, thus, concentrates upon the changes to that methodology with the reader referred to the preceding publication.
In the earlier approach, the Prandtl number, the ratio of the momentum to the ion heat diffusivity (Pr = χ φ /χ i ), was fixed to gyrokinetic predictions.Herein, a considerably less constrained linear ansatz is taken: where ρ φ is the normalized toroidal flux coordinate given by ρ φ = √ (Φ 0 − Φ)/(Φ 0 − Φ sep ), Φ 0 the toroidal flux at the magnetic axis (r = 0) and Φ sep at the separatrix (definition of ρ ψ is similar with poloidal fluxes Ψ).The Prandtl number is constrained to be positive for 0 < ρ φ < 1, as a negative diffusivity is unphysical.In general, a more flexible ansatz could be implemented for cases where this does not represent the data well.However, more degrees of freedom come at the expense of increased computation and cross-correlations between regression parameters resulting in larger statistical uncertainties of the solutions.Therefore, it is paramount to retain as few regression parameters as possible while still matching the experimental data.
In the earlier work [33], the convective velocity V c was taken to scale with logarithmic density gradient R/L ne = − R ne ∂ r n e and the magnetic shear s.At ASDEX Upgrade (AUG), these quantities have large experimental uncertainties.Although there is strong theoretical motivation in their use to constrain V c [34], in this work, to allow an independent validation of theory, a less constrained free choice is taken for the pinch number in the form of a cubic polynomial, with the boundary condition set to vanish on the magnetic axis: The pinch number is constrained to be positive for 0 < ρ φ < 1, as an inward convective velocity is robustly predicted by theory for ion temperature gradient mode (ITG) and trapped electron mode (TEM) dominated conditions, as expected in these plasmas [34][35][36].
The modeled intrinsic torque is implemented as an additional flux, whose dominant contribution is thought to be turbulence driven residual stress.This flux emerges, especially within the plasma core, from turbulence [36][37][38][39][40], but, particularly at the plasma edge, other effects such as ion orbit losses [41], magnetic ripple [42,43], neutral drag [44], and neoclassical toroidal viscosity [45,46] contribute.Together, the sum of these fluxes translates to an intrinsic torque via τ int = −m i n i R 2 ∂V/∂ρ Π int with V the enclosed plasma volume at a given normalized radius ρ, see formalism in [33].Similar to the ansatz taken for the pinch number, the shape of the stress term is prescribed by with c s the sound speed.In contrast to the earlier work [33], the intrinsic torque now explicitly includes the time dependencies of χ i through its coupling with χ φ , which will be shown to be crucial in the time-dependent perturbation analysis.The same turbulence intensity that scales the diffusivity also dictates the amplitude of the residual stress, motivating this approach, see [39,47,48].The transport equations are implemented and solved with the 1.5D transport code ASTRA [49,50].The heat fluxes and NBI torque are calculated with the TRANSP code [51].The regression scalars C 1 . . .C 8 are subject to the statistical based global optimization algorithm 'Differential Evolution' [52], which scans a multi-dimensional, non-differentiable parameter space.They are varied to reproduce the experimentally measured steady-state, amplitude, and phase profiles of the plasma toroidal rotation.In contrast to the earlier work, no temporal smoothing is applied in the input to retain as much temporal information from the data as possible.Here, optimization is focused on reproducing only the first harmonic of the time dependence, but no role of higher harmonics was observed in the analysis.Higher harmonics could be included, for example in the case of low-frequency, asymmetric modulation.The chosen cost function is a reduced chi-squared measure with error weighting.In this cost function, the steady-state, amplitude and phase profiles have equal influence.The corresponding uncertainties in the regression parameters are displayed by error bars that reflect a one standard deviation in the underlying reduced chi-squared statistics.
The flexibility and time dependence of this model make it applicable to a wide variety of experimental conditions.Together with a robust error analysis, meaningful comparisons and validations of gyrokinetic predictions can be performed.Finally, the ability to extract an intrinsic torque in the presence of a simultaneous variation of the other momentum transport channels allows the systematic study of any parametric dependencies.The presented methodology is a new and powerful tool to analyze momentum transport in tokamak plasmas when a sufficient reaction to momentum modulation (here from a powerful neutral beam) can be diagnozed.

Experimental description
The objective of these experiments was to run comparable D and H discharges to isolate possible isotope effects in the deduced momentum transport.The analyzed discharges (H #41550 from 1.8 to 4.2 s and D #40076 from 2 to 4.2 s) were performed in well-established conditions, namely, type-I ELMy H-mode plasmas in a toroidal field of B φ = −2.5 T, a plasma current of I p = 0.8 MA, a line averaged core electron density of 6.5 • 10 19 m −3 , with a lower-single-null, favorabledrift, configuration.
The main observational diagnostics for this work are the charge exchange recombination spectroscopy (CXRS) systems [53,54], which provide the plasma impurity ion temperature and rotation on a 10 ms integration time base, that is, however, insufficient to resolve the sawteeth (f D ST ≈ 23 Hz and f H ST ≈ 46 Hz with an inversion radius of ρ φ ≈ 0.15) or the ELMs (f ELM ≈ 65 . . .70 Hz).The impact of these perturbations is effectively averaged out for the key diagnostic measurements.More importantly, neither the sawteeth nor the ELM frequencies were harmonics of the applied NBI modulation frequency (5 Hz), and so do not engender spurious higher harmonic artifacts in the modulation response.To minimize the impact of these perturbations on the momentum transport simulation, the region of interest in this work was limited to 0.2 < ρ φ < 0.8, clearly avoiding the radial regions which may be affected.
The plasma stored energy of the two discharges was comparable, W H MHD = 0.30 MJ and W D MHD = 0.36 MJ, as was total applied heating power, P H tot = 5.2 MW and P D tot = 6.6 MW.The main isotope fraction was ≈95% in the H discharge and ≈100% in the D discharge.For D, CXRS measured the main impurity boron and Z D eff ≈ 1.2 with the D gas flux set to Γ D ≈ 1 • 10 22 atoms s −1 .For the H discharge, the intrinsic B density was insufficient to provide high quality CXRS data, so nitrogen (N) was puffed with Γ N ≈ 0.35 • 10 21 atoms s −1 , which resulted in a Z H eff ≈ 1.38.In AUG, simultaneous CXRS measurements of N and B are common and show good agreement in both rotation and temperature in the plasma core, well within the measurement uncertainties.The use of two different impurity species to diagnoze the ion parameters is not expected to impact the analysis due to the very low impurity concentration of c N = 0.4% (in H) and c B = 0.3% (in D), which were calculated via the CHICA code [55].Previous experiments have shown that a significantly higher gas flux is required for H compared to D to achieve the same core density, due to significantly increased particle transport in the pedestal [56,57].Herein, the H gas puff was set to Γ H ≈ 2 • 10 22 atoms s −1 to obtain the same core electron and ion densities in both discharges.Figure 1(a) shows a comparison of the density profiles obtained.
As discussed earlier [33], NBI modulation is the method of choice to assess momentum transport coefficients in Hmode plasmas on AUG.The modulation of a single NBI source with reduced power (P H mod ≈ 1.1 MW and P D mod ≈ 0.7 MW) was a good compromise between providing a well resolved rotation modulation and limiting the perturbation in the other transport channels.In addition, off-axis NBI modulation was preferred, due to a lower induced ion temperature perturbation that, in turn, diminishes the time-dependent changes in all transport parameters.Using this NBI modulation scenario, the resultant perturbation of the electron density in the experiments was ≈0.5% at mid-radius in both cases.The T e (T i ) perturbations at mid-radius were ≈2(3)% and ≈4(4)% for D and H, respectively, while the stored energy, W MHD , varied by ≈5% (D) and ≈7% (H).The modulation of the rotation v φ is of the order of ≈8.5% (D) and ≈10% (H).Aside from their modulation, all these quantities' average values remain stable over the analyzed time range.Therefore, a time average over the analyzed time windows provides an appropriate comparison of the kinetic profiles and applied torque and heat fluxes between the discharges, see figures 1 and 2 respectively.Furthermore, as the density is near constant across the modulation, the change in the plasma angular momentum L φ = m i n i R v φ is only driven by changes in the rotation.This justifies choosing to perform the model-experiment optimization against the rotation velocity rather than the angular momentum, which is, of course, the physically conserved quantity.The temporal evolution of the density is tracked in the calculations.The electron density and temperature are inferred using the Integrated Data Analysis [58], based on lithium beam emission spectroscopy [59], laser interferometry [60], the measurement of the electron cyclotron emission radiometry [61] and the Thomson scattering diagnostic [62].
It is noteworthy that the error bars shown on n e and T e in this work represent the statistical uncertainty from the data reconstruction using Bayesian inference.The role of uncertainties in that analysis method is further discussed in [63].The error bars shown on T i and v φ are based on the statistical uncertainties from the fitting of the charge exchange spectra.These uncertainties are small due to the high cleanliness of the studied plasmas and the excellent conditions provided by N puffing in the H case for data fitting.The uncertainties on v φ additionally include systematic errors from calibration, which can occasionally result in jumps or dips in otherwise smooth velocity profiles.For T i , the systematic errors are negligible, as the dispersion relation of the spectrometer is very well known.It is important to mention that for all four quantities, the plotted uncertainties correspond to the average of the uncertainties for each radial position over time.
On-axis electron cyclotron resonance heating (ECRH) with P ECRH ≈ 0.6 MW was applied to both discharges to avoid impurity accumulation [64,65].In the D case, 4.8 MW of steady NBI heating with an extraction voltages of U D ex = 53 keV was applied.For the H case, the TRANSP code was then used to determine the NBI settings required to match the torque and heat fluxes to the D discharge.This approach was motivated by previous studies of the isotope dependence of heat transport [15,16,[66][67][68][69].In general, stronger ion heat transport is seen in H than in D discharges, for the same applied heating power.In [70][71][72], it was noted that the equipartition term's strong mass dependence influences the ion heat transfer.For the studied discharges, ion heating via electron collisions is a substantial contribution for H, while it is nearly negligible for D, see figure 2(c).This leads to larger Q i for H for the same heating recipe, which translates to χ H i > χ D i .However, the transport is similar when Q i is matched rather than the NBI power.
It was chosen to reproduce comparable torque and ion heat fluxes at the same density in both, H and D discharges.Based on scenario optimization with TRANSP, 3.4 MW of NBI was thus applied to the H discharge.This was achieved by two full voltage beams from box 1 (U D ex,1 = 52 keV).The modulating beam from box 2 was operated at reduced voltage (U D ex,2 = 63 keV).Together, this combination resulted in similar surface integrated ion heat fluxes as shown in figure 2(b).In addition, similar ion and electron temperature gradients were observed and measured in the two discharges, see figures 3(b) and (e).This is again suggestive of similar heat diffusivities in the core of the plasma, see figure 2(a), consistent with previous studies [70][71][72].The electron heat flux was not priorized in the scenario-optimization, but matched with a deviation of ≈20% at mid-radius.However, as this does not influence this analysis, it is not shown here for the sake of brevity.Following this procedure, the torque on the plasmas was successfully matched, see figure 2(d).
This matched external torque acted on plasmas with comparable target densities, see figure 1 plasma toroidal velocity, shown in figure 1(b), however, is only slightly higher for H. From the mass difference between H and D, for the same applied torque, one would expect a factor of two higher rotation velocity in the H discharge.This is not observed, as the torque in the H discharge is slightly lower, see figure 2(d), while the density is slightly higher, both of which contribute to reducing the observed rotation.As shown later, the remaining rotation value proximity is accounted for by small differences in momentum transport.

(a). The time-averaged
The ion and electron temperatures are slightly higher in the D discharge, figures 1(c) and (d).Despite comparable ion heat fluxes, shown in figure 2(b), the confinement at the pedestal is lower for H, because the inter-ELM transport increases with lower plasma particle mass [72].This is consistent with lower pedestal top temperatures in the H discharge, see figures 4(b) and (c), which then map to the lower temperatures in the core of the H discharge.The shapes of the n e pedestal profiles are also different, see figure 4(a), but have same pedestal top density.The resulting ion pressure profiles in figure 4(d) are higher for the D case and show steeper gradients, see figure 4(g).The error bars shown on gradients are derived from the regression of Gaussian processes, where the statistical uncertainties of the experimental profiles of panels (a)-(d) were used to approximate Gaussian errors.The error bars on the pressure, see figure 4(d), are calculated from analytical error propagation.On the kinetic profiles, the standard deviation over time for each radial position is shown with band structures.Considering that the standard deviation is comparable to or even smaller than the uncertainties associated with the kinetic profiles, it is worth noting that the propagation of the uncertainties also serves as an indicator of the stability of the provided profiles.
In these experiments, radial electric field E r measurements were not available, but the E r well in the pedestal can be estimated from the main ion diamagnetic term in the electric force balance [73,74].This analysis suggests a deeper E r well for the D case where the stronger E r gradients generate higher E × B velocity shear, which augment edge turbulence suppression [12], and E × B driven intrinsic torque [75].In figure 4, the edge CXRS T i data were artificially displaced radially inwards by 15 mm for the H discharge to match the region of steepest gradients with respect to the T e data from the Thomson scattering diagnostic.This displacement corresponds to a shift of 0.02 in the used coordinate ρ ψ .The need for this shift is assumed to be due to small errors in the equilibrium reconstruction, which affect the relative alignment of diagnostic measurements made at different toroidal and poloidal locations [76].Similarly, for the D discharge, the T i data was shifted radially inwards by 10 mm, corresponding to 0.015 in ρ ψ .As the gradients are not very sharp in these discharges, the relative alignment carries some uncertainty.
To isolate effects such as the isotope dependence, other, dimensionless, parameters should also be matched [77].A selection of such parameters is compared in figure 3.As for figure 4, the error bars shown on gradients in this figure are derived from the regression of Gaussian processes with the Gaussian errors being approximated with the statistical uncertainties of the input profiles.In the other dimensionless  parameters, the error bars result from analytical uncertainty propagation.In panel (a), the logarithmic density gradients are shown.The density gradients agree within the measurement uncertainties, due to large error bars in the underlying density profile, see figure 1(a).In both cases, the gradients increase in the inner half of the plasma radius while remaining flat outside mid-radius.The match inside mid-radius is not optimal, while reasonable agreement is obtained in the outer-half of the plasma.As will be shown in the next section, the convective momentum flux only contributes significantly to the transport towards the outermost region of the plasma core, where good agreement was obtained.
The logarithmic ion temperature gradient plays a key role in the ion heat diffusivity.The match is reasonable at most radii, as shown in figure 3(b), where the collisionality ν * is comparable up to ρ φ = 0.6.The difference results from the lower T e , as ν * ∼ n e /T 2 e .The normalized gyro-radius ρ * = ρ i /a = √ m i T i /eBa (with a the minor radius, ρ i the ion Larmor radius, the electron charge e) is lower for H.This is expected, as the ion mass m i and temperature T i are smaller for H.To match ρ * , the magnetic field and the plasma current would have to be modified to simultaneously keep the q profile similar.This would have made achieving a density match considerably more challenging.Due to limitations in experimental run time, it was decided to accept ρ * difference in favor of comparable densities and q profiles.No dependence in momentum diffusion or convection upon ρ * is expected, while a number of theoretical works predict the intrinsic torque to depend weakly on ρ * [48,78].The logarithmic electron temperature gradients are shown in figure 3(e).They are within error bars due to experimental uncertainties in the T e profiles, see figure 1(c), between 0.3 < ρ φ < 0.6, whereas, outside this region, a slight difference is found.The plasma β e ∼ n e T e /B 2 matches well within experimental uncertainties, see figure 3(f ), and the magnetic safety factor q is comparable in both discharges, as expected, from the similar equilibria.Finally, the temperature ratio T e /T i , which has a strong influence on the turbulence regime, is shown in figure 3(g).It agrees well and inherits rather large uncertainties around mid-radius from the experimental electron temperature profiles.
Together, this scenario offers a data set in which most of the relevant parameters, which govern transport and the underlying turbulence, agree in the core (where this analysis is concentrated) within experimental uncertainties.This allows the momentum transport analysis framework to be applied to this data with a reasonable expectation of isolating any large effects of the mass dependence on momentum transport.

Momentum transport modeling
The momentum transport analysis described in section 2 was applied to the D and H discharges.In both cases, it was possible to well reproduce the experimental steady-state, and modulation amplitude and phase profiles of the toroidal plasma rotation.Figure 5 shows an example for the D discharge (the H case is not shown as it is extremely similar).The experimental data is shown by brown solid lines with the modeled regression over-plotted with green dashed lines.All three profiles are well reproduced within the error bars.The 'noise' within the modeling seen in panel 5(d) is due to the absence of temporal smoothing of the experimental ion temperature and density profiles in the simulation.Beyond the modulation, also the temporal dependence of χ i and the torque in the model result from these variations.The rotation modulation can, however, be cleanly distinguished outside these variations.The momentum transport coefficients associated with the modeling are presented in figure 6, with the D results shown by orange dashed lines.In figure 7(a), the influence of diffusion, convection, intrinsic torque, and the externally applied NBI torque are shown in the same flux units for the D discharge to gauge the mechanisms' relative strengths.The diffusive flux is strongly balanced by the torque flux from NBI and the intrinsic flux.The convective flux plays a minor role in the plasma core, but becomes significant towards the outermost region of the core.The modulation amplitudes of these fluxes are shown in figure 7(b).The H case shows similarly consistent results that are again not shown for brevity.
In figure 6, the H results are displayed by blue solid lines for comparison with the D results.The Prandtl numbers, figure 6(a), and diffusion profiles, figure 6(b), are nearly identical, which is consistent with the ion heat diffusivities presented in figure 2. From theory, the momentum diffusivity is predicted to be closely related to the ion heat diffusivity, as both quantities are different moments of the velocity distribution [79,80].With little difference in χ i between the H and D case, little difference in χ φ is expected, as embodied by the general concept of limited variations of the Prandtl number, for different main ion masses.
The pinch number profiles are shown in figure 6(c) and mostly match within uncertainties.The tendency of higher core values for the H discharges may be connected to changes of the different density and temperature gradients in this region, see figures 3(a), (b) and (e).This would agree with theory: the pinch number is expected to scale most strongly with R/L ne , and, to lesser extent, with the temperature gradients [36].Within the accuracy of this methodology, no fundamental differences are thus observed between the D and H cases.The core deviation between H and D experiments will be discussed further in section 6, during the comparison with gyrokinetic predictions.
Finally, the intrinsic torque profiles are shown in figure 6(d).They also match, within uncertainties.In the plasma core, the intrinsic torque is found to be small with a strong co-current directed intrinsic torque required towards the outermost region of the plasma core to model the experimental data.The shape and magnitude of the intrinsic torque profiles are consistent with results from DIII-D [27].The intrinsic torque in the outermost region of the plasma core required to model the D discharge is slightly larger than for the H plasma. Referring back to figure 4(g), this correlates with the gradients of the pedestal ion pressure in these discharges.This is consistent with a scaling of the pedestal intrinsic torque with the ion pressure pedestal gradient presented by Solomon et al [19,27,81].There, it was suggested the co-current intrinsic torque could result from E × B driven residual stress in the inner shear layer of the edge E r well [75] and that the observed scaling of this torque with ∇p i is a proxy for E × B shear driven residual stress.Similar results have been also obtained in JFT-2M [82], Alcator C-Mod [29], NSTX [83], and LHD [84] where strong correlations with R/L Ti were also identified.Interestingly, herein, the normalized ion temperature gradients are quite similar such that changes in the ion pressure gradient (and consequently E × B shear) should be ascribed to the differences in the density pedestals.This is consistent with an E × B shear driven residual stress mechanism, rather than an ion temperature gradient main drive.However, two cases are insufficient to draw firm conclusions, particularly as other parameters such as ρ * and ν * , see figures 3(c) and (d), differ slightly.Rather, the scaling of the edge intrinsic torque as a function of local edge plasma parameters will be the subject of future work, engaging a wider experimental data set.
From the modeling performed here, it is clear that the experimental rotation modulation data is best reproduced when an edge localized intrinsic torque is included.However,  Assessed transport coefficients agree within error bars.The intrinsic torque shows slight differences for ρφ > 0.6, consistent with a variation in the edge pressure gradient.Pinch number profiles show variation inside error bars for ρφ < 0.4 which can be correlated to differences in the temperature and density gradients.the simulation boundary was set to ρ φ = 0.8 corresponding to ρ ψ ≈ 0.88, so no statements can be made about the pedestal of the plasma.To investigate a possible edge localized nature in the intrinsic torque, i.e. whether modeling is successful with a pedestal localized intrinsic torque alone, a numerical experiment was performed for the D discharge.The modeling boundary was set to ρ ψ = 0.97 and the residual stress in the model, polynomial part in equation ( 3), was substituted by a single Gaussian shaped profile with the position, width, and amplitude of that Gaussian left as free parameters.The diffusion and convection were extrapolated from the initial solution and maintained constant to isolate the effect in the modeled residual stress.The edge values of the Prandtl number, diffusion, and pinch number are shown in figures 8(a)-(c).This model yields an intrinsic torque with a maximum at the pedestal top, ρ ψ ≈ 0.93, shown in figure 8(d) by a green dotted line.For comparison, the extrapolated values from the cubic residual stress model are also shown.The modeled core steady-state rotation, phase, and amplitude profiles (not shown) are well reproduced within error bars with this residual stress ansatz showing that the Gaussian residual stress profile is an equally valid solution.This model again demonstrates that an edge localized intrinsic torque is really required to match the experimental data, but also that this technique cannot identify a unique radial distribution of that torque in the plasma edge region.The required integrated values, however, are very similar, suggesting that the technique has provided insight into the absolute magnitude of such a torque.The plotted intrinsic torque is not a symmetric Gaussian, as the scaling with χ φ in equation (3) must still be applied.This exercise did not succeed for H case, as the phase profile fluctuated too strongly at the edge to allow the radial boundary of the analysis to be set to a sufficiently high value.
In summary, the experimental analysis finds extremely similar core momentum transport coefficients and intrinsic torque between the matched H and D discharges.This suggests that, as expected, no fundamental isotope dependence exists for momentum transport within the inner plasma core.In the outermost region of the plasma core, the diffusion and convection remain similar out to a normalized radius of 0.8, while the intrinsic torque tends to displays differences, inside error bars, around ρ ϕ = 0.7, which are likely connected to differences in pedestal gradients between the discharges.
The D discharge used in this work is a repeat of the reference discharge of the data set examined in the earlier publication [33], but with a NBI modulation frequency of 5 Hz and using off-axis NBI torque modulation.The assessed transport coefficients for the D discharge agree, however, well with those assessed in the earlier publication.The diffusion and intrinsic torque also match quantitatively with the earlier results, whereas the pinch number was slightly lower in this work.This can be, potentially, ascribed to the changed convection ansatz employed or to the inclusion of an intrinsic torque, together with its time dependence, within the optimization.In general, this also demonstrates the robustness of the applied methodology, which is now less constrained than in the earlier work.
The results of these analyses contribute to a better understanding of the experimentally measured rotation profiles shown in figure 1(b).It was expected that, for the same transport levels, the H discharge would rotate at twice the D velocity due to the lower ion mass.As mentioned above, part of the explanation is the slightly higher torque from the NBI in the D case, compare figure 2(d).Additionally, the assessed momentum transport coefficients, with slightly lower intrinsic torque in the outermost region of the plasma core, contribute to explaining why the H rotation is found only ≈20% higher, demonstrating the sensitivity of the rotation to edge transport effects.
A robust global error analysis is an important feature of this methodology.As already described, the error bars in figures 5 and 6 reflect all the solutions with χ 2 red,norm < 1.5, corresponding to a ±σ environment of a Gaussian error distribution.For a better understanding of the shape of the parameter space spanned by the three interacting transport mechanisms, in figure 9 the logarithmic value of the cost function log(χ 2 red ) The fitting quality of the experimental profiles is nearly identical.There is a maximum of the intrinsic torque around ρ ψ ≈ 0.93.The Gaussian shape of the residual stress is nested with the diffusivity before the resulting intrinsic torque is calculated.This results in the slight asymmetry of the green dotted curve. is plotted.Here, the cost-function value depending on two of the three transport coefficients is shown with the third fixed to the best solution and, again for brevity, this is only shown for the D discharge at ρ φ = 0.5.In the first plot, due to experimental uncertainty in the fitted rotation data, there is a range of Pr that still reproduce the experimental data reasonably well.This uncertainty within the modeling results, to some degree, from experimental uncertainty.Moreover, the effects of the transport terms are not entirely orthogonal.This can lead, for example, to a possible balance between the outward diffusion and the inward convection.This can be seen as a valley-like blue structure in the cost function contour plots that is a reflection of the size of the modeling error bars.However, as the convective and the residual stress fluxes scale with χ φ , this uncertainty becomes coherent and cumulates.It increases further due to a possible balance the inward, co-current convection and a co-current intrinsic torque at higher radii.Therefore, the error bars shown on the modeling and the transport coefficients reflect the uncertainty derived from the uncertainties in the experimental rotation data and uncertainties in the methodology.As can be seen in the contour plots, there are no indications of multiple global minima, but only a single, well defined, best solution marked in these plots with a cross.

Validation of transport models
As this is the first work of its kind to include a time-dependent intrinsic torque in the modeling, it is of interest to probe the change in the result when neglecting the intrinsic torque and/or its time dependence.This is particularly informative as most previous studies of momentum transport applied at least one of these simplifications.Three numerical tests were performed on the same experimental data set to study this.First, only constant diffusion and convection were used with the intrinsic torque neglected.The resulting phase profiles are too flat and the obtained amplitude profiles were overstrong (not shown in figure 10 for clarity of the plot).In overview, the regressed fit strongly decreases the Prandtl number and diffusivity.The pinch, in turn, must decrease strongly to low, unphysical, values, even approaching zero over large parts of the radial profile.The corresponding cost-function Figure 10.Different modeling types of the D discharge #40076 (2.0-4.2 s).As shown in figure 5, the modeling with time-dependent intrinsic torque nicely agrees with experimental data, it is not again plotted for the sake of clarity.The solution without intrinsic torque does not reproduce the experimental data well.The solution with a constant intrinsic torque results in a comparable, though slightly worse modeling, especially for the phase profile (c) and a significantly different Prandtl number (d).
value and modeling quality are clearly inferior.This, together with previous theoretical arguments [39,47,85], demonstrates the need for time dependencies in the transport coefficients.
Second, the simulations included time-dependent diffusion and convection, but still neglect intrinsic torque.This analysis is shown in dotted dark-purple lines in figures 10(a)-(c).Green solid lines show the changed solution using the time-dependent intrinsic torque, that was shown earlier to accurately reproduce the experimental values.Modeling without intrinsic torque underestimates the steady-state profile of the rotation until ρ φ = 0.4 suggesting a missing cocurrent intrinsic rotation.In fact, based on the fit-quality criteria, this would not have been an acceptable modeling solution, χ 2 red ≈ 5.The Prandtl and pinch numbers are shown in figures 10(d) and (e) by dotted dark-purple lines.Neglecting the intrinsic torque causes the Prandtl number profile to be flat and lower and the pinch number to be slightly lower.The missing co-current intrinsic rotation cannot entirely be mimicked by the convection, as this would increase the modeled modulation amplitude profile even further.To compensate, the diffusion lowers to recover the steady-state rotation.The same observations are found for the H case, not shown in the interest of brevity.
Third, a temporally constant, but radially variable intrinsic torque was included.The model regression is shown by figures 10(a)-(c) with dashed brown lines.It is already better than without an intrinsic torque, but results in clear deviations in the amplitude and phase profiles, with the latter outside the experimental uncertainties.For the D case, the corresponding cost-function value is more than double that of the complete time-dependent result and the assessed transport coefficients result in modeled profiles outside experimental uncertainties.For the H discharge, the cost-function value is only slightly higher that of the complete modeling, as the intrinsic torque is smaller in the H discharge, and has less importance in the modeling.For D, it is the Prandtl number that is affected significantly.The assessment with constant intrinsic torque show similar profiles, but is higher for D and H.These differences are generally stronger towards the outermost region of the plasma core, where the effect of the time dependencies of what was found to be a more edge located intrinsic torque are strongest, compare figure 7(b).
It is probable that the over-prediction of Prandtl and pinch numbers results from the missed time dependence of the cocurrent intrinsic torque, which would contribute to higher modeled modulation amplitude profiles.The convection, which has to mimic this contribution, then becomes too large leading to over-steep phase profiles.As a result, the Prandtl number then also has to increase to re-flatten the phase profile again, as observed here.Finally, the constant intrinsic torque has to be increased over the full regression to compensate the increased diffusion, see figure 10(f ).Although modeling a constant intrinsic torque is clearly better than without any intrinsic torque, it remains inferior to the regression that includes its time dependence.As explained above, from a theoretical standpoint, there are strong arguments to invoke a time dependence in the residual stress from the modulated temperature associated with the neutral beam heating [39,47,48].
In general, it is possible that in earlier works, e.g.Tala et al [30], neglecting the intrinsic torque, and/or its time dependence, contributed to a systematic mismatch between gyrokinetic predictions and experimentally estimated transport coefficients.This comparison of different modeling approaches was expanded to a number of other AUG discharges (#29216 2-4.5 s, #34042 6.1-7.4 s, #39015 7.0-9.3s, #41551 6.4-7.6 s, #40076 4.2-6.4s) that cover a range of plasma parameters.Including time-dependent intrinsic torques led, in five out of seven cases, to a better match to the experimental data.No particularities were identified in the two cases where fitting with a constant intrinsic torque sufficed.Furthermore, modeling without intrinsic torque in general yielded lower Prandtl and pinch numbers and a clearly poorer match with experiment in all investigated cases.As described above, an artificially lower Prandtl number can result from the modeled diffusion compensating for the missing co-current contribution from the intrinsic torque.However, in other works, or for other discharges, the distortion of the transport coefficients could behave differently, e.g.convection could be used to compensate for the co-current intrinsic torque which would then lead to an over-estimation of the pinch number and, to compensate the effect upon the phase-and amplitude profiles, an over-estimation of the Prandtl number, as seen in [30] or observed for #32869 by the authors.However, all these possibilities depend upon the exact shape of the intrinsic torque and its precise time dependence.It should be noted that the exact modeling procedure, such as the choice of fitting radial range and the choice of cost function, can also influence the assessed results.In general, however, the generalized inclusion of time dependencies has made the regression less sensitive to these choices.
Together, this exercise has shown that the fitting with a time-dependent intrinsic torque reproduces the experimental data best and is, therefore, the model of choice for the following comparison with theory.

Comparison with gyrokinetic predictions
The experimentally assessed transport coefficients will be compared here with gyrokinetic predictions.The gyrokinetic flux tube code GKW [86,87] was used to predict the Prandtl and the pinch number, as in the earlier work [33].The fluxes were weight-averaged over a spectrum of five binormal wave numbers between 0.2 < k y ρ i < 0.9.These local, quasi-linear calculations were carried out for five radial positions to resolve a radial dependence of the Prandtl and pinch numbers.The residual stress predictions obtained from linear calculations arise from up-down asymmetry of the equilibrium and are very small in the plasma core [78].Higher order terms in ρ * , which can become important in the presence of strong background gradients, cannot be calculated from quasi-linear simulations, but require CPU-costly, global, non-linear calculations that are beyond the scope of this work [88][89][90][91].With nominal parameters from the experimental measurements, nearly all the calculations converged to well resolved unstable modes.The fastest growing instabilities were found to be ITG, except for the innermost radial point for H that is characterized by a real frequency in the electron drift direction, and could be identified as also influenced by TEMs. Figure 11 shows the variation of the predicted Prandtl and pinch numbers over the analyzed time window.For D, 10 randomly distributed time-points were sampled and resulting average and standard-deviation of the transport coefficients were plotted.Error bars reflect the numerical standard-deviation.This provides an estimation of the expected variance of the transport during the modulation as well as the potential variation in such calculations when averaged over this time window.A further calculation, with time-averaged profiles, agreed with the results obtained from sampling and averaging single time-points, indicating that the result is quite robust.Figure 12 shows the time-averaged calculations for the H case.
Taken together, the gyrokinetically assessed Prandtl numbers follow closely the experimental results agreeing within uncertainties in most radial positions for both isotope species, see figures 11(a) for D and 12(a) H.In both cases, experiment yields values slightly higher than the gyrokinetic predictions.However, one should note that the error bars on the experimentally assessed Prandtl numbers appear asymmetric at lower radii.This results from the positive Prandtl number fitting constraint even at small radii outside the range of this analysis.Combined with the linear ansatz for Pr, this constraint puts a lower limit on the diffusion towards the plasma center.Slightly lower values of the Prandtl number would, however, have been acceptable solutions based on the described fitting criteria, which is not reflected in the drawn error bars.This especially applies to the experimentally assessed Prandtl number for H.With this in mind, a corresponding gyrokinetic prediction for Pr for H is clearly recovered by the experiment within the uncertainties and assumptions intrinsic to the experimental modeling.In the earlier publication [33], the Prandtl number was fixed to gyrokinetic predictions where modeling were only shown to be consistent with the experimental data.Herein, Pr could choose any linear slope and offset, enabling direct comparison with, and validation of theoretical predictions.
The predictions for the pinch numbers are also recovered within uncertainties for most radial positions, as shown in figures 11(b) (D) and 12(b) (H).There remain slight differences in the experimental results between H and D in the core, as discussed in the previous section: the pinch number for H is higher than for D. It is encouraging to note that these differences are reflected in the gyrokinetic simulations.While for D, the GKW predicted pinch number for the innermost radial position is −R • V c /χ φ ≈ 0.77, it increases to ≈3 for H with only the innermost radial channel being TEM dominated for H.A first step in understanding the root of this difference comes by comparing the input of the gyrokinetic calculations at this radial position.Two main parameters appear.First, for D the ratio of the temperature gradients is smaller, R/L Te ≈ 10 < R/L Ti ≈ 6.0, while for H R/L Te ≈ 8.5 < R/L Ti ≈ 4.15, also shown in figures 1(b) and (e).Secondly, the density gradients are different, for D R/L ne ≈ 2.4, while for H R/L ne ≈ 3.2.
In three numerical experiments, the H calculations for the innermost radial point were repeated, but now replacing first the temperature gradients, then the density gradient, and finally both together, with their values from the D discharge.The aim is to separate the effects of the temperature and density gradients on the gyrokinetic predictions.One concludes that the TEM found in the core of the H case is temperature gradient driven, as replacing only the R/L Te,i values is sufficient for this region to remain ITG in H.This result has a  GKW prediction and experimental modeling for H discharge #41550 (1.8-4.2 s).Error bars on the Prandtl number appear to be asymmetric, because lower possible solutions were neglected as they would have lead to a negative Prandtl number for small radii due to the linear ansatz.However, one can expect that within the fitting domain, also slightly lower solutions would have been an acceptable solution.
significant effect on the pinch number, which decreased dramatically (−R • V c /χ φ ≈ 1.1) close to its value for D, where the dominant mode is always ITG.In the second experiment, replacing the density gradient still predicted a weak TEM core with a pinch number of −R • V c /χ φ ≈ 1.9, closer to the D result, but still nearly double the size.In the third experiment, where both temperature and density gradients were replaced, the D pinch was nearly recovered (−R • V c /χ φ ≈ 0.93 vs. 0.77 from D).This demonstrates that the differences between D and H can be attributed to differences in the kinetic profile gradients of their respective discharges and is not a result of an isotope dependence.
A comparison of gyrokinetic simulations for D and H, see figures 11 and 12, shows that similar Prandtl and pinch numbers are predicted, consistent with the experimental conclusion that there is no discernible isotope dependence in the core momentum transport, as expected from gyrokinetic theory.In a simple gyrokinetic picture, one neglects collisions and treats electrons adiabatically.These adiabatic electrons then follow the fluctuations in the parallel electric pressure directly.In this formalism, the gyrokinetic equations are naturally normalized to the ion mass and, for example, the diffusivity follows a gyro-Bohm dependence χ GB ∼ √ m i [92].This gyro-Bohm mass dependence is slightly modified when the electrons are treated kinetically.The ion to electron mass ratio is then taken into account, introducing the isotope dependent effects.The gyro-Bohm scaling also is broken when collisions are taken into account.The mass ratio then plays a substantial role in the calculation of the parallel wave-function.Therefore, in calculations with kinetic electrons and collisions, the gyro-Bohm dependence in the diffusivity is no longer cleanly recovered.This is the case for the calculations in this work.The differences seen between the prediction for H and D mainly result from changes in the plasma gradients whereas mass effects are small.To demonstrate this, a further numerical experiment was performed with GKW using D input data with a H electron to ion mass ratio.These results are indicated by circles in figure 11.Both the Prandtl and pinch number predictions are seen to increase slightly for the H mass ratio, but the effect remains significantly smaller than experimental uncertainties.It would, thus, not be possible to resolve such an effect with the present experimental methodology.
Overall, the good agreement found between experimentally inferred and theoretically predicted transport coefficients is very promising.In reference to the benchmark from fitting reduced different models in the last section, the inclusion of the intrinsic torque with its time dependence not only provides better experimental data modeling, but also transport coefficients that are in close agreement with theoretical predictions.

Summary
In this work, an advanced methodology was presented that can separate the three time-dependent components of momentum transport from NBI modulation experiments.A statisticalbased error analysis is included that demonstrates that a clear global optimization minimum in the spanned parameter space of the transport coefficients exists.This is, to the authors' knowledge, the first time that momentum diffusion, convection, and residual stress fluxes simultaneously are extracted from experimental data by maintaining the time dependence in all transport coefficients, as prescribed by the measured changes in the ion heat conductivity.This was used to study the mass dependence of momentum transport in the core plasma of the AUG tokamak.A H and D discharge pair was created as matched ion heat fluxes, NBI torque, and as comparable dimensionless parameters as possible.As shown, this resulted in similar heat transport and gradient scale lengths.
The experimental rotation data was modeled and the inferred momentum transport coefficients were found to be identical within experimental uncertainties for both isotopes.Therefore, one can conclude that no strong isotope effect, i.e. little difference in the momentum transport associated with the ion mass, is present in these experiments where observed differences were shown to be ascribable to the changes in the kinetic profiles of the discharges.This is consistent with the gyrokinetic calculations, which predict a weak mass dependence whose result would be smaller than the experimental uncertainties.GKW predictions for the Prandtl and the pinch numbers agree quantitatively with the experimental results, providing a direct validation of the gyrokinetic theory for these quantities over the experimental parameter range explored thus far.For the first time, a clear match of experiment and theory has been found in the context of momentum transport research.This work provides confidence that presentday gyrokinetic models are sufficient for modeling momentum diffusivity and convection in future devices, provided an appropriate boundary condition can be set.More work remains to understand momentum transport at the plasma periphery.
This analysis shows that neglecting the intrinsic torque or its time dependence results in regressed parameters that generate profiles outside experimental uncertainties.When not included, the resulting transport coefficients are distorted by trying to mimic the missing intrinsic torque.The intrinsic torque is found to be localized near the plasma edge and the two isotopes presented herein suggest that this torque is correlated with the edge ion pressure gradient and/or E × B shear.A better understanding of the generation and the parametric dependencies of the intrinsic torque at the edge will be goal of a future work.

Figure 1 .
Figure 1.Comparison of core kinetic profiles of the studied H and D discharges.The error bars reflect the statistical uncertainty of the data reconstruction, the band structures display the standard deviation over the analyzed time windows for each radial point.

Figure 2 .
Figure 2. The total ion heat fluxes (b), electron to ion equipartition term (c), and total torque (d) in the compared discharges resulting from the TRANSP calculations.The ion heat diffusivity (a) is based on ASTRA calculations.The band structures reflect the standard deviation over the analyzed time window for each radial point, mainly resulting from the modulation.

Figure 3 .
Figure 3.Comparison of dimensionless parameters in the plasma core.

Figure 4 .
Figure 4. Comparison of selected pedestal parameters.As radial coordinate, ρ ψ is chosen, which is defined outside the last close flux surface and resolves the plasma edge better.The arrows in panel (f ) illustrate the shift of the T i data to match the region of steepest gradients of Te.The standard deviation of the kinetic profiles over time is shown with band structures.

5 .
Experimental data and modeling of the D discharge #40076 (2.0-4.2 s) with a time-dependent intrinsic torque.

Figure 6 .
Figure 6.Comparison of the assessed momentum transport coefficients for both discharges, modeling with time-dependent intrinsic torque.Assessed transport coefficients agree within error bars.The intrinsic torque shows slight differences for ρφ > 0.6, consistent with a variation in the edge pressure gradient.Pinch number profiles show variation inside error bars for ρφ < 0.4 which can be correlated to differences in the temperature and density gradients.

Figure 7 .
Figure 7. Panel (a) shows a comparison of the estimated contributions to the total momentum flux from various transport mechanisms and sources, for the D discharge #40076 (2.0-4.2 s) with a time-dependent intrinsic torque.The quantities have units of momentum flux (shown in the y-label).A negative value corresponds to an inwards-directed flux, resulting in co-current rotation.The diffusive flux is mainly balanced by the torque flux from the NBI and the intrinsic flux.The black symbols show the sum of all four components which balance out all together.This demonstrates that the fluxes are consistent and lead to a stationary profile when being time-averaged.In panel (b), the amplitude of the modulation of the fluxes are shown.Effects of a time dependence are absolutely strongest towards the outermost region of the plasma core.

Figure 8 .
Figure 8. Numerical experiments investigating the edge localized nature of the assessed intrinsic torque.For fixed and extrapolated Pr and pinch numbers, the model for the residual stress was replaced by a Gaussian.While in the original fitting, the boundary was at ρ ψ ≈ 0.88, here it was set to ρ ψ ≈ 0.97, the boundaries are shown with corresponding vertical lines.The fitting result of the numerical experiment is shown by dotted line in panel (d).The fitting quality of the experimental profiles is nearly identical.There is a maximum of the intrinsic torque around ρ ψ ≈ 0.93.The Gaussian shape of the residual stress is nested with the diffusivity before the resulting intrinsic torque is calculated.This results in the slight asymmetry of the green dotted curve.

Figure 9 .
Figure 9. Cost-function value depending on fitted parameters at ρφ = 0.5 for the D discharge #40076 (2.0-4.2 s).The cross marks the best solution.

Figure 11 .
Figure11.GKW prediction (brown crosses) and experimental modeling for D discharge #40076 (2.0-4.2 s).Error bars on the predictions show the standard deviation of the coefficients over the analyzed time windows of the experiment.In addition, the GKW calculation from a numeric experiment in which the H mass ratio, instead of the D mass ratio, was used is shown by circles.This indicates the predicted effect of the isotope dependence.

Figure 12 .
Figure12.GKW prediction and experimental modeling for H discharge #41550 (1.8-4.2 s).Error bars on the Prandtl number appear to be asymmetric, because lower possible solutions were neglected as they would have lead to a negative Prandtl number for small radii due to the linear ansatz.However, one can expect that within the fitting domain, also slightly lower solutions would have been an acceptable solution.