Rotation and Neoclassical Ripple Transport in ITER

TL;DR

This study models neoclassical ripple transport and torque in ITER, showing ferritic inserts significantly reduce toroidal viscosity and impact plasma rotation, with detailed calculations accounting for various plasma and magnetic field factors.

Contribution

It provides the first comprehensive calculation of NTV in ITER considering realistic 3D magnetic equilibria, ripple effects, and ferromagnetic structures, improving understanding of plasma rotation damping.

Findings

NTV torque is dominated by the 8 ripple component.

Ferritic inserts reduce NTV torque by about 75% near the edge.

Ripple-trapping may be a significant NTV mechanism in ITER.

Abstract

Neoclassical transport in the presence of non-axisymmetric magnetic fields causes a toroidal torque known as neoclassical toroidal viscosity (NTV). The toroidal symmetry of ITER will be broken by the finite number of toroidal field coils and by test blanket modules (TBMs). The addition of ferritic inserts (FIs) will decrease the magnitude of the toroidal field ripple. 3D magnetic equilibria with toroidal field ripple and ferromagnetic structures are calculated for an ITER steady-state scenario using the Variational Moments Equilibrium Code (VMEC). Neoclassical transport quantities in the presence of these error fields are calculated using the Stellarator Fokker-Planck Iterative Neoclassical Conservative Solver (SFINCS). These calculations fully account for $E_r$, flux surface shaping, multiple species, magnitude of ripple, and collisionality rather than applying approximate analytic NTV formulae. As NTV is a complicated nonlinear function of $E_r$, we study its behavior over a plausible range of $E_r$. We estimate the toroidal flow, and hence $E_r$, using a semi-analytic turbulent intrinsic rotation model and NUBEAM calculations of neutral beam torque. The NTV from the $\rvert n \rvert = 18$ ripple dominates that from lower $n$ perturbations of the TBMs. With the inclusion of FIs, the magnitude of NTV torque is reduced by about 75% near the edge. We present comparisons of several models of tangential magnetic drifts, finding appreciable differences only for superbanana-plateau transport at small $E_r$. We find the scaling of calculated NTV torque with ripple magnitude to indicate that ripple-trapping may be a significant mechanism for NTV in ITER. The computed NTV torque without ferritic components is comparable in magnitude to the NBI and intrinsic turbulent torques and will likely damp rotation, but the NTV torque is significantly reduced by the planned ferritic inserts.

Equations

$${\delta }_B=(B_{\max }-B_{\min })/(B_{\max }+B_{\min }),$$

$$\nabla \cdot {\Pi }_{\zeta }^{\mathrm{turb}}({\Omega }_{\zeta })+\nabla \cdot {\Pi }_{\zeta }^{\mathrm{NC}}({\Omega }_{\zeta })={\tau }^{\mathrm{NBI}},$$

$${\Pi }_{\zeta }^{\mathrm{turb}}=-m_i{n}_i{\chi }_{\zeta }⟨R^2⟩\frac{\partial {\Omega }_{\zeta }}{\partial r}+{\Pi }_{\mathrm{int}}.$$

$$⟨...⟩=\frac{1}{V^{\prime }}{\int }_0^{2\pi }d\theta {\int }_0^{2\pi }d\zeta \sqrt{g}(...)$$

$$V^{\prime }={\int }_0^{2\pi }d\theta {\int }_0^{2\pi }d\zeta \sqrt{g},$$

$$-m_i\frac{1}{V^{\prime }}\frac{\partial }{\partial r}\left(V^{\prime }{n}_i{\chi }_{\zeta }⟨R^2⟩\frac{\partial {\Omega }_{\zeta }}{\partial r}\right)=-\frac{1}{V^{\prime }}\frac{\partial }{\partial r}\left(V^{\prime }{\Pi }_{\mathrm{int}}\right)+{\tau }^{\mathrm{NBI}}.$$

$${\tau }^{\mathrm{NBI}}=-\frac{1}{V^{\prime }}\frac{\partial }{\partial r}\left(V^{\prime }{m}_i{n}_i{\chi }_i⟨R^2⟩\frac{\partial {\Omega }_{\zeta }}{\partial r}\right).$$

$${\Omega }_{\zeta }(r)=-{\int }_r^a\frac{v_{ti}{\rho }_{\ast ,\theta }}{2P_r{L}_T^2}\tilde{\Pi }({\nu }_{\ast })d{r}^{\prime },$$

$$\tilde{\Pi }({\nu }_{\ast })=\frac{({\nu }_{\ast }/{\nu }_c-1)}{1+({\nu }_{\ast }/{\nu }_c)},$$

$$V_{||}^a=\left(\frac{1}{n_a}\right)\int d^{3}v{v}_{||}{f}_a,$$

$$(v_{||}\mathbf{b}+{\mathbf{v}}_E+{\mathbf{v}}_{\mathrm{m}a})\cdot (\nabla f_{a1})-C(f_{a1})=-{\mathbf{v}}_{\mathrm{m}a}\cdot \nabla \psi \left(\frac{\partial f_{a0}}{\partial \psi }\right)+\frac{Z_{a}e{v}_{||}B⟨E_{||}B⟩}{T_a⟨B^2⟩}{f}_{a0}$$

$${\mathbf{v}}_E=\frac{1}{B^2}\mathbf{B}\times \nabla \Phi$$

$${\mathbf{v}}_{\mathrm{m}a}\cdot \nabla \psi =\frac{1}{{\Omega }_{a}B}\left(v_{||}^2+\frac{v_{\perp }^2}{2}\right)\mathbf{b}\times \nabla B\cdot \nabla \psi ,$$

$$⟨V_{||}^{a}B⟩=-\frac{G}{Z_{a}e{n}_a}\left[\frac{1}{n}\frac{d(nT)}{d{\psi }_P}+Z_{a}e\frac{d\Phi }{d{\psi }_P}-k_{||}\frac{dT}{d{\psi }_P}\right],$$

$$J_{BS}=\sum _a{n}_a{Z}_{a}e⟨V_{||}^{a}B⟩,$$

$${\Gamma }_{\psi ,a}=⟨\int d^{3}v({\mathbf{v}}_{\mathrm{m}a}\cdot \nabla \psi )f_a⟩,$$

$${\tau }^{\mathrm{NTV}}=-B^{\theta }\sqrt{g}\sum _a{Z}_{a}e{\Gamma }_{\psi ,a},$$

$${\tau }^{\mathrm{NTV}}=-{\nu }_{\zeta }⟨R^2⟩mn({\Omega }_{\zeta }-{\Omega }_{\zeta ,\mathrm{offset}}),$$

$$B=\sum _{m,n}{b}_{mn}^c\cos (m\theta -n\zeta )+b_{mn}^s\sin (m\theta -n\zeta ),$$

$${\mathbf{v}}_{\mathrm{m}a}=\frac{v^2}{2{\Omega }_a{B}^2}(1+{\xi }^2)\mathbf{B}\times \nabla B+\frac{v^2}{{\Omega }_{a}B}{\xi }^2\nabla \times \mathbf{B},$$

$${\mathbf{v}}_{\mathrm{m}a}^{\perp }=\frac{\nabla \psi \times ({\mathbf{v}}_{\mathrm{m}a}\times \nabla \psi )}{|\nabla \psi {|}^2}=\frac{v^2}{2{\Omega }_a{B}^2}\frac{(\mathbf{B}\times \nabla \psi )}{|\nabla \psi {|}^2}\nabla \psi \cdot \left[(1-{\xi }^2)\nabla B+2B{\xi }^2(\mathbf{b}\cdot \nabla \mathbf{b})\right].$$

$${\mathbf{v}}_{\mathrm{m}a}^{\perp }=\frac{v^2}{2{\Omega }_a{B}^2}(\mathbf{B}\times \nabla \psi )(1-{\xi }^2)\frac{(\nabla \psi \cdot \nabla B)}{|\nabla \psi {|}^2}.$$

$${\int }_0^{r}dV(r^{\prime })H(r^{\prime })={\int }_0^{r}dV(r^{\prime })\left(P_{\mathrm{LH}}+P_{\mathrm{NBI}}+P_{\mathrm{ECH}}+R_{\mathrm{DT}}(3.5\mathrm{MeV})\right).$$

$$Q(r)=\frac{{\int }_0^{r}dV(r^{\prime })\left(P_{\mathrm{LH}}+P_{\mathrm{NBI}}+P_{\mathrm{ECH}}+R_{\mathrm{DT}}(3.5\mathrm{MeV})\right)}{A(r)},$$

$${\Pi }_{\mathrm{int}}\sim \frac{{\rho }_{\theta }\tilde{\Pi }({\nu }_{\ast })Q⟨R⟩}{v_{ti}{L}_T}.$$

Full text

I Introduction

Toroidal rotation is critical to the experimental control of tokamaks: the magnitude of rotation is known to affect resistive wall modes Bondeson and Ward (1994); Garofalo et al. (2002), while rotation shear can decrease microinstabilities and promote the formation of transport barriers Burrell (1997); Terry (2000). As some ITER scenarios will be above the no-wall stability limit Liu et al. (2004), it is important to understand the sources and sinks of angular momentum for stabilization of external kink modes. One such sink (or possible source) is the toroidal torque caused by 3D non-resonant error fields, known as neoclassical toroidal viscosity (NTV). Dedicated NTV experiments have been conducted in the Mega Amp Spherical Tokamak (MAST) Hua et al. (2010), the Joint European Tokamak (JET) Lazzaro et al. (2002); de Vries et al. (2008), Alcator C-MOD Wolfe et al. (2005), DIII-D Garofalo et al. (2008); Reimerdes et al. (2009), JT-60U Honda et al. (2014), and the National Spherical Tokamak Experiment (NSTX) Zhu et al. (2006).

In addition to the ripple due to the finite number (18) of toroidal field (TF) coils, the magnetic field in ITER will be perturbed by ferromagnetic components including ferritic inserts (FIs) and test blanket modules (TBMs). TBMs will be installed in three equatorial ports to test tritium breeding and extraction of heat from the blanket. The structural material for these modules is ferritic steel and will produce additional error fields in response to the background field. The TBMs will be installed during the H/He phase in order to test their performance in addition to their possible effects on confinement and transport Chuyanov et al. (2010). It is important to understand their effect on rotation during the early phases of ITER. Experiments at DIII-D using mock-ups of TBMs found a reduction in toroidal rotation by as much as 60% due to an $n=1$ locked mode Schaffer et al. (2011). Here $n$ is the toroidal mode number. Further experiments showed compensation by $n=1$ control coils may enable access to low NBI torque (1.1 Nm) regimes relevant to ITER without rotation collapse Lanctot et al. (2017). In addition to TBMs, ferritic steel plates (FIs) will be installed in each of the TF coil sections in order to mitigate energetic particle loss due to TF ripple Tobita et al. (2003). Experiments including FIs on JT-60U Urano et al. (2007) and JFT-2M Kawashima et al. (2001) have found a reduction in counter-current rotation with the addition of FIs. As FIs will decrease TF ripple, they may decrease the NTV in ITER.

While the bounce-averaged radial drift vanishes in a tokamak, trapped particles may wander off the flux surface in the presence of non-axisymmetric error fields. Particles trapped poloidally (bananas) can drift radially as the parallel adiabatic invariant, $J_{||}=\oint dl\,v_{||}$, becomes a function of toroidal angle in broken symmetry. Here $v_{||}$ is the velocity coordinate parallel to $\bm{b}=\bm{B}/B$ and integration is taken along the field between bounce points. If local ripple wells exist along a field line and the collisionality is small enough that helically trapped particles can complete their collisionless orbits, these trapped particles may grad-$B$ drift away from the flux surface Stringer (1972). The TF ripple in ITER causes local wells along the field line, corresponding to $\alpha=\epsilon/(qn\delta_{B})<1$ Stringer (1972). Here $\epsilon=r/R$ is the inverse aspect ratio, $r$ is the minor radius, $R$ is the major radius, $q$ is the safety factor, and $\delta_{B}$ is a measure of the amplitude of the ripple. Because of ITER’s low collisionality, $\nu_{*}\ll(\delta_{B}/\epsilon)^{3/2}$, ripple-trapped particles can complete their collisionless orbits Shaing (2003). Here the normalized collision frequency is $\nu_{*}=qRv_{ti}/(\nu_{ii}\epsilon^{3/2})$ where the ion-ion collision frequency is $\nu_{ii}$. The ion thermal velocity is $v_{ti}=\sqrt{2T_{i}/m_{i}}$ where $T_{i}$ is the ion temperature and $m_{i}$ is the ion mass. Therefore, both ripple trapping and banana diffusion should be considered for NTV in ITER. For a general electric field, the neoclassical electron and ion fluxes are not necessarily identical in broken symmetry. The resulting radial current induces a $\bm{J}\times\bm{B}$ torque which is often counter-current.

Analytic expressions for neoclassical fluxes in several rippled tokamak regimes have been derived by various authors, making assumptions about the magnitude of the perturbing field, electric field, magnetic geometry, collisionality, and the collision operator. Multiple regimes are typically needed to describe all radial positions, classes of particles, and helicities of the magnetic field for a single discharge. When collisions set the radial step size of trapped particles, the transport scales as $1/\nu$ where $\nu$ is the collision frequency. The $1/\nu$ regime can be relevant for both ripple trapped and banana particles with small radial electric field. With a non-zero radial electric field, transport from the collisional trapped-passing boundary layer leads to fluxes that scale as $\sqrt{\nu}$. When the collisionality is sufficiently low, the collisionless detrapping/trapping layer becomes significant, where fluxes scale as $\nu$. Here bananas can become passing particles due to the variation of $B_{\max}$ along their drift trajectories Shaing et al. (2009a), and ripple trapped particles can experience collisionless detrapping from ripple wells to become bananas Shaing and Callen (1982a, b). If the collisionality is small compared with the typical toroidal precession frequency of trapped particles, the resonant velocity space layer where the bounce-averaged toroidal drift vanishes can dominate the neoclassical fluxes, leading to superbanana-plateau transport Shaing et al. (2009b). In the presence of a strong radial electric field, the resonance between the parallel bounce motion and drift motion of trapped particles can also result in enhanced transport, known as the bounce-harmonic resonance Linsker and Boozer (1982); Park et al. (2009). The $1/\nu$ and $\sqrt{\nu}$ stellarator regimes for helically-trapped particles have been formulated by Galeev and Sagdeev Galeev and Sagdeev (1969), Ho and Kulsrud Ho and Kulsrud (1987), and Frieman Frieman (1970). These results were generalized to rippled tokamaks in the $1/\nu$ regime by Stringer Stringer (1972), Connor and Hastie Connor and Hastie (1973), and Yushmanov Yushmanov (1982). Kadomtsev and Pogutse Kadomtsev and Pogutse (1971) and Stringer Stringer (1972) presented the scaling of ripple diffusion including trapping/detrapping by poloidal rotation, where fluxes scale as $\nu$. This regime is likely to be applicable for ITER’s low collisionality and strong radial electric field. Banana diffusion in the $1/\nu$ regime has been evaluated by Davidson Davidson (1976), Linkser and Boozer Linsker and Boozer (1982), and Tsang Tsang (1977). The corresponding $\nu$ transport was studied by Tsang Tsang (1977) and Linsker and Boozer Linsker and Boozer (1982). Shaing emphasized the relationship between nonaxisymmetric neoclassical transport and toroidal viscosity Shaing and Callen (1983). The theory for NTV torque due to banana diffusion has been formulated in the $1/\nu$ Shaing (2003), $\nu-\sqrt{\nu}$ Shaing et al. (2008), $\nu$ Shaing et al. (2009a), and superbanana-plateau Shaing et al. (2009b) regimes in addition to an approximate analytic formula which connects these regimes Shaing et al. (2010).

The calculation of NTV torque requires two steps: (i) determine the equilibrium magnetic field in the presence of ripple and (ii) solve a drift kinetic equation (DKE) with the magnetohydrodynamic (MHD) equilibrium or apply reduced analytic formulae. The first step can be performed using various levels of approximation. The simplest method is to superimpose 3D ripple vacuum fields on an axisymmetric equilibrium, ignoring the plasma response. A second level of approximation is to use a linearized 3D equilibrium code such as the Ideal Perturbed Equilibrium Code (IPEC) Park et al. (2009) or linear M3D-C1 Jardin et al. (2008). A third level of approximation is to solve nonlinear MHD force balance using a code such as the Variational Moments Equilibrium Code (VMEC) Hirshman et al. (1986a) or M3D-C1 Ferraro et al. (2010) run in nonlinear mode. In this paper we use free-boundary VMEC to find the MHD equilibrium in the presence of TF ripple, FIs, and TBMs.

Many previous NTV calculations Zhu et al. (2006); Hua et al. (2010); Cole et al. (2011); Park et al. (2009) have been performed using reduced analytic models with severe approximations. Solutions of the bounce-averaged kinetic equation have been found to agree with Shaing’s analytic theory except in the transition between regimes Sun et al. (2010). However, the standard bounce-averaged kinetic equation does not include contributions from bounce and transit resonances. Discrepancies have been found between numerical evaluation of NTV using the Monte Carlo neoclassical solver FORTEC-3D and analytic formulae for the $1/\nu$ and superbanana-plateau regimes Satake et al. (2011a, b). NTV calculations with quasilinear NEO-2 differ from Shaing’s connected formulae Shaing et al. (2010), especially in the edge where the large aspect ratio assumption breaks down Martitsch et al. (2016). Rather than applying such reduced models, in this paper a DKE is solved using the Stellarator Fokker-Planck Iterative Neoclassical Conservative Solver (SFINCS) Landreman et al. (2014) to calculate neoclassical particle and heat fluxes for an ITER steady-state scenario. The SFINCS code does not exploit any expansions in collisionality, size of perturbing field, or magnitude of the radial electric field (beyond the assumption of small Mach number). It also allows for realistic experimental magnetic geometry rather than using simplified flux surface shapes. All trapped particle effects including ripple-trapping Stringer (1972), banana diffusion Linsker and Boozer (1982), and bounce-resonance Linsker and Boozer (1982) are accounted for in these calculations. The DKE solved by SFINCS ensures intrinsic ambipolarity for axisymmetric or quasisymmetric flux surfaces in the presence of a radial electric field while this property is not satisfied by other codes such as DKES Hirshman et al. (1986b); van Rij and Hirshman (1989). This prevents spurious NTV torque density, which is proportional to the radial current. As SFINCS makes no assumption about the size of ripple, it can account for non-quasilinear transport, such as ripple trapping, rather than assuming that the Fourier modes of the ripple can be decoupled. For TF ripple, the deviation from the quasilinear assumption has been found to be significant in benchmarks between SFINCS and NEO-2 Martitsch et al. (2016).

In addition to NTV, neutral beams will provide an angular momentum source for ITER. As NBI torque scales as $P/E^{1/2}$ for input power $P$ and particle energy $E$, ITER’s neutral beams, with $E=1$ MeV and $P=33$ MW, will provide less momentum than in other tokamaks such as JET, with $E=125$ keV for $P=34$ MW Ćirić et al. (2011). NBI-driven rotation will also be smaller in ITER because of its relatively large moment of inertia, with $R=6$ m compared to 3 m for JET. However, spontaneous rotation may be significant in ITER. Turbulence can drive significant flows in the absence of external momentum injection, known as intrinsic or spontaneous rotation. This can be understood as a turbulent redistribution of toroidal angular momentum to produce large directed flows. For perturbed tokamaks this must be in the approximate symmetry direction. According to gyrokinetic orderings and inter-machine comparisons by Parra et al Parra et al. (2012), intrinsic toroidal rotation is expected to scale as $V_{\zeta}\sim T_{i}/I_{p}$ where $I_{p}$ is the plasma current, and core rotations may be on the order of 100 km/s (ion sonic Mach number $M_{i}\approx 8\%$) in ITER. Scalings with $\beta_{N}=\beta_{T}aB_{T}/I_{P}$ by Rice et al Rice et al. (2007) predict rotations of a slightly larger scale, $V_{\zeta}\approx 400$ km/s ($M_{i}\approx 30\%$). Here $\beta_{T}=2\mu_{0}p/B_{T}^{2}$, $B_{T}$ is the toroidal magnetic field in tesla, $a$ is the minor radius at the edge in meters, and $p$ is the plasma pressure. Co-current toroidal rotation appears to be a common feature of H-mode plasmas and has been observed in electron cyclotron (EC) DeGrassie et al. (2007), ohmic DeGrassie et al. (2007), and ion cyclotron range of frequencies (ICRF) Noterdaeme et al. (2003) heated plasmas. Gyrokinetic GS2 simulations with H-mode parameters find an inward intrinsic momentum flux, corresponding to a rotation profile peaked in the core toward the co-current direction Lee et al. (2014). In an up-down symmetric tokamak, the radial intrinsic angular momentum flux can be shown to vanish to lowest order in $\rho_{*}=\rho_{i}/a$, but neoclassical departures from an equilibrium Maxwellian can break this symmetry and cause non-zero rotation in the absence of input momentum Barnes et al. (2013). Here $\rho_{i}=v_{ti}m_{i}/{Z_{i}eB}$ is the gyroradius where $Z_{i}$ is the ion species charge.

In section II we present the ITER steady state scenario and free boundary MHD equilibrium in the presence of field ripple. In section III we estimate rotation driven by NBI and turbulence. This flow velocity is related to $E_{r}$ in section IV. The NTV torque due to TF ripple, TBMs, and FIs is evaluated in section V. In section VI the scaling of transport calculated with SFINCS with ripple magnitude is compared with that predicted by NTV theory, and in section VII neoclassical heat fluxes in the presence of ripple are presented. In section VIII, we assess several tangential magnetic drift models on the transport for this ITER scenario and a radial torque profile is presented. In section IX we summarize the results and conclude.

II ITER Steady State Scenario and Free Boundary Equilibrium Calculations

We consider an advanced ITER steady state scenario with significant bootstrap current and reversed magnetic shear Poli et al. (2014). The input power includes 33 MW NBI, 20 MW EC, and 20 MW lower hybrid (LH) heating for a fusion gain of $Q=5$. This 9 MA non-inductive scenario is achieved with operation close to the Greenwald density limit. The discharge was simulated using the Tokamak Simulation Code (TSC) in the IPS Elwasif et al. (2010) framework for the calculation of the free-boundary equilibrium and the RF calculations, and TRANSP for calculations of the NBI heating and current and torque. The discharge was simulated using the Tokamak Simulation Code (TSC) Jardin et al. (1986) and TRANSP Hawryluk (1979) using a current diffusive ballooning mode (CDBM) Fukuyama et al. (1995, 1998) transport model and EPED1 Snyder et al. (2011) pedestal modeling. The NBI source is modeled using NUBEAM Goldston et al. (1981); Pankin et al. (2004) with 1 MeV particles. The beams are steered with one on-axis and one off-axis, which avoids heating on the midplane wall gap and excess heat deposition above or below the midplane. Further details of the steady state scenario modeling can be found in table 1 of Poli et al. (2014).

The density ($n$), temperature ($T$), safety factor ($q$), total plasma current, and bootstrap current profiles are shown in figure 1. Neoclassical transport will be analyzed in detail at the radial locations indicated by dashed horizontal lines ($r/a=0.5,0.7,0.9$). Throughout we will use the radial coordinate $r/a\propto\sqrt{\Psi_{\mathrm{T}}}$ where $\Psi_{\mathrm{T}}$ is the toroidal flux.

The VMEC free boundary Hirshman et al. (1986a) magnetic equilibrium was computed using the TRANSP profiles along with filamentary models of the toroidal field (TF), poloidal field (PF), and central solenoid (CS) coils and their corresponding currents. The vacuum fields produced by the three TBMs and the FIs have been modeled using FEMAG Shinohara et al. (2009). The equilibrium is computed for four geometries: (i) including only the TF ripple, (ii) including TF ripple, TBMs, and FIs, (iii) TF ripple and FIs, and (iv) axisymmetric geometry. We define the magnitude of the magnetic field ripple to be

[TABLE]

where the maximum and minimum are evaluated at fixed radius and VMEC poloidal angle $\theta$. In figure 2, $\delta_{B}$ is plotted on the poloidal plane for the three rippled VMEC equilibria. A fourth case is also shown in which the component of $\bm{B}$ with $\mathinner{\!\left\lvert n\right\rvert}=18$ was removed from the geometry with TBMs and FIs in order to consider the $\,\mathinner{\!\left\lvert n\right\rvert}<18$ ripple from the TBMs (bottom right). When only TF ripple is present, significant ripple persists over the entire outboard side, while in the configurations with FIs the ripple is much more localized in $\theta$. When TBMs are present, the ripple is higher in magnitude near the outboard midplane ($\delta_{B}\approx 1.4\%$), while in the other magnetic configurations $\delta_{B}\approx$ 1% near the outboard midplane. For comparison, the TF ripple during standard operations is $0.08\%$ in JET de Vries et al. (2008) and $0.6\%$ in ASDEX Upgrade Martitsch et al. (2016). In JT-60U the amplitude of TF ripple is reduced from $\delta_{B}\approx 1.7\%$ to $\delta_{B}\approx 1\%$ by FIs Urano et al. (2007).

In figure 3, the magnitude of $\bm{B}$ is plotted as a function of toroidal angle $\zeta$ at $\theta=0$ and $\theta=\pi/4$. Away from the midplane ($\theta=\pi/4$) the FIs greatly decrease the magnitude of the TF ripple. Near the midplane the FIs do not decrease the magnitude of the toroidal ripple as strongly, as the number of steel plates is reduced near the midplane Shinohara et al. (2009). The ferromagnetic steel of the TBMs concentrates magnetic flux and locally decreases $B$ in the plasma near their location. This causes enhancement of $\delta_{B}$ near $\theta=0$.

III Estimating Toroidal Rotation

In order to predict the ripple transport in ITER, the radial electric field, $E_{r}=-\Phi^{\prime}(r)$, must be estimated, as particle and heat fluxes are nonlinear functions of $E_{r}$. This is equivalent to predicting the parallel flow velocity, $V_{||}$, which scales monotonically with $E_{r}$. As we simply wish to determine a plausible value of $E_{r}$, the difference between $V_{||}$ and $V_{\zeta}$, the toroidal flow, will be unimportant for our estimates. We define $V_{\zeta}$ in terms of the toroidal rotation frequency, $V_{\zeta}=\Omega_{\zeta}R$, where $\Omega_{\zeta}\approx\Omega_{\zeta}(r)$. As $I_{P}$ and the toroidal magnetic field are both directed clockwise when viewed from above, $V_{\zeta}$ and $V_{||}$ will point in the same direction. Here we use the convention that positive $V_{\zeta}$ corresponds to co-current rotation. For this rotation calculation, angular momentum transport due to neutral beams and turbulence will be considered. There is an additional torque caused by the radial current of orbit-lost alphas Rosenbluth and Hinton (1996), but it will be negligible ($\approx 0.006$ Nm/m3). The following time-independent momentum balance equation is considered in determining $\Omega_{\zeta}(r)$,

[TABLE]

where $\Pi^{\mathrm{turb}}_{\zeta}$ and $\Pi^{\mathrm{NC}}_{\zeta}$ are the toroidal angular momentum flux densities due to turbulent and neoclassical transport and $\tau^{\mathrm{NBI}}$ is the NBI torque density. For this paper the feedback of $\Pi_{\zeta}^{\mathrm{NC}}$ on $\Omega_{\zeta}$ will not be calculated. Determining the change in rotation due to NTV would require iteratively solving this equation for $\Omega_{\zeta}$, as $\Pi_{\zeta}^{\mathrm{NC}}$ is a nonlinear function of $\Omega_{\zeta}$.

The quantity $\Pi_{\zeta}^{\mathrm{turb}}$ consists of a diffusive term as well as a term independent of $\Omega_{\zeta}$ which accounts for turbulent intrinsic rotation,

[TABLE]

For simplicity, an angular momentum pinch, $P_{\zeta}$, will not be considered for this analysis. As $RP_{\zeta}/\chi_{\zeta}\approx 2$, there would be a factor of 2 difference in rotation peaking at the core due to the turbulent momentum source at the edge Lee (2013). Here $\chi_{\zeta}$ is the toroidal ion angular momentum diffusivity. The flux surface average is denoted by $\langle...\rangle$,

[TABLE]

where $\sqrt{g}$ is the Jacobian. Ignoring NTV torque, we will solve the following angular momentum balance equation,

[TABLE]

Equation 6 is a linear inhomogeneous equation for $\Omega_{\zeta}$, as the right hand side is independent of $\Omega_{\zeta}$. We can therefore solve for the rotation due to each of the source terms individually and add the results to obtain the rotation due to both NBI torque and turbulent intrinsic torque.

The NBI-driven rotation profile is evolved by TRANSP assuming $\chi_{\zeta}=\chi_{i}$, the ion heat diffusivity. The total beam torque density, $\tau^{\mathrm{NBI}}$, is calculated by NUBEAM including collisional, $\bm{J}\times\bm{B}$, thermalization, and recombination torques. The following momentum balance equation is solved to compute $\Omega_{\zeta}$ driven by NBI,

[TABLE]

We consider a semi-analytic intrinsic rotation model to determine the turbulent-driven rotation Hillesheim et al. (2015),

[TABLE]

where $\rho_{*,\theta}=v_{ti}m_{i}/(Z_{i}eB_{\theta}\langle R\rangle)$ is the poloidal normalized gyroradius, $B_{\theta}=\bm{B}\cdot\partial\bm{r}/\partial\theta$, and $L_{T}=-\left(\partial\ln T_{i}/\partial r\right)^{-1}$ is the temperature gradient scale length. The Prandtl number $P_{r}=\chi_{\zeta}/\chi_{i}$ is again taken to be 1. Equation 8 is obtained assuming that $\Pi_{\mathrm{int}}$ balances turbulent momentum diffusion in steady state, $\Pi_{\mathrm{int}}=m_{i}n_{i}\chi_{\zeta}\langle R^{2}\rangle\partial\Omega_{\zeta}/\partial r$. This model considers the intrinsic torque driven by the neoclassical diamagnetic flows, such that $\Omega_{\zeta}\sim\rho_{*,\theta}v_{ti}/\langle R\rangle$ and $\Omega_{\zeta}\Pi_{\mathrm{int}}/Q_{i}\sim\rho_{*,\theta}$ where $Q_{i}=n_{i}T_{i}\chi_{i}L_{T}$ is the turbulent energy flux. We also take $\Omega_{\zeta}(a)=0$. The quantity $\widetilde{\Pi}(\nu_{*})$ is an order unity function which characterizes the collisionality dependence of rotation reversals, determined from gyrokinetic turbulence simulations Barnes et al. (2013),

[TABLE]

where $\nu_{c}=1.7$. Because of ITER’s low collisionality, we do not expect a rotation reversal, which is correlated with transitioning between the banana and plateau regimes. Equation 8 was integrated using profiles for the ITER steady state scenario.

The flux-surface averaged toroidal rotation, $\langle V_{\zeta}\rangle=\Omega_{\zeta}(r)\langle R\rangle$, predicted by these models is shown in figure 4. NBI torque contributes to significant rotation at $r/a\lesssim 0.4$ where the torque density also peaks (see figure 13), while turbulent torque produces rotation in the pedestal due to the $L_{T}^{-2}$ scaling of our model. The intrinsic rotation calculated is comparable to that predicted from theoretical scaling arguments by Parra et al Parra et al. (2012), $V_{\zeta}\approx 100$ km/s. At the radii that will be considered for neoclassical calculations (indicated by dashed vertical lines), intrinsic turbulent rotation may dominate over that due to NBI. However, we emphasize that it is an estimate based on scaling arguments, as much uncertainty is inherent in predicting turbulent rotation. The volume-averaged toroidal rotation due to both NBI and turbulent torques, 113 km/s, is slightly larger than that predicted from dimensionless parameter scans on DIII-D, 87 km/s Chrystal et al. (2017).

For stabilization of the resistive wall mode (RWM) in ITER, it has been estimated Liu et al. (2004) that a critical central Mach number, $M_{A}=\Omega_{\zeta}(0)/\omega_{A}\gtrsim 5\%$, must be achieved given a peaked rotation profile. Here $\omega_{A}=B/(\langle R\rangle\sqrt{\mu_{0}m_{i}n_{i}})$ is the Alfvèn frequency. With a central rotation frequency $\Omega_{\zeta}(0)\approx 2\%\,\omega_{A}$ as shown in figure 4, it may be difficult to suppress the RWM in ITER with rotation alone. As this calculation does not take into account NTV torque, $M_{A}$ is likely to be smaller than what is shown. Additionally, the TBM are known to increase the critical rotation frequency as they have a much shorter resistive time scale than the wall Liu et al. (2004). More recent analysis has shown that even above such a critical rotation value, the plasma can become unstable due to resonances between the drift frequency and bounce frequency Berkery et al. (2010); Liu et al. (2009).

IV Relationship Between $E_{r}$ and $V_{||}$

Neoclassical theory predicts a specific linear-plus-offset relationship between $V_{||}$ and $E_{r}$, but it does not predict a particular value for either $V_{||}$ or $E_{r}$ in a tokamak. Neoclassical calculations of $V_{||}$ are made in order to determine an $E_{r}$ profile consistent with our estimate of $V_{\zeta}\approx\langle V_{||}B\rangle/\langle B^{2}\rangle^{1/2}$ made in section III. The parallel flow velocity for species $a$ is computed from the neoclassical distribution function,

[TABLE]

which we calculate with the SFINCS Landreman et al. (2014) code. SFINCS is used to solve a radially-local DKE for the gyro-averaged distribution function, $f_{a1}$, on a single flux surface including coupling between species.

[TABLE]

Here $a$ indicates species, $f_{a0}$ is an equilibrium Maxwellian, $\psi=\Psi_{\mathrm{T}}/2\pi$, $Z$ indicates charge, and $C$ is the linearized Fokker-Planck collision operator. Gradients are performed at constant $W=m_{a}v^{2}/2+Z_{a}e\Phi$ and $\mu=v_{\perp}^{2}/(2B)$. The $\bm{E}\times\bm{B}$ drift is

[TABLE]

and the radial magnetic drift is

[TABLE]

where $v_{\perp}$ is the velocity coordinate perpendicular to $\bm{b}$. The quantity $\Omega_{a}=Z_{a}eB/m_{a}$ is the gyrofrequency. Transport quantities have been calculated using the steady state scenario ion and electron profiles and VMEC geometry. We consider a three species plasma (D, T, and electrons), and we assume that $n_{D}=n_{T}=n_{e}/2$. The second term on the right hand side of equation 11 proportional to $E_{||}$ is negligible for this non-inductive scenario with loop voltage $\approx 10^{-4}$ V. For the calculations presented in sections IV, V, VI, and VII, $\bm{v}_{\mathrm{m}a}\cdot\nabla f_{a1}$ is not included. The effect of keeping this term is shown to be small in section VIII.

The relationship between $E_{r}$ and $\langle V_{||}B\rangle/\langle B^{2}\rangle^{1/2}$ for electrons and ions at $r/a=0.9$ is shown in figure 5. Only one curve is shown for each species as the addition of ripple fields does not change the dependence of $V_{||}$ on $E_{r}$ significantly ($\lesssim 5\%$). While radial transport of heat and particles changes substantially in the presence of small ripple fields (see sections V, VI, and VII), the parallel flow is much less sensitive to the perturbing field. Note that the parallel flow is non-zero in axisymmetry while the radial current vanishes without symmetry-breaking.

In a tokamak we can write $\langle V_{||}^{a}B\rangle$ in terms of a dimensionless parallel flow coefficient, $k_{||}$, and thermodynamic drives,

[TABLE]

where $2\pi\psi_{P}$ is the poloidal flux, $G(\psi)=RB_{\zeta}$, and $B_{\zeta}=\bm{B}\cdot\partial\bm{r}/\partial\zeta$. The low collisionality, large aspect ratio limit Hinton and Hazeltine (1976); Hirshman and Sigmar (1981) $k_{||}\approx 1.17$ is often assumed in NTV theory Callen (2011); Sun et al. (2011) in relating analytic expressions of torque density to toroidal rotation frequency. The value of the ion $k_{||}$ calculated by SFINCS for ITER parameters varies between 0.5 near the edge and 0.9 near the core. The bootstrap current computed with SFINCS,

[TABLE]

is consistent with that computed by TRANSP within 10% for $r/a\geq 0.5$. Though there is some discrepancy in the core, they have the same qualitative behavior and similar maxima. The bootstrap current in TRANSP is computed using a Sauter model Sauter et al. (1999), an analytic fit to numerical solutions of the Fokker-Planck equation.

V Torque Calculation

The NTV torque density, $\tau^{\mathrm{NTV}}$, is calculated from radial particle fluxes, $\Gamma_{\psi}$,

[TABLE]

using the flux-force relation,

[TABLE]

where $B^{\theta}=\bm{B}\cdot\nabla\theta$ and the summation is performed over species. This expression relates radial particle transport to a toroidal angular momentum source caused by the non-axisymmetric field. This relationship can be derived from action-angle coordinates Albert et al. (2016), neoclassical moment equations Shaing (1986), or from the definition of the drift-driven flux Shaing (2006).

The calculation of $\tau^{\mathrm{NTV}}$ for three geometries at $r/a=0.9$ is shown in figure 6. Here positive corresponds to the co-current direction. The numerically computed NTV torque is found to vanish in axisymmetric geometry, as expected. Overall, the magnitude of $\tau^{\mathrm{NTV}}$ with only TF ripple is larger than that with the addition of both the FIs and the TBMs. In figure 8 we show that the $\mathinner{\!\left\lvert n\right\rvert}<18$ TBM ripple produces much less torque than the $\mathinner{\!\left\lvert n\right\rvert}=18$ ripple, so the decrease in $\tau^{\mathrm{NTV}}$ magnitude with both FIs and TBMs can be attributed to the decrease in ripple in the presence of FIs. As will be discussed in section VI, neoclassical ripple transport in most regimes scales positively with $\delta_{B}$. The addition of FIs significantly decreases the magnitude of $\delta_{B}$ across most of the outboard side, and as a result the magnitude of $\tau^{\mathrm{NTV}}$ is reduced. The dashed vertical line indicates the value of $\langle V_{||}B\rangle/\langle B^{2}\rangle^{1/2}$ and $E_{r}$ predicted from the intrinsic and NBI rotation model. At this value of $E_{r}$ the presence of ferritic components decreases the magnitude of the torque density by about $75\%$.

The circle indicates the offset rotation at the ambipolar $E_{r}$. If no other angular momentum source were present in the system, $\tau^{\mathrm{NTV}}$ would drive the plasma to rotate at this velocity. Although $\tau^{\mathrm{NTV}}$ differs significantly between the two geometries they have similar offset rotation velocities, $V_{\zeta}$ = -10 km/s with TF ripple only and -6 km/s with TBMs and FIs. Note that for $E_{r}$ greater than this ambipolar value, $\tau^{\mathrm{NTV}}$ is counter-current while neutral beams and turbulence drive rotation in the co-current direction, so $\tau^{\mathrm{NTV}}$ is a damping torque. The NTV torque due to TF ripple only is larger in magnitude than $\tau^{\mathrm{NBI}}$ and $\tau^{\mathrm{turb}}$ while that with TBMs and FIs is of similar magnitude (see figure 13). Therefore, NTV torque may be key in determining the edge rotation in ITER.

The magnitude of $\tau^{\mathrm{NTV}}$ peaks at $E_{r}=0$ where $1/\nu$ transport becomes dominant. Although $\nu_{*}$ is sufficiently small such that the superbanana-plateau regime becomes relevant, the physics of superbanana formation is not accounted for in these SFINCS calculations which do not include $\bm{v}_{\mathrm{m}}\cdot\nabla f_{1}$. Superbanana-plateau transport will be considered in section VIII. At $r/a=0.9$, the $1/\nu$ regime applies for $\mathinner{\!\left\lvert E_{r}\right\rvert}\lesssim 0.2$ kV/m where the effective collision frequency of trapped particles is larger than the $E\times B$ precession frequency. The peak at small $\mathinner{\!\left\lvert E_{r}\right\rvert}$ also corresponds to the region of $1/\nu$ transport of particles trapped in local ripple wells. Much NTV literature is based on banana diffusion and ripple trapping in the $1/\nu$ regime Stringer (1972); Connor and Hastie (1974), which is not applicable for the range of $E_{r}$ predicted for ITER. For the range of applicable $E_{r}$, bounce-harmonic resonance may occur. The $l=1$, $n=18$ resonance condition, $\omega_{b}-n(\omega_{E}+\omega_{B})=0$, will be satisfied for $v_{||}\approx v_{ti}$ at $E_{r}\approx 7$ kV/m. Here $\omega_{b}$ is the bounce frequency, $\omega_{E}$ is the $E\times B$ precession frequency, and $\omega_{B}$ is the toroidal magnetic drift precession Park et al. (2009). Note that here $\bm{v}_{\text{m}}\cdot\nabla f_{1}$ is not included in the kinetic equation ($\omega_{B}=0$), but the physics of the bounce harmonic resonance between $\omega_{E}$ and $\omega_{b}$ is still accounted for in our calculation. However, we see no evidence of enhanced $\tau^{\mathrm{NTV}}$ near this $E_{r}$ that would be indicative of a bounce-harmonic resonance.

NTV torque is often expressed in terms of a toroidal damping frequency, $\nu_{\zeta}$,

[TABLE]

where $\Omega_{\zeta,\mathrm{offset}}$ is the offset rotation frequency. We note that $\tau^{\mathrm{NTV}}$ does appear to scale linearly with $E_{r}$ (and thus $\Omega_{\zeta}$) for $\mathinner{\!\left\lvert E_{r}\right\rvert}\gtrsim 30$ kV/m. However, $\tau^{\mathrm{NTV}}$ is a complicated nonlinear function of $\Omega_{\zeta}$ for $\mathinner{\!\left\lvert E_{r}\right\rvert}\lesssim 30$ kV/m at the transition between collision-limited $1/\nu$ transport and $\nu-\sqrt{\nu}$ transport, so equation 18 is not a very useful representation in this context.

In figure 7 we present $\tau^{\mathrm{NTV}}$ at $r/a=0.9$ due to the electron and ion radial current in the presence of TF ripple only (left) and TF ripple with ferromagnetic components (right). The $E_{r}$ corresponding to the offset rotation frequency for the electrons is positive while that of the ions is negative. At the predicted $E_{r}$, $\tau^{\mathrm{NTV}}$ due to the electron particle flux is positive while that due to ion particle flux is negative. At all radial locations the electron contribution to $\tau^{\mathrm{NTV}}$ is less than 10% of the total torque density.

In order to decouple the influence of the FI ripple and the TBM ripple, $\tau^{\mathrm{NTV}}$ at $r/a=0.9$ is calculated for toroidal modes (i) $\mathinner{\!\left\lvert n\right\rvert}\leq 18$, (ii) $\mathinner{\!\left\lvert n\right\rvert}=18$, and (iii) $\mathinner{\!\left\lvert n\right\rvert}<18$, shown in figure 8. For $\mathinner{\!\left\lvert n\right\rvert}\leq 18$ and $\mathinner{\!\left\lvert n\right\rvert}=18$, VMEC free boundary equilibria were computed including these toroidal modes. For $\mathinner{\!\left\lvert n\right\rvert}<18$, the SFINCS calculation was performed including the desired $n$ from the VMEC fields. Here $B$ is decomposed as,

[TABLE]

where $\theta$ and $\zeta$ are VMEC angles. The covariant and contravariant components of $B$ along with their partial derivatives and $\sqrt{g}$ are similarly decomposed such that the DKE can be solved for the desired toroidal modes.

The TBM produces a wide spectrum of toroidal perturbations, including $\mathinner{\!\left\lvert n\right\rvert}=1$ and $\mathinner{\!\left\lvert n\right\rvert}=18$. While the FIs decrease the magnitude of the $\mathinner{\!\left\lvert n\right\rvert}=18$ ripple, the TBM contributes most strongly to low mode numbers. As SFINCS is not linearized in the perturbing field, the torque due to $\mathinner{\!\left\lvert n\right\rvert}\leq 18$ is the not the sum of the torques due to $\mathinner{\!\left\lvert n\right\rvert}=18$ and $\mathinner{\!\left\lvert n\right\rvert}<18$. We find that the $\mathinner{\!\left\lvert n\right\rvert}=18$ ripple drives about 100 times more torque than the lower $n$ ripple. This result is in agreement with most relevant rippled tokamak transport regimes, which feature positive scaling with $n$ Shaing (2003); Shaing et al. (2008). For tokamak banana diffusion, in the $\sqrt{\nu}$ boundary layer Shaing et al. (2008) ion transport scales as $\Gamma_{\psi}\sim\sqrt{n}$ and in the $1/\nu$ regime Shaing (2003) $\Gamma_{\psi}\sim n^{2}$. Moreover, it is more difficult to form ripple wells along a field line from low-$n$ ripple, so ripple trapping cannot contribute as strongly to transport. This matches our findings that the higher harmonic $\mathinner{\!\left\lvert n\right\rvert}=18$ ripple contributes more strongly to $\tau^{\mathrm{NTV}}$ than the $\mathinner{\!\left\lvert n\right\rvert}<18$ ripple.

In figure 9, the SFINCS calculation of $\tau^{\mathrm{NTV}}$ with TF ripple only is shown at $r/a$ = 0.5, 0.7, and 0.9. For these three radii the maximum $\delta_{B}=0.26\%$, 0.51%, and 0.82% respectively. As $\tau^{\mathrm{NTV}}$ scales with a positive power of $\delta_{B}$ in most rippled tokamak regimes, it is reasonable to expect that the magnitude of $\tau^{\mathrm{NTV}}$ would decrease with decreasing radius. On the other hand, transport scales strongly with $T_{i}$. In the $\sqrt{\nu}$ banana diffusion regime Shaing et al. (2008) $\Gamma_{\psi}\sim v_{ti}^{4}\sqrt{\nu_{ii}}\sim T_{i}^{5/4}$. The combined effect of decreased ripple and increased temperature with decreasing radius leads to comparable torques with decreasing radius in the presence of significant $E_{r}$. The scaling with $T_{i}$ is even stronger in the $1/\nu$ regime Stringer (1972); Shaing (2003), where $\Gamma_{\psi}\sim v_{ti}^{4}/\nu_{ii}\sim T_{i}^{7/2}$. Indeed, we find that the magnitude of $\tau^{\mathrm{NTV}}$ at $E_{r}=0$ increases with decreasing radius.

VI Scaling with Ripple Magnitude

In figure 10, the NTV torque density calculated by SFINCS is shown as a function of the magnitude of the ripple, $\delta_{B}$, for TF only geometry. The additional ferromagnetic ripple is not included, while the $\mathinner{\!\left\lvert n\right\rvert}=18$ components of $\bm{B}$, its derivatives, and $\sqrt{g}$ are rescaled as described above. The quantity $\tau^{\mathrm{NTV}}$ is calculated at $r/a=0.9$ with $E_{r}=30$ kV/m, corresponding to the intrinsic rotation estimate. The color-shaded background indicates the approximate regions of applicability of the collisional boundary layer ($\nu-\sqrt{\nu}$) and the collisionless detrapping/trapping ($\nu$) regimes. The boundary between these regimes corresponds to the $\delta_{B}$ for which the width in pitch angle of the detrapping/retrapping layer is similar to the width of the collisional boundary layer, $(\delta_{B}/\epsilon)\sim(\nu/(\epsilon\omega_{E}))^{1/2}$. The $1/\nu$ regime Shaing (2003) does not apply at this $E_{r}$, as $\omega_{E}\gg\nu/\epsilon$ where $\omega_{E}=E_{r}/B^{\theta}$ is the $E\times B$ precession frequency. The radial electric field is also large enough that the resonance between $\bm{v}_{E}$ and $\bm{v}_{\mathrm{m}}$ cannot occur, so the superbanana-plateau Shaing et al. (2009b) and superbanana Shaing et al. (2009c) regimes are avoided. This significant $E_{r}$ may also allow the bounce-harmonic resonance to occur Park et al. (2009). Transport from ripple-trapped particles in the $\nu-\sqrt{\nu}$ regime may also be significant for these parameters.

The observed scaling appears somewhat consistent with ripple trapping in the stellarator $\sqrt{\nu}$ regime Ho and Kulsrud (1987) which predicts $\Gamma_{\psi}\sim\delta_{B}^{3/2}$. However, this result is inconsistent with predictions for tokamak ripple transport in the $\nu$ regime, $\Gamma_{\psi}\sim\delta_{B}^{0}$ Tsang (1977); Linsker and Boozer (1982). Contributions from other transport regimes may also influence the observed scaling. In the banana diffusion $\sqrt{\nu}$ regime $\tau^{\mathrm{NTV}}\sim\delta_{B}^{2}$ and in the $\nu$ regime $\tau^{\mathrm{NTV}}\sim\delta_{B}$. Bounce-harmonic resonant fluxes scale as $\delta_{B}^{2}$ Park et al. (2009). A scaling between $\delta_{B}^{0}$ and $\delta_{B}^{2}$ has been predicted for plasmas close to symmetry with large gradient ripple in the absence of $E_{r}$ Calvo et al. (2014). For $\delta_{B}$ smaller than $0.82\%$, the actual value of ripple at $r/a=0.9$ for ITER geometry, the scaling of $\tau^{\mathrm{NTV}}$ with $\delta_{B}$ appears similar to $\delta_{B}^{3/2}$. The disagreement between the SFINCS calculations and the quasilinear prediction, $\Gamma_{\psi}\sim\delta_{B}^{2}$, indicates the presence of nonlinear effects such as local ripple trapping and collisionless detrapping. The departure from quasilinear scaling increases with $\delta_{B}$, which is consistent with comparisons of SFINCS with quasilinear NEO-2 Martitsch et al. (2016). We see that $\tau^{\text{NTV}}$ shows very shallow scaling between $\delta_{B}=0.05$ and $\delta_{B}=0.2$. This could be in agreement with a scaling of $\delta_{B}^{0}$ predicted for $\nu$ regime ripple-trapping in tokamaks Kadomtsev and Pogutse (1971); Stringer (1972). In this region the collisionless detrapping boundary layer and collisional boundary layer are of comparable widths, so it is possible that the transport here is not described well by any of the displayed scalings. Furthermore, near the collisionless detrapping-trapping regime, $\delta_{B}$ becomes comparable to the inverse aspect ratio and the assumptions made for rippled tokamak theory are not satisfied.

VII Heat Flux Calculation

As well as driving non-ambipolar particle fluxes, the breaking of toroidal symmetry drives an additional neoclassical heat flux. In figure 11, the SFINCS calculation of the heat flux, $Q^{\mathrm{NC}}$, is shown for three magnetic geometries: (i) axisymmetric (blue solid), (ii) with TF ripple only (red dash-dot), and (iii) TF ripple with TBMs and FIs (green dashed). In the presence of TF ripple, the ripple drives an additional heat flux that is comparable to the axisymmetric heat flux. However, with the addition of the FIs the heat flux is reduced to the magnitude of the axisymmetric value, except near $E_{r}=0$ where $1/\nu$ transport dominates.

While the radial ripple-driven particle fluxes will significantly alter the ITER angular momentum transport, the neoclassical heat fluxes are insignificant in comparison to the turbulent heat flux. Note that the neoclassical heat flux is $\lesssim 5\%$ of the heat flux calculated from heating and fusion rate profiles (see appendix A), $Q\approx 0.2$ MW/m2. Thus we can attribute $\gtrsim 95\%$ of the heat transport to turbulence. If ITER ripple were scaled up to $\delta_{B}\gtrsim 30\%$, the neoclassical ripple heat transport would be comparable to the anomalous transport at this radius.

VIII Tangential Magnetic Drifts

Although $(\bm{v}_{E}+\bm{v}_{\mathrm{m}})\cdot\nabla f_{1}$ is formally of higher order than the other terms in equation 11, it has been found to be important when $\nu_{*}\lesssim\rho_{*}$ Calvo et al. (2017); Matsuoka et al. (2015) and has been included in other calculations of 3D neoclassical transport. In the SFINCS calculations shown in sections IV, V, VI, and VII, $\bm{v}_{\mathrm{m}}\cdot\nabla f_{1}$ has not been included, but now we examine the effect of including parts of this term. As SFINCS does not maintain radial coupling of $f_{1}$, only the poloidal and toroidal components of this magnetic drift term can be retained while the radial component cannot. Note that the radial magnetic drift is retained in $\bm{v}_{\mathrm{m}}\cdot\nabla f_{0}$. We first implement $\bm{v}_{\mathrm{m}}\cdot\nabla\theta$ and $\bm{v}_{\mathrm{m}}\cdot\nabla\zeta$ using the following form of the magnetic drifts,

[TABLE]

where $\xi=v_{||}/v$. However, a coordinate-dependence can be introduced as we simply drop one component of $\bm{v}_{\mathrm{m}}$. For a coordinate-independent form, one must project $\bm{v}_{\mathrm{m}}$ onto the flux surface. Additionally, when poloidal and toroidal drifts are retained, the effective particle trajectories do not necessarily conserve $\mu$ when $\mu=0$. The drifts can be regularized in order to satisfy $\dot{\xi}(\xi=\pm 1)=0$. Regularization also eliminates the need for additional particle and heat sources due to the radially local assumption and preserves ambipolarity of axisymmetric systems Sugama et al. (2016). To this end, we also implement a coordinate-independent magnetic drift perpendicular to $\nabla\psi$,

[TABLE]

Note that the $\nabla B$ drift term is regularized while the curvature drift term is not. As tangential drifts are important for the trapped portion of velocity space, we can consider $\xi^{2}\ll 1$. For this reason we drop the curvature drift for regularization,

[TABLE]

This is similar to the form presented by Sugama Sugama et al. (2016), but we have chosen a different form of regularization. This choice for $\bm{v}_{\mathrm{m}a}^{\perp}$ does not alter the conservation properties shown by Sugama, as it remains in the $\bm{B}\times\nabla\psi$ direction and vanishes at $\xi=\pm 1$. We note that the phase space conservation properties rely on the choice of a modified Jacobian in the presence of tangential magnetic drifts. In SFINCS we have not implemented such a modification. However, as Sugama shows, the correction to the Jacobian is an order $\rho_{*}$ correction. As particle and heat sources have been implemented in SFINCS, we have confirmed that the addition of tangential magnetic drifts does not necessitate the use of appreciable source terms. We note that this form of the tangential magnetic drifts we have chosen does not include a magnetic shear term which is present in the bounce-averaged radial drift. This non-local modification has been found to significantly alter superbanana transport Albert et al. (2016); Shaing (2015) and the drift-orbit resonance Martitsch et al. (2016).

An $E_{r}$ scan at $r/a=0.7$, where $\rho_{*}$ becomes comparable to $\nu_{*}$, is shown in figure 12. When $\bm{v}_{\mathrm{m}}\cdot\nabla f_{1}$ is added to the kinetic equation, the $1/\nu$ peak at $E_{r}=0$ is shifted toward a slightly negative $E_{r}$, corresponding to the region where $(\bm{v}_{E}+\bm{v}_{\mathrm{M}})\cdot\nabla\zeta\approx 0$, where superbanana-plateau transport takes place. For ITER parameters at this radius, the collisionality is large enough that superbananas cannot complete their collisionless trajectories but small enough that non-resonant trapped particles precess, $\nu_{*}^{\mathrm{SB}}\ll\nu_{*}\ll\nu_{*}^{\mathrm{SBP}}$, where $\nu_{*}^{\mathrm{SBP}}=\rho_{*}q^{2}/\epsilon^{1/2}$ and $\nu_{*}^{\mathrm{SB}}=\rho_{*}\delta_{B}^{3/2}q^{2}/\epsilon^{2}$ Shaing et al. (2009c, b), thus superbanana-plateau transport is relevant.

When the in-surface magnetic drifts are present, the depth of the resonant peak is diminished. In the absence of tangential drifts, the bounce-averaged toroidal drift vanishes at $E_{r}=0$ for all particles regardless of pitch angle and energy. When tangential drifts are added to the DKE, the resonant peak will occur at the $E_{r}$ for which thermal trapped particles satisfy the resonance condition. However, only particles above a certain energy and at the resonant pitch angle will participate in the superbanana-plateau transport, thus the depth of the peak is diminished. Note that local ripple trapping might also contribute to the $1/\nu$ transport at small $\mathinner{\!\left\lvert E_{r}\right\rvert}$. For $\mathinner{\!\left\lvert E_{r}\right\rvert}>20$ kV/m, the range relevant for ITER, the addition of $\bm{v}_{\mathrm{m}}\cdot\nabla f_{1}$ has a negligible effect on $\tau^{\mathrm{NTV}}$. The addition of tangential magnetic drifts would not dramatically change the results in previous sections.

We compute a radial profile of $\tau^{\mathrm{NTV}}$ due to TF ripple including $\bm{v}_{\text{m}}^{\perp}\cdot\nabla f_{1}$. The intrinsic rotation model and NBI rotation model are used to estimate $E_{r}$ at each radius, as shown in the $E_{r}$ profile in figure 13. As $E_{r}$ crosses through 0 at $r/a=0.56$ for the NBI rotation model, tangential drifts will affect the transport. In figure 13, we compare the magnitude of $\tau^{\mathrm{NTV}}$ due to TF ripple with $\tau^{\mathrm{NBI}}$ and $\tau^{\mathrm{turb}}=-\nabla\cdot\Pi_{\mathrm{int}}$, the turbulent momentum source causing intrinsic rotation. The $\tau^{\mathrm{NBI}}$ profile was computed by NUBEAM as used in section III, and $\tau^{\mathrm{turb}}$ is estimated using $\Pi_{\mathrm{int}}\sim(\rho_{\theta}/L_{T})\widetilde{\Pi}(\nu_{*})Q(\langle R\rangle/v_{ti})$ (see appendix A). At $r/a=0.62$ superbanana-plateau transport dominates when NBI rotation is considered, and $\tau^{NTV}$ is about 6 times larger than when the higher-rotation turbulent torque $E_{r}$ is considered. For both $E_{r}$ estimates $\mathinner{\!\left\lvert\tau^{\text{NTV}}\right\rvert}$ increases with decreasing radius due to the scaling with $T_{i}$ as discussed in section V. Note that the turbulent torque produces much rotation in the pedestal according to this model as $\tau^{\mathrm{turb}}\propto 1/L_{T}$. The integrated NTV torque, -45.6 Nm with the turbulent rotation model and -71 Nm with the NBI rotation model, is larger in magnitude than the NBI torque, 35 Nm, but smaller than the turbulent torque, 93 Nm. Here the integrated $\tau^{\text{turb}}$ is significantly larger than that obtained from dimensionless parameter scans on DIII-D, 33 Nm Chrystal et al. (2017). This is possibly due to the assumed scaling in our turbulent rotation model, which may not be physical near the edge.

In the region $0.5\lesssim r/a\lesssim 0.9$, the magnitude of $\tau^{\mathrm{NTV}}$ is comparable to $\tau^{\mathrm{turb}}$ and greater than $\tau^{\mathrm{NBI}}$. The NTV torque will likely significantly damp rotation in the absence of inserts, decreasing MHD stability. However, the resulting rotation profile may be sheared because of the significant counter-current NTV source at the edge and co-current NBI source in the core. We estimate the rotation shear, $\gamma=\Delta V_{\zeta}/\Delta r\approx 0.4(v_{ti}/R)$, using the neoclassical offset at $r/a=0.7$ and the NBI-driven rotation at $r/a=0.4$. Assuming the maximum linear growth rate for drift wave instabilities, $\gamma_{\mathrm{Lin}}\approx v_{ti}/R$ Connor et al. (2004), this rotation shear may be large enough to suppress microturbulence. In concert with reversed magnetic shear sustained by heating and current drive sources Poli et al. (2014), rotation shear may support the formation of an ITB Waltz et al. (1994) for this steady state scenario.

IX Summary

We calculate neoclassical transport in the presence of 3D magnetic fields, including toroidal field ripple and ferromagnetic components, for an ITER steady state scenario. We use an NBI and intrinsic turbulent rotation model to estimate $E_{r}$ for neoclassical calculations. We find that without considering $\tau^{\mathrm{NTV}}$, toroidal rotation with $M_{A}\lesssim 2\%\,$ is to be expected, which is likely not large enough to suppress resistive wall modes Liu et al. (2004). We use VMEC free boundary equilibria in the presence of ripple fields to calculate neoclassical particle and heat fluxes using the drift-kinetic solver, SFINCS. At large radii $r/a\gtrsim 0.5$, $\tau^{\mathrm{NTV}}$ due to TF ripple without ferritic components is comparable to $\tau^{\mathrm{NBI}}$ and $\tau^{\mathrm{turb}}$ in magnitude but opposite in sign, which may result in flow damping at the edge and a decrease in MHD stability. As the integral NTV torque is similar in magnitude to the NBI torque, non-resonant magnetic braking cannot be ignored in analysis of ITER rotation. The torque profile may also result in a significant rotation shear which could suppress turbulent transport. While the addition of FIs significantly reduces the transport ($\approx 75\%$ reduction at $r/a=0.9$), the low $n$ perturbation of the TBM produces very little NTV torque. The neoclassical heat flux caused by ripple is insignificant in comparison to the turbulent heat flux. Though NTV torque has been shown to be important for ITER angular momentum balance, iteratively solving for the rotation profile with $\tau^{\mathrm{NTV}}$ will be left for future consideration.

Several transport regimes must be considered for ITER NTV: the $\nu-\sqrt{\nu}$ banana diffusion, bounce-resonance, and $\nu-\sqrt{\nu}$ ripple trapping regimes. The calculated scaling of $\tau^{\mathrm{NTV}}$ with $\delta_{B}$ is between the $\delta_{B}^{3/2}$ scaling of the ripple trapping $\sqrt{\nu}$ regime and the $\delta_{B}^{2}$ scaling predicted in the $\nu-\sqrt{\nu}$ regime at small $\delta_{B}$. There is room for further comparison between SFINCS calculations of $\tau^{\mathrm{NTV}}$ and analytic fomulae. However, we note that the analytic theory for transport of ripple-trapped particles in a tokamak close to axisymmetry in the presence of $E_{r}$ is not fully developed.

Appendix A Approximate Turbulent Heat Flux and Torque

As $\Pi_{\mathrm{int}}$ is proportional to $Q$ in our model, we must estimate $Q$ using the input heating power and D-T fusion rates calculated with TRANSP. The LH, NBI, and ECH power densities ($P_{\mathrm{LH}}$, $P_{\mathrm{NBI}}$, and $P_{\mathrm{ECH}}$) are integrated along with the fusion reaction rate density ($R_{\mathrm{DT}}$) to calculate the total integrated heating source, $H(r)$,

[TABLE]

As $\int QdS=\int HdV$,

[TABLE]

where $A(r)$ is the flux surface area and $V(r)$ is the volume enclosed by flux surface $r$. We have shown in section VII that the neoclassical heat flux is insignificant in comparison to $Q(r)$, so we can attribute $Q(r)$ to turbulent heat transport. The calculated $Q$ is shown in figure 14.

We estimate $\tau^{\mathrm{turb}}=-\nabla\cdot\Pi_{\mathrm{int}}\sim-\Pi_{\mathrm{int}}/a$ using

[TABLE]

The quantity $\tau^{\mathrm{turb}}$ is shown in figure 13.

Acknowledgements

The authors would like to thank I. Calvo, F. Parra, J. Hillesheim, J. Lee, G. Papp, S. Satake, and J. Harris for helpful input and discussions. This work was supported by the US Department of Energy through grants DE-FG02-93ER-54197 and DE-FC02-08ER-54964. The computations presented in this paper have used resources at the National Energy Research Scientific Computing Center (NERSC).