Physical insights from the aspect ratio dependence of turbulence in negative triangularity plasmas

The discovery of the High confinement mode (H-mode) in ASDEX [1] opened the gates to a new era of research in the field of magnetically confined fusion plasmas. After exceeding a certain threshold of heating power, a tokamak plasma spontaneously creates a transport barrier at the edge, thereby enabling the formation of a region of steep pressure gradient. This is the so-called pedestal which greatly increases the temperature and density in the core of the plasma. The performance of H-mode plasmas is substantially better than standard Low confinement mode (L-mode) plasmas, which do not display a noticeable pedestal. However, the promise of H-mode is tempered by the appearance of powerful instabilities called edge localized modes (ELMs) [2], which are quasi-periodic relaxations of the pedestal that release large amounts of particles and energy into the plasma-facing components of the tokamak. In addition, the power released by an ELM is directly proportional to the power injected in the plasma, making ELMs extremely dangerous for the integrity of a power plant. A single ELM could seriously damage the first wall of the machine and compromise the steady production of energy.

For these reasons, the European roadmap towards DEMO (DEMOnstration power plant) identified the ELMs as one of the most important issues to be solved [3] and recommends searching for methods to reduce the size of ELMs or for scenarios where they can be avoided without compromising confinement. A solution may come from plasma shaping.

Decades of plasma physics research have proven that giving a tokamak plasma a non-circular poloidal cross-sectional shape can greatly modify the stability of the plasma and its confinement of particles and energy. The poloidal cross-sectional shape of a plasma can be modified by changing the two most important shaping factors: elongation

$$\begin{align} \kappa = \frac{Z_\textrm{max}-Z_\textrm{min}}{2a},\end{align} \tag{ 1 }$$

and triangularity

$$\begin{align} \delta = \frac{2R_\textrm{geo}-\left(R_\textrm{top}+R_\textrm{bot}\right)}{2a},\end{align} \tag{ 2 }$$

where $Z_\textrm{max}$ and $Z_\textrm{min}$ are the maximum and minimum vertical positions of the flux surface, $R_\textrm{top}$ and $R_\textrm{bot}$ are the radial locations of the uppermost and lowermost points, $R_\textrm{geo}$ is the geometric center of the flux surface and a is the minor radius.

The configuration that has become the standard scenario in the majority of modern tokamaks is the elongated dee-shaped one, also called positive triangularity (PT) because the value of δ is larger than zero. The term dee comes from the resemblance of the cross-section to the letter 'D', when shown to the right of the tokamak axis of symmetry (as is typical). PT plasmas became popular because they have good MHD stability and can enter H-mode with less heating power. Of course, the downside is the aforementioned ELMs, the seriousness of which has been widely appreciated only more recently. An alternative that has gained interest is negative triangularity (NT), where the plasma is shaped such that δ is less than zero. A large number of experiments on several medium-size tokamaks (on TCV for several decades [4–6], followed more recently by DIII-D [7, 8] and AUG [9]) as well as many numerical studies [10–12], find that NT reduces turbulent transport, and can inhibit the transition to H-mode. These features allow a NT plasma to have an L-mode-like edge pressure profile and avoid ELMs, but with sufficiently increased pressure gradient at the edge and in the outer core to recover H-mode-like values of density and temperature in the core.

Despite the great number of studies on the topic, the physical processes underlying the beneficial effect of NT on turbulent transport are not fully understood yet. Moreover, as observed by Merlo and Jenko in [11], it is key to study how other plasma parameters influence the impact of NT on confinement so we can use this knowledge to optimize. In this work, we will study the dependence of transport on the aspect ratio $A = R_0/r$ (the ratio of major radius R₀ to minor radius r) in NT plasmas. The motivation for such a study is twofold. First, studying how a NT plasma behaves at low A gives insight into the feasibility of a NT Spherical Tokamak (ST), i.e. a tokamak characterized by a small aspect ratio. Second, studying the large aspect ratio limit greatly simplifies the physics and therefore the analysis. Moreover, as we will see, the conclusions will turn out to be valid also at the aspect ratios typical for a standard tokamak. As the aspect ratio is not a parameter that can be easily changed and scanned in real experiments, we make use of high-fidelity flux tube gyrokinetic GENE [13] simulations, which allow us to explore the parameter space and capture the physical processes controlling turbulence.

In this paper, we will see that a non-trivial picture will emerge from the aspect ratio scan, showing that NT has better confinement than PT at large and standard aspect ratio values, regardless of the nature of the turbulent regime. However, at tight A NT becomes detrimental in trapped electron mode (TEM) dominated regimes, while staying beneficial in ion temperature gradient (ITG) mode dominated regimes. This complex dependency on the aspect ratio triggered a deeper study to understand the physical processes underlying the stabilization/destabilization of certain types of turbulence by a NT geometry. While an explanation of the behaviour at small A is still missing, we will present a simple but compelling physical picture of why NT weakens ITG instability at large and standard aspect ratios. We refer to this understanding as simple because it relies only on two quantities, i.e. magnetic drifts and Finite Larmor Radius (FLR) effects, which in turn are easy to calculate from the shape of the plasma. Motivated by similarities between our model and the work done in [12] for TEM turbulence, we also tried to apply our model to a scenario with pure TEM drive, finding that FLR effects and magnetic drifts can only explain part of the physics for this kind of scenario.

Finally, to achieve a better understanding of the small A regime and also use more realistic equilibria, we will consider the ST SMART [14–17]. This is a ST that only very recently started operation at the University of Seville, characterized by strong geometric flexibility, which should enable the creation of single-null and double-null NT and PT configurations. The input parameters for GENE simulations are taken from predictive simulations for SMART made with TRANSP [18]. We will perform only local linear simulations. They will reveal which kind of turbulence dominates and how it responds to a change in various parameters. We will see that, because of the high β typical of STs, large electron temperature gradient and large magnetic shear, turbulence is dominated by micro-tearing modes (MTM). Moreover, we will find that NT is more sensitive to MTMs, with growth rates larger than those computed for PT. However, we will see that it should be possible to operate SMART in an ITG-dominated regime, where it is possible to exploit the beneficial effect of NT (consistently with previous simulations). We were not able to perform reliable NL simulations because of the challenges commonly encountered when simulating scenarios mainly driven by MTMs (e.g. lack of converged heat fluxes when using standard resolutions).

The remainder of this paper is organized as follows. In section 2, we describe briefly how GENE works and how it has been used to carry out the present study. In section 3, we present how the aspect ratio and turbulent regime can modify the beneficial effect of NT on confinement. In section 4, we present the first linear simulations performed for the SMART tokamak, providing predictions for the turbulent regime and the performance of the NT option with respect to PT. In section 5, we motivate and propose a simple physical picture of why NT weakens ITG turbulence in large and standard aspect ratio tokamaks. In section 6, we apply the model described in the previous section to a scenario with pure TEM drive. Finally, section 7 draws some conclusions and provides a final discussion.

This work consists exclusively of numerical simulations performed with the flux tube gradient-driven version of the GENE code [13]. GENE is a physically comprehensive Eulerian gyrokinetic code that solves the Vlasov–Maxwell equations discretized on a 5-dimensional (5D) grid. This leads to a large set of coupled ordinary differential equations that can be solved as an initial value problem or as an eigenvalue problem. The code can be run in a linear mode if the nonlinear terms of the GK equation are neglected, or in a nonlinear mode when the nonlinear terms are retained. The former is considerably faster and is especially useful to quickly gain an understanding of the nature of turbulence. The latter is computationally more expensive but needed to physically model micro-turbulence. The real space is parametrized by a 3D system of non-orthogonal field-aligned coordinates $(x,y,z)$, which correspond respectively to the radial, binormal, and parallel (to the magnetic field B) coordinates. They are defined by

$$\begin{align} \begin{cases} x = x\left(\psi\right) \\ y = C_y\left(q\left(\psi\right)\theta-\varphi\right), \\ z = \theta \end{cases}\end{align} \tag{ 3 }$$

where ψ is the poloidal magnetic flux, ϕ the toroidal angle, θ the straight field-line poloidal angle, q is the safety factor and C_y is a constant. The remaining two coordinates of the 5D discretized grid, $v_\parallel$ and µ, describe velocity space and correspond to the parallel velocity and the magnetic moment (defined as $\mu = (m_\sigma v_\perp^2)/(2B)$, where m_σ is the mass of the spiecies σ and $v_\perp$ the perpendicular velocity). In a flux tube, the background values of density, temperature, and flow (as well as their gradients) are evaluated at the flux surface of interest and kept constant across the simulation domain. This choice is justified if the radial scale of turbulence is much smaller than the radial scale of variation of equilibrium quantities (machine scale), which is fairly well satisfied for TCV and better satisfied in larger devices like DIII-D or reactor-scale devices. The real space simulation domain corresponds to a small rectangle extended in x and y, that follows the magnetic field lines along z. Periodic boundary conditions are used for the perpendicular coordinates x and y, while pseudo-periodic boundary conditions are applied to z. For more details, the interested reader is referred to [19].

To carry out a comprehensive analysis, we will consider scenarios characterized by different turbulent regimes, i.e. regimes with different dominant turbulent drive. Moreover, to maintain realism, we will use discharges produced in TCV and DIII-D, and scenarios predicted for SMART. In total, we will consider six different cases: a realistic TEM dominated scenario (rTEM), a realistic ITG dominated scenario (rITG), a pure density gradient-driven TEM scenario (pTEM), two pure ITG-driven scenarios (pITG-1 and pITG-2), and an MTM dominated scenario (the one relative to SMART).

GENE can use different methods or codes to reconstruct the magnetic equilibrium at a specific flux surface. These can be divided into numerical reconstruction codes (CHEASE [20], gist, and tracer-efit) and analytical models ($s-\alpha$ and Miller). For all the simulations presented in this work, we used the local equilibrium Miller model, because, as an analytical model, it allows us to easily change the geometry by changing a few parameters. As shown in [21], a Miller equilibrium is completely defined by a total of 14 scalar quantities. As shown in equation (4), the parametrization of the flux surface shape requires 6 parameters: the geometric axis R₀, the elevation Z₀, the aspect ratio $A = R_0/r$, the elongation κ, the triangularity δ and the squareness ζ.

$$\begin{align} \begin{cases} R\left(\theta\right) = R_0\left[1+A^{-1}\cos{\left(\theta+\arcsin{\left(\delta\sin{\theta}\right)}\right)}\right] \\ Z\left(\theta\right) = Z_0+\kappa r\sin{\left(\theta+\zeta\sin{\left(2\theta\right)}\right)}, \end{cases}\end{align} \tag{ 4 }$$

where R and Z are the radial and axial cylindrical coordinates as a function of the poloidal angle θ (θ = 0 corresponds to the outboard midplane, while $\theta = \pm\pi$ the inboard midplane). The computation of the poloidal field necessitates 8 additional parameters: the elongation shear s_κ , the triangularity shear s_δ , the squareness shear s_ζ , the safety factor q₀, the magnetic shear $\hat{s}$, the pressure gradient p', and the Shafranov Shift given by $\partial_rR_0$, and the variation of elevation $\partial_rZ_0$.

To perform our simulations, we used Miller to model relevant flux surfaces in actual experimental magnetic equilibria. We flipped the triangularity (from NT to PT or vice versa) and changed the aspect ratio varying the minor radius. During these operations, the values of density and temperature, and their logarithmic gradients, were kept fixed. This procedure is commonly used in GK simulations and enables one to isolate the effect of geometry on transport. The numerical convergence of all simulations has been tested by resolution studies. More details will be given in the next sections, but, depending on the simulation, we used $L_x\sim[110-180]\rho_i$ and $L_y\sim[125-250]\rho_i$ for the radial and binormal widths of the simulation box with $[128-256]$ k_x and $[32-64]$ k_y modes. Depending on the aspect ratio, we used between 32 to 48 grid points in z. The size and resolution of velocity space have been kept fixed: $L_v = 3v_{\textrm{th},\sigma}$, $L_\mu = 9\mu B_0/T_\sigma$, $n_v = 42$ and $n_\mu = 12$, where $v_{\textrm{th},\sigma} = \sqrt{2T_\sigma/m_\sigma}$ is the thermal velocity, B₀ the toroidal magnetic field at R₀ on the flux surface of interest, m_σ the reference mass and T_σ the reference temperature.

To carry out the most realistic analysis possible, we based our parameters on an electron-cyclotron-resonance (ECR) heated TCV NT discharge (#68 145), a neutral beam injection (NBI) heated TCV NT discharge (#69 682), and a high power NBI-heated DIII-D PT discharge (#176 492). We reconstructed the magnetic equilibria with Miller, choosing a radius of analysis (defined here as the normalized square root of the toroidal flux) $\rho_\textrm{tor} = 0.8$ for TCV and $\rho_\textrm{tor} = 0.75$ for DIII-D. Then we flipped the sign of δ and s_δ and we artificially changed the aspect ratio keeping fixed the major radius and changing the minor radius. At the chosen radial locations, the logarithmic gradients were computed from the experimental data and then were kept fixed while the geometry was changed.

We chose discharges #68 145 and #176 492 to have cases dominated by TEM and ITG turbulence respectively. Originally, #69 682 was chosen to provide an ITG-dominated case, but linear GK simulations revealed that it was TEM-dominated. Additionally, we modified the logarithmic gradients of #68 145 and #69 682 to create idealized cases, namely two pure ITG-driven scenarios (pITG-1 and pITG-2) and a pure density gradient-driven TEM one (pTEM). Some of the most important input parameters used in the simulations are reported in table 1, while the shapes of the three original scenarios with their flipped counterparts are shown in figure 1. For the realistic TCV and DIII-D simulations collisions, electromagnetic effects and carbon as the main impurity were retained. The pITG and pTEM simulations neglected impurities, but collisions and electromagnetic effects were kept. All simulations used a kinetic treatement for the electrons.

Zoom In Zoom Out Reset image size

Table 1. Logarithmic gradients of electron temperature $R/L_{Te}$, ion temperature $R/L_{Ti}$, electron density $R/L_{ne}$, carbon density $R/L_{nC}$, as well as the ion-electron $T_e/T_i$ temperature ratio, electron density n_e , the carbon to electron density $n_C/n_e$, the local electron β, the electron-ion collisional frequency ν_ei and geometric parameters needed to specify the Miller equilibrium.

	rTEM	pTEM	pITG-1	pITG-2	rITG
$R/L_{Te}$	14.54	0	0	0	6.37
$R/L_{Ti}$	8.48	0	8.48	12.89	8.98
$R/L_{ne}$	6.01	7.81	0	0	1.47
$R/L_{nC}$	9.16	—	−	—	1.30
Te [keV]	0.41	0.41	0.41	0.31	1.47
$T_e/T_i$	0.98	0.98	0.98	1.37	1.17
$n_e\,[10^{19}m^{-3}]$	1.95	1.95	1.95	1.11	4.25
$n_C/n_e$	0.013	—	−	—	0.018
$\beta [\%]$	0.16	0.16	0.16	0.07	0.71
νei [s−1]	0.84	0.84	0.84	0.838	0.45
A	5	5	5	5	3.7
q	2.01	2.01	2.01	2.20	3.41
$\hat{s}$	1.6	1.6	1.6	2.02	2.43
κ	1.37	1.37	1.37	1.31	1.47
$|\delta|$	0.16	0.16	0.16	0.25	0.21
ζ	−0.07	−0.07	−0.07	0.24	−0.31
sκ	0.19	0.19	0.19	0.16	0.35
$|s_\delta|$	0.33	0.33	0.33	0.59	0.54
sζ	−0.21	−0.21	−0.21	0.12	−0.14

We started by performing linear simulations for each scenario to establish the turbulence drives stated above. These results are presented in appendix A. We proceeded with nonlinear simulations where we scanned the aspect ratio by varying the minor radius of the flux surface, while keeping the major radius fixed. Figure 2 shows the ratio between the nonlinear heat fluxes from the NT cases ($Q^\textrm{NT}$) and the PT cases ($Q^\textrm{PT}$). When the ratio of fluxes is below one, the NT configuration has better confinement, while the opposite is true when the ratio is above one. The shaded areas correspond to the standard deviation of the averaged values computed from the time traces of the heat fluxes.

Zoom In Zoom Out Reset image size

Figure 2 clearly shows that the influence of NT on transport has a complex and non-trivial dependency on aspect ratio. Regardless of the turbulent regime, at large and conventional values of A, NT always has better confinement than PT. In this range of aspect ratio values, the reduction of the heat flux by NT varies between 1.5 for the TEM-dominated regime to a factor of 2.5 for the pITG scenarios. For STs (i.e. A < 2.5), we find a bifurcation. When turbulence is purely ITG or ITG-dominated, NT still greatly stabilizes turbulence, reducing the heat flux by a factor of ${\sim}2$. Instead, for a pure-TEM or TEM-dominated regime, NT becomes detrimental, increasing the heat flux by a factor of ${\sim}2.5$.

Heat flux spectra are plotted in figure 3 for different regimes using two values of aspect ratio: conventional and spherical. We are not showing the spectra for the large aspect ratio cases because they are similar to conventional A. We show only the most important component, i.e. electron or ion channel, depending on the considered regime. From figures 3(a)–(c), which correspond to ITG regimes, we can see that at conventional A the reduction of heat flux is equally distributed among all the considered k_y modes. Instead, when we move to a smaller A, the reduction is more pronounced for modes with $k_y\rho_s$ bigger than 0.2. If we now look at figure 3(d), which corresponds to the realistic TEM case, we can see that the detrimental effect of NT comes from small k_y modes at tight A. These are the modes that are most strongly destabilized when one decreases the aspect ratio for TEM turbulence.

Zoom In Zoom Out Reset image size

To achieve a deeper insight, it is important to study the stiffness of transport, i.e. how the heat fluxes change when the value of the driving gradient is changed. Since the driving gradient is easiest to identify for pITG and pTEM, we have performed simulations for these cases with increased or decreased $R/L_{Ti}$ and $R/L_{ne}$, respectively. The study of the stiffness is extremely important to make PT/NT comparisons at fixed heat flux (i.e. gradient driven simulations).

Figure 4 displays the results. At conventional aspect ratio, flipping the triangularity does not change the stiffness but rather increases the critical gradient, i.e. the gradient at which modes become unstable. This results in reduced turbulent transport and better confinement properties. It also means that, at fixed heat flux, one can produce a plasma with increased gradients. The observation that, in a standard A tokamak, NT primarily affects the critical gradient and not the stiffness is in agreement with many other works [11, 22]. In contrast, at small A we can see the opposite. Indeed, flipping the triangularity does not alter the critical gradient much but does change the stiffness. In the pITG regime NT reduces the stiffness, while in the pTEM regime NT increases it. All these observations indicate that turbulence responds to a change in shape very differently in conventional tokamaks and STs.

Zoom In Zoom Out Reset image size

In this section, we present the first gyrokinetic simulations of the PT and NT scenarios foreseen for SMART.

SMART is a novel, compact, spherical tokamak ($1.4\leq A\unicode{x2A7D}3.0$) that has recently begun operations at the University of Seville. SMART has a strong geometric flexibility that allows it to produce plasmas with an elongation of $\kappa\unicode{x2A7D}3.0$ and a triangularity between $-0.6\unicode{x2A7D}\delta\unicode{x2A7D}0.6$. Both single-null and double-null configurations can be produced. For more technical details, the interested reader is referred to [14–17].

Given the results of the previous section, experimental results from SMART will be very interesting. All the previous simulations are case studies and can partially predict how a NT ST will behave in the limit of certain idealized regimes. To more realistically model SMART, the input parameters for GK simulations have been obtained from TRANSP predictive simulations for the PT scenario [18]. In figure 5 we show the magnetic equilibria and in figure 6 we show the kinetic profiles predicted by TRANSP for the PT scenario. This is an ohmic heated scenario and, since intrinsic rotation is hard to predict, rotation and flow shear were not considered in the following gyrokinetic simulations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

All linear simulations have been performed using Miller parametrization to construct a local magnetic equilibrium at $\rho_\textrm{tor} = 0.88$. We chose this location because it is sufficiently close to the edge to have large $|\delta|$, but sufficiently far from the edge to avoid the region with steep gradients. The simulations retained kinetic electrons, collisions and electromagnetic effects. Other important parameters used in the simulations are listed in table 2. As before, to isolate the effect of the geometry, we kept these parameters fixed and flipped the shape. In this section, we will only show linear simulations as we were not able to achieve reliable NL simulations, which are subject to a future publication.

Table 2. GENE input parameters at $\rho_{tor} = 0.88$.

$R/L_{Te}$	15.1
$R/L_{Ti}$	9.2
$R/L_{n}$	2.0
Te [keV]	0.10
$T_i/T_e$	0.99
$n_{e}\;[10^{19}$m$^{-3}]$	2.83
$\beta_e [\%]$	0.56
νei	8.62
A	2.2
q	3.33
$\hat{s}$	3.71
κ	1.76
$|\delta|$	0.16
ζ	−0.04
sκ	0.66
$|s_\delta|$	0.96
sζ	−0.28

As already mentioned, one has to be careful in assessing the performance of a scenario using only linear simulations. Indeed, the physics of nonlinear saturation is missing and we cannot compute the final heat fluxes without it. However, a linear analysis can allow one to determine the dominant types of instability as well as to study the linear stability of the modes. We started with a k_y scan (at $k_x = 0$) to determine the type of turbulence that dominates at the ion scale, i.e. $k_y\rho_s = [0.05-2.0]$. To be consistent with the work shown in previous sections, we did not consider electron scale turbulence (i.e. $k_y\rho_s\gg2.0$). However, Electron Temperature Gradient (ETG) modes can play an important role in spherical tokamaks [23, 24]. This is left as future work. The results of the simulations at the ion scale are displayed in figure 7. Comparing the magnitude of the growth rates, we see that NT has larger values than PT. This does not necessarily imply larger nonlinear heat fluxes, but does mean that the modes grow faster. However, the most interesting observation is the nature of turbulence. For the entire ion scale, the spectrum is dominated by MTMs. Indeed, looking at the electrostatic potential φ and the parallel vector potential $A_\parallel$ in figure 8, we can see that these modes have tearing parity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is well-known that STs can be dominated by MTM turbulence [25, 26]. MTM turbulence was not significant in the previous sections, as β and the ratio of electron to ion temperature gradient were smaller. Therefore, to better characterize MTM turbulence in STs and how it behaves depending on triangularity, we performed parameter scans for the $k_y\rho_i = 0.4$ mode.

Figure 9 shows the results of a β scan. First, we can identify two distinct regions depending on the value of β. For $\beta \lt 0.2\%$, both PT and NT are dominated by ITG turbulence. It is interesting to notice that we recover the same result that we observed for the TCV and DIII-D low A cases: NT has ITG modes with smaller growth rates. At larger values of β we can identify a transition to the MTM regime. Importantly, NT's transition to MTM-dominated regime occurs at $\beta\sim0.2\%$, while PT's transition is at $\beta\sim0.3\%$. This suggests that NT is more sensitive to electromagnetic instabilities. At even larger β, we can see that MTMs are always stronger in NT than PT. The reason for this behaviour is still unclear and is a topic for future work.

Zoom In Zoom Out Reset image size

Next, we carried out a scan in collisionality for the $k_y\rho_i = 0.4$ mode, which is shown in figure 10. Again in this parameter scan we see a mode transition. At exactly zero collisionality, as predicted by previous theoretical work [27], we see that MTM cannot exist and instead the most unstable mode is TEM. For this TEM instability, NT has a lower growth rate than PT. However, as soon as the collisionality is finite, MTM dominates and NT has a faster growth rate. If we then increase collisionality, it reduces the gap between the growth rates of the two geometries. In figure 11, we show a scan of the electron temperature gradient, which is another important drive of the MTM instability. Once again, we can see that PT and NT transition to MTM at somewhat different values, and that once the MTMs become dominant, the modes are more unstable in NT.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Another parameter that influences MTMs is the magnetic shear $\hat{s}$. Figure 12 shows the growth rate and frequency of the $k_y\rho_i = 0.4$ mode as the magnetic shear is changed. We see that NT and PT have very different thresholds for the onset of the MTM regime. Indeed, NT transitions at a lower value of $\hat{s}$ for reasons not currently understood.

Zoom In Zoom Out Reset image size

This linear GK analysis allowed us to identify the dominant type of turbulence in the PT and NT equilibria expected for SMART. In these equilibria both shapes should be dominated by MTM at the ion scale, and NT should have a stronger linear drive. However, we cannot make strong statements about the real performance of the NT option in SMART. Our linear simulations do not predict the nonlinear saturation of these modes and they did not consider the effect of impurities. However, based on the results obtained in this work, we believe that NT will be most beneficial for SMART in a regime where MTMs are less important and ITGs dominate. For this reason, we performed a triple scan in electron and ion temperature gradients and magnetic shear. The results are shown in figures 13 and 14.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

First, we observe that an ITG region exists in both PT and NT and increases as the magnetic shear is reduced. In this region, where ITG turbulence dominates, linear growth rates are always lower in NT (the green regions in figure 14). The opposite is true when we have MTMs. Therefore, it should be possible to operate in ITG-dominated conditions, where NT should be beneficial, but the size of this ITG region is considerably larger in PT than NT, making the requirements more strict. To conclude, if we want to operate SMART in a regime where we are confident to have a beneficial effect from NT, we should reduce the magnetic shear, increase the ion temperature gradient and decrease the electron temperature gradient.

As shown in sections 3 and 4, the aspect ratio dependence of transport in NT and PT plasmas is complicated. To better understand it, we started by studying the scenario that we expected to be the simplest: the large A pITG case. Indeed, in the limit of $A\gg$1 ITG modes depend only on two main geometry-related quantities: FLR effects and magnetic drifts. We then applied this picture to the conventional and tight A cases, where many other factors arising from geometry can affect turbulence.

To make the derivation as clear as possible, we give here an outline of the contents of this section. First, we will illustrate how magnetic drifts and FLR effects influence the strength of ITG modes. Then, we will show how the geometry of the cross-sectional shape enters the FLR effects and magnetic drifts. This will allow us to show how NT weakens ITG. Lastly, we will show the numerical results that support this simple theory.

5.1. Impact of magnetic drifts and FLR effects on ITG

In this subsection, we will employ [28] and [29] to show how magnetic drifts and FLR effects contribute to the strength of ITG turbulence.

Starting from the standard physical picture of the toroidal ITG instability [28], one can derive the relation

$$\begin{align} \omega_\textrm{D} = \frac{\left(\tau+b\right)\eta}{1+3\left(\tau+b\right)^2}\omega_*,\end{align} \tag{ 5 }$$

that determines the condition for which the growth rate γ of an ITG mode will be maximized. Here ω_D is the magnetic drift, $\tau = T_i/T_e$ is the ion to electron temperature ratio, $b = k_\perp^2\rho_i^2$, ρ_i is the ion Larmor radius, $k_\perp$ is the perpendicular wavenumber, $\eta = L_{ni}/L_{Ti}$, where $L_{ni} = -\textrm{d}(\ln{n_i})/\textrm{d}x$ is the ion density logarithmic gradient, $L_{Ti} = -\textrm{d}(\ln{T_i})/\textrm{d}x$ is the ion temperature logarithmic gradient, $\omega_*$ is the diamagnetic drift frequency. All the details regarding this derivation, including the definition of all the mentioned quantities, are reported in appendix B.

If we set τ = 1 and b = 0 and we recast the formula in terms of velocities, the condition becomes

$$\begin{align} \textbf{v}_\textrm{D} = \frac{\eta}{4}\textbf{v}_*.\end{align} \tag{ 6 }$$

This equation shows that to have the maximum possible growth rate, the magnetic drift and the diamagnetic drift velocities have to be of the same order. Since this is a very simplified model, one should not focus on the actual value of the constant that relates the two velocities. Rather, we should focus on the fact that if the magnetic drift velocity of the ions $\textbf{v}_\textrm{D}$ is too high or too low with respect to the diamagnetic velocity $\textbf{v}_*$, the growth rate will be small. In simple terms, if the depletion/accumulation zones of ions are traveling too fast or too slow with respect to the ion drift wave, a mismatch between the two develops, making the mode less unstable. The need to satisfy this 'resonance' is key in the explanation of NT beneficial effect on ITGs.

Apart from magnetic drifts, the only other way that the flux surface shape enters into the GK model in the $A\gg1$ limit is through the FLR effects. Specifically, in the local GK equations, the radial dependence of equilibrium quantities is dropped by the application of a gyroaverage, i.e. an integral average along the fast gyromotion of particles. It can be shown [30] that the gyroaverage of a scalar field $\mathcal{F}$ reduces to a multiplication of the field with the zeroth order Bessel function J₀, i.e.:

$$\begin{align} <\mathcal{F}\left(\textbf{X}\right)> = \sum_{k_y,k_x}\mathcal{F}\left(k_x,k_y\right)\textrm{e}^{\textrm{i}\textbf{k}_\perp\cdot\textbf{X}}J_0\left(\rho_\sigma k_\perp\right),\end{align} \tag{ 7 }$$

where k_x and k_y are the radial and binormal Fourier modes and ρ_σ is the Larmor radius of the σ species. Therefore, the FLR effects correspond to a damping of the mode. The larger $\textbf{k}_\perp$, the stronger the damping.

5.2. Impact of plasma shape on magnetic drifts and FLR effects

FLR effects enter the GK model through the zeroth order Bessel function $J_0(k_\perp^2\rho_s^2)$. Since we are interested in the geometric dependence of this quantity, it is important to note that $\textbf{k}_\perp$ is the perpendicular wavenumber of a mode and is defined as

$$\begin{align} \textbf{k}_\perp = k_x\boldsymbol{\nabla}x+k_y\boldsymbol{\nabla}y.\end{align} \tag{ 8 }$$

Therefore, using the quantities defined in appendix C, we can write $k_\perp^2$ as

$$\begin{align} k_\perp^2 = g^{xx}k_x^2+2g^{xy}k_xk_y+g^{yy}k_y^2.\end{align} \tag{ 9 }$$

Hence, the geometry of the cross-sectional shape enters equation (9) through g^xx , g^xy and g^yy . On the other hand, the magnetic drift frequency is $\omega_\textrm{D} = \textbf{v}_\textrm{D}\cdot\textbf{k}_\perp$. Using equation (B3), we see it can be written as

$$\begin{align} \omega_\textrm{D}& = \left[\frac{1}{B\Omega}\left(\frac{v_\perp^2}{2}+v_\parallel^2\right)\textbf{b}\wedge\boldsymbol{\nabla}B\right]\cdot\left[ k_x\boldsymbol{\nabla}x+k_y\boldsymbol{\nabla}y\right] \nonumber\\ & = \Lambda\left(k_x\underbrace{\left(-\partial_y B-\frac{\gamma_2}{\gamma_1}\partial_z B\right)}_{v_{Dx}}+k_y\underbrace{\left(\partial_x B -\frac{\gamma_3}{\gamma_1}\partial_z B\right)}_{v_{Dy}}\right), \end{align} \tag{ 10 }$$

where the scalar Λ comprises all constant terms, v_Dx is the radial component of the magnetic drift velocity and v_Dy the binormal one. The derivation of this equation is reported in appendix D. Therefore, the geometry of the plasma shape determines the drift frequency through the derivatives of the magnetic field and the elements of the metric tensor. Note that, with the GENE sign convention, a negative value of ω_D corresponds to the bad curvature region and positive values to the good curvature region.

Now that we have expressions for the FLR effects and magnetic drifts, we can see how triangularity affects them. It must be pointed out that these two quantities, contrary to the simplified model presented in the previous subsection, depend on the poloidal angle, and are thus a function of z. Moreover, both FLR effects and magnetic drifts depend on k_x and k_y . To simplify our analysis, we will consider the $k_x = 0$ mode as it is typically the fastest growing linear mode for standard turbulence. Doing so allows us to test k_y as a prefactor and will use only g^yy as a proxy for the FLR effects and v_Dy as a proxy for the magnetic drifts.

By using Miller local equilibrium geometry specification [31], it is straightforward to calculate g^yy and v_Dy for any value of triangularity. In figure 15(a) we take the equilibrium of the pITG-2 scenario and progressively change the triangularity. Figure 15(a) shows the flux surface shapes, figures 15(b) and (c) show v_Dy and g^yy as functions of the poloidal angle for different values of triangularity.

Zoom In Zoom Out Reset image size

We see that, as the triangularity is decreased, the FLR effects become stronger in the vicinity of the outboard midplane, while the magnetic drifts become faster, i.e. more negative. In light of the observations made in the previous section, it is straightforward to say that one mechanism contributing to the stabilization of ITG in NT geometry comes from stronger FLR effects at the outboard midplane, where ITGs are stronger. However, the impact of magnetic drift is not clear. As mentioned before, negative values of v_Dy correspond to bad curvature. However, more negative values do not necessarily mean stronger instability. Indeed, as we showed in the previous section, ITG is destabilized by matching the magnetic drift velocity and the diamagnetic drift velocity. Therefore, more negative values of v_Dy can imply weaker ITG if v_Dy is already larger than the condition set by equation (6). In the next section we will test this hypothesis with linear and nonlinear simulations.

5.3. Numerical results

We performed linear and nonlinear simulations to test if the faster magnetic drifts and stronger FLR effects of the NT geometry are responsible for weaker ITG turbulence.

5.3.1. Linear results—Large A.

We started with linear GK simulations. The main goal was to verify that the simple gyrofluid model presented in sections 5.1 and 5.2 holds.

We considered the large A pITG-2 scenario and modified GENE to artificially change the poloidal dependence of v_Dy , keeping everything else untouched. We modified only v_Dy because in linear simulations we consider only $k_x = 0$ modes as they are assumed to grow the fastest. We took the original NT v_Dy profile and modified it according to

$$\begin{align} v_{Dy}\left(z\right) = \sigma v_{Dy}^{NT}\left(z\right)+v_{Dy0},\end{align} \tag{ 11 }$$

where σ and $v_{Dy0}$ are arbitrary scalar quantities that allowed us to change the variation of the profiles with poloidal angle and the offset respectively. We considered four values of σ (i.e. 0.08, 0.3, 1.0 and 1.5) and for each of them we changed the offset to shift the profile up and down. The considered v_Dy profiles are shown in figure 16. For each modified profile, we performed a linear simulation for the $(k_x,k_y) = (0,0.4)$ mode, which is the fastest growing one and usually one of the modes that contributes the most to the NL heat fluxes.

Zoom In Zoom Out Reset image size

This procedure allows us to isolate the effect of the magnetic drift and disentangle the contribution of the modulation from the value at the outboard midplane (where turbulence is strongest). Figure 17 shows the results of the simulations. The different curves correspond to different modulations σ, each of which is plotted against the value of v_Dy at z = 0 in order to take into account the shift of the profile. We can see that the linear simulations confirm and deepen the simple picture that we gave at the beginning of the section. First, regardless of the modulation of the drift profile, a maximum point is always present at $v_{Dy}\sim-1$. Given the input parameters of our simulations $v_*\sim-13$ (in normalized units), so this optimal point is located at $v_{Dy}{\sim}v_*/13$ (as highlighted by the light blue box in figure 17). On the other hand, we see that the variation of v_Dy with poloidal angle tends to stabilize the mode.

Zoom In Zoom Out Reset image size

With this information we can better understand how geometry influences ITG via magnetic drifts. In figure 18, we plot the ideal $v_{Dy}(z)$ to enable the resonance together with the actual $v_{Dy}(z)$ in the NT and PT cases. We see that v_Dy varies less with poloidal angle in a PT geometry than in a NT one, and it is always closer to the resonance. Thus, the proximity to the resonance destabilizes ITG modes in PT geometry, explaining why NT is beneficial. However, linear simulations cannot be taken as a final proof. We also need nonlinear simulations to test our understanding in a realistic quasi-stationary turbulent state.

Zoom In Zoom Out Reset image size

5.3.2. Nonlinear results—Large A.

In this subsection, we show a series of nonlinear simulations to confirm that, at large A, stronger FLR effects and faster magnetic drifts fully explain the beneficial effect of NT on pure-ITG turbulence.

Nonlinear GK simulations are much more computationally expensive than linear simulations, which prevents an analogous scan to figure 17. Instead, we artificially swap all the geometric coefficients responsible for FLR effects and magnetic drifts between NT and PT. Indeed, from equations (9) and (10), we can easily identify all the coefficients that one needs to exchange between the two geometries. Clearly, this procedure creates unrealistic equilibria, but it is perfect to isolate the contributions coming from FLR effects and magnetic drifts. Apart from these two effects, everything else is left untouched in the GK equations. A similar exercise was carried out by Merlo and Jenko in [11], but only for linear cases.

We will consider the large A pITG-2 scenario and perform new simulations where FLR effects or magnetic drifts $v_\textrm{D}$ (or both) are changed. All simulations were performed with a kinetic treatment of electrons. Figure 19(a) displays the results. The fourth and fifth color bars of figure 19(a) correspond to the self-consistent reference NT and PT cases. From the third color bar, we can see that, by imposing PT's v_D in the NT geometry, we can increase the heat flux by a factor of 1.3. From the second column, we can see that by changing only the FLR effects we can increase the heat flux by a factor of 1.5. This is a very important result, because it confirms our simple physical picture showing that that FLR and magnetic drifts are destabilizing in PT and they have similar importance. If we then look at the first column, i.e. the one where both coefficients were changed, we see that the new heat fluxes closely match the ones of the self-consistent PT case. As expected, this confirms that in the large aspect ratio limit turbulence can be described just in terms of FLR effects and magnetic drifts. The same arguments are supported by columns six and seven. We also performed the same exercise for the large A pITG-1 scenario, which is still characterized by pITG turbulence, but has different input parameters and magnetic equilibria. The results are shown in figure 19(b) and they show the very same findings observed before. These observations confirm that, in the large A limit, the beneficial effect of NT on ITG comes from a stabilizing effect arising from stronger FLR effects and faster magnetic drifts with respect to the resonance condition set by the diamagnetic drift velocity and thus by the ITG.

Zoom In Zoom Out Reset image size

Encouraged by these results, we tried to apply our physical picture to the conventional and tight A pITG cases and to other large A pITG cases where other geometric parameters were changed apart from triangularity. In the following subsection, we show the results for the conventional and tight A pITG cases, where we applied the very same methodology showed before. In appendix E we show the results for scans in other geometric parameters (i.e. elongation and magnetic shear).

5.3.3. Nonlinear results—Conventional and small A.

We performed the same numerical study from the previous section, but with conventional A. The results of nonlinear simulations for pITG-2 and pITG-1 cases are shown in figure 20. We see that stronger FLR effects and faster magnetic drifts are also stabilizing at this aspect ratio. However, at finite aspect ratio the geometry enters via many more effects (e.g. particle trapping, parallel streaming), so it is not necessarily true that these two effects alone can explain the beneficial effect of NT. We see that for both scenarios, FLR effects and magnetic drifts have a dominant effect on the heat fluxes. For the pITG-1 case, they can entirely explain the stabilizing effect of NT. For the pITG-2 case, they account for approximately two-thirds of the stabilization. Therefore, also in the conventional aspect ratio cases, stronger FLR effects and faster magnetic drifts explain the beneficial effect of NT.

Zoom In Zoom Out Reset image size

As a final remark, we also investigated to what extent our physical picture can explain the effect of NT at tight A. We can see from figure 21 that magnetic drifts and FLR effects are no longer sufficient to explain the stabilization of ITG turbulence by NT. Contrary to what was shown before, we see complicated and counter-intuitive results. When we impose the FLR effects and magnetic drifts from NT into the PT geometry, turbulence is reduced as in previous cases with A = 50 and A = 5. On the other hand, when we impose the FLR effects and magnetic drifts from PT into the NT geometry, we still see a stabilization. These trends are clearly in contrast with the physical picture that we held for large and conventional A, suggesting that at tight aspect ratio other physical processes are involved. The investigation of these is left as future work.

Zoom In Zoom Out Reset image size

TEMs are driven unstable by a resonance in phase space between the diamagnetic drift velocity of the underlying drift wave and the toroidal precessional drift of trapped electrons. Based on this fundamental drive, one can construct a similar physical picture as we did for ITGs. This has already been done in [12]. It can be shown that the toroidal precessional drift is closely related to the bounce average of the magnetic drift velocity of trapped particles [32].

In this section, we will evaluate the NT TEM physical picture considered in the aspect ratio scan. We had to consider a conventional aspect ratio scenario to have sufficiently strong TEM instability. At large aspect ratio, the fraction of trapped particles is very low, so the strength of TEM turbulence is very weak. This complicates the physical picture because, as already mentioned, only at large aspect ratio we can reduce the effect of geometry just to FLR effects and magnetic drifts. However, even if we cannot explain everything in terms of these two factors, it is important to understand the significance of their role in the stabilization of TEMs by NT.

We proceeded with a similar methodology as used for the NL simulations of the pITG cases: we swapped the magnetic drifts, FLR effects and other effects between NT and PT.

Figure 22 shows the results. The fifth and sixth color bars show the NT and PT self-consistent cases, i.e. all the physical factors are computed self-consistently with the respective geometry. The fourth column shows the case where only magnetic drifts have been swapped between NT and PT. We can see that the heat fluxes rise, but the increase is mild. If we move further to the left, the next bar corresponds to the case where the FLR effects have been swapped along with magnetic drifts. The heat fluxes are increased further, but the total heat flux is still much lower than in the PT case. The same trend can be observed if we do the same exercise starting from the PT scenario, as shown by the seventh bar. Therefore, we can conclude that, in TEM-dominated regimes, faster magnetic drifts and stronger FLR effects are only some of the factors contributing to the reduction of turbulence by NT. To find the missing pieces of the puzzle, we swapped other geometric factors that could play a role in the stabilization of TEM turbulence. The results correspond to the first and second bar of figure 22. First, we swapped the geometric parameters responsible for the parallel dynamics, i.e. the parallel advection term and the trapping term [33]. The parallel dynamics can play an important role in TEM turbulence because they are related to the bouncing points of the trapped particles. Swapping the parallel dynamics (in addition to the magnetic drifts and FLR effects) is showed by the second column. We observe that this increases the total heat flux, but it is still lower than in the PT case. Therefore, the parallel dynamics are part of the story, but not the most important factor for the stabilization of TEM by NT. We also tried to swap the geometric parameters that enter the computation of the collisional operator (when expressed in parallel velocity $v_\parallel$ and magnetic moment µ coordinates), namely the magnetic field strength and the Jacobian. This effect is important because it regulates the diffusion of particles across the trapping-passing boundary, influencing the fraction of trapped particles and thus TEM turbulence. Accordingly, if we look at the first column, where we have also swapped the collisional diffusion between NT and PT, we can see that it is possible to recover the same heat flux of the self-consistent PT case. Moreover, the change in the heat flux is larger than that caused by all the previous geometrical effects.

Zoom In Zoom Out Reset image size

For the pTEM case, we can conclude that the physics underlying the confinement benefits of NT is not due to a single physical effect but instead a synergistic combination. This exercise is similar to previous work [12], but the results are slightly different, potentially due to the different equilibria considered. Indeed, [12] was able to recover the heat flux of the opposite triangularity scenario just by swapping magnetic drifts and FLR effects. Here, they explain only part of the picture.

To conclude, in this section we showed that FLR effects and the resonance between magnetic drifts and diamagnetic drifts can be applied, to some extent, to explain how NT stabilizes TEM turbulence. However, they alone are not sufficient to explain the entire effect. Parallel dynamics and collisional diffusivity also play a crucial role.

This work consisted of a thorough numerical investigation of the interplay between aspect ratio and the beneficial effect of NT on turbulent transport. It consisted of three distinct, but complementary studies.

• In the first part, we considered five scenarios characterized by different types of turbulence: a realistic TCV-inspired TEM-dominated scenario, a realistic DIII-D-inspired ITG-dominated scenario, two idealized pure-ITG scenarios, and an idealized pure-TEM scenario. We observed that at large and conventional aspect ratio, NT had better confinement than PT regardless of the turbulent regime. At tight aspect ratios, for ITG-dominated turbulence, NT remained beneficial with respect to PT, while NT was detrimental in TEM-dominated scenarios, leading to heat fluxes larger by a factor of up to 2.5. An analysis of the stiffness showed that at conventional aspect ratio, flipping the sign of triangularity from PT to NT increases the critical gradient and leaves the stiffness unaffected. For small A, flipping the triangularity does not influence the critical gradient, but can change the stiffness. Specifically, for TEM NT is stiffer, while for ITG scenarios PT is stiffer.
• In the second part, we performed the first GK simulations of the NT and PT equilibria predicted for SMART. Linear analysis showed that MTMs dominate at ion scales. Moreover, in NT they have larger linear growth rate when compared to PT with the same kinetic parameters. A scan in electron β showed that MTMs appear at lower β and $R/L_{Te}$ in NT and have larger growth rates. However, when ITG is the dominant type of instability, NT has weaker linear growth rates, which is compatible with the findings of the previous section. Finally, a scan of the magnetic shear showed that large values are crucial to drive MTMs at large β and electron temperature gradients. A triple scan in ion temperature gradient, electron temperature gradient and magnetic shear showed that by lowering the magnetic shear and the ratio of electron temperature gradient to ion temperature gradient it is possible to extend the space where SMART NT option can be operated in ITG-dominated regime and exploit the beneficial effect of NT.
• In the third part, we started from a simplified model of the toroidal ITG instability and were able to explain the beneficial effect of NT on confinement for ITG at large and conventional A. Using linear and nonlinear simulations, we showed that the stabilization of ITG in NT comes from faster magnetic drifts with respect to a resonance condition set by the diamagnetic drift, and from stronger FLR effects. However, even though the beneficial effect of NT is still present at small A in ITG regimes, the physical picture developed in this paper does not apply as many other physical effects are important.

We conclude that the interplay between aspect ratio and NT is complicated, but must be taken into account for the successful optimization of a possible NT fusion reactor. Based on our observations, the strong beneficial effect of NT can be exploited at any aspect ratio if we operate in an ITG-dominated regime. However, based on the result for SMART, we conclude that in a ST, MTMs may well be the dominant instability. If so, NT will likely have more unstable turbulence than PT. However, NL simulations are needed to confirm this picture. Future work will focus on understanding if the observations made for MTM turbulence in spherical tokamaks hold in standard tokamaks as well.

We acknowledge Jason Parisi for the extremely insightful conversation on curvature drive for ITG and ETG turbulence, which paved the way for the explanation of why NT stabilizes ITGs. We acknowledge Stephan Brunner for all the invaluable suggestions and for pointing us to M. Beer's Ph.D thesis. We also thank Giovanni Di Giannatale, Alessandro Geraldini, Antoine Hoffmann, Fabien Jeanquartier, Oleg Krutkin, Haomin Sun and Arnas Volcokas for the constant help and inspiring discussions. This work has been carried out within the framework of the EUROfusion Consortium, via the Euratom Research and Training Programme (Grant Agreement No 101052200 - EUROfusion) and funded by the Swiss State Secretariat for Education, Research and Innovation (SERI). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union, the European Commission, or SERI. Neither the European Union nor the European Commission nor SERI can be held responsible for them. This work was also carried out within the Theory, Simulation, Verification and Vailidation (TSVV) task 2. The authors acknowledge the CINECA award under the ISCRA initiative and the EUROfusion WP-AC, for the availability of high-performance computing resources and support. This work was also carried out (partially) using supercomputer resources provided under the EU-JA Broader Approach collaboration in the Computational Simulation Centre of International Fusion Energy Research Centre (IFERC-CSC). This work was supported in part by the Swiss National Science Foundation.

All data that support the findings of this study are included within the article (and any supplementary files).

In this appendix, we show the results of the linear simulations performed for the different scenarios shown in figure 2 to determine the turbulent regime as the aspect ratio is changed. Figure 23 shows the growth rates and the frequencies as functions of the binormal wavenumber for the $k_x = 0$ mode. We see that, as expected, the real TEM and pTEM cases are dominated by TEM turbulence while the pITG-1, pITG-2 and real ITG are dominated by ITG turbulence.

Zoom In Zoom Out Reset image size

In this appendix we report the derivation that leads to equation (5). Toroidal ITG modes are driven by the magnetic field line curvature and gradients in the magnetic field strength. They can be unstable where the curvature is 'bad', meaning the magnetic field curvature (or equivalently gradient) has a component parallel to the temperature gradient.

Let's consider a simplified slab domain, where the magnetic field points along z and the temperature gradient points in the −x direction. This would correspond to the outboard midplane in a tokamak and is shown in figure 24(a). A small perturbation of the ion pressure will generate a neutrally stable ion drift wave that moves downwards with a velocity $\textbf{v}_*=\frac{\textbf{B}\wedge{\boldsymbol{\nabla}} P}{qNB^2}{qNB^2}$, also called diamagnetic drift velocity, where P is the pressure, q the charge of the ion, N the density and B is the magnetic field. Now, if we include the presence of the magnetic field curvature κ and gradient $\boldsymbol{\nabla}B$, the ions will drift with a velocity $\textbf{v}_\textrm{D}$, which is called the magnetic drift velocity. It is the sum of two terms, one from the curvature v_c and one from the gradient $\textbf{v}_{\nabla B}$, and is expressed as

$$\begin{align} \textbf{v}_\textrm{D}& = \textbf{v}_{c}+\textbf{v}_{\nabla B} \nonumber\\ & = \frac{v_\parallel^2}{\Omega}(\boldsymbol{\nabla}\wedge \textbf{b})_\perp+\frac{v_\perp^2}{2B\Omega}\textbf{b}\wedge\boldsymbol{\nabla}B \nonumber\\ & = \frac{v_\parallel^2}{\Omega B}\textbf{b}\wedge\left(\boldsymbol{\nabla}B+\frac{4\pi}{B}\boldsymbol{\nabla}P\right)+\frac{v_\perp^2}{2B\Omega}\textbf{b}\wedge\boldsymbol{\nabla}B, \end{align} \tag{ B1 }$$

where $\textbf{b} = \textbf{B}/B$ is the unit vector parallel to B, Ω is the Larmor frequency, and $v_\parallel$ and $v_\perp$ are the components of the velocity parallel and perpendicular to the magnetic field respectively. If we neglect pressure gradient (i.e. in the $\beta = P/B^2/2\mu_0\ll1$ limit), we can write for simplicity

$$\begin{align} \textbf{v}_\textrm{D}& = \frac{1}{B\Omega}\left(\frac{v_\perp^2}{2}+v_\parallel^2\right)\textbf{b}\wedge\boldsymbol{\nabla}B. \end{align} \tag{ B2 }$$

We can estimate the characteristic value of this velocity as

$$\begin{align} \textbf{v}_\textrm{D} = \frac{v_\textrm{th}^2}{B\Omega}\textbf{b}\wedge\boldsymbol{\nabla}B,\end{align} \tag{ B3 }$$

where $v_\textrm{th} = \sqrt{T/m}$ is the ion thermal speed. Therefore, at the outboard midplane where κ and $\boldsymbol{\nabla}B$ point in the −x direction the ions will start drifting downwards as pictured in figure 24(b).

Zoom In Zoom Out Reset image size

The temperature dependence in the expression for $\textbf{v}_\textrm{D}$ will make the ions in the cold regions drift slower than those in the warmer regions, creating zones of accumulation and depletion of ions. This mechanism is showed in figure 24(c). The subsequent charge separation creates an electric field that leads to another type of drift, the ExB drift, given by

$$\begin{align} v_{E\wedge B} = \frac{\textbf{E}\wedge\textbf{B}}{B^2}.\end{align} \tag{ B4 }$$

The direction of the ExB drift is crucial to determining stability. Indeed, at the outboard midplane, the ExB drift will reinforce the initial perturbation making the mode unstable, as illustrated in figure 24(d). Instead, if we are on the inboard midplane, the direction of the drift will be such to counter the initial perturbation, stabilizing the mode. For this reason, these two regions are called regions of 'bad' and 'good' curvature respectively.

One can derive a dispersion relation for the ITG modes, which corresponds to the qualitative physical picture [28]. In gyrofluid theory, it is given by.

$$\begin{align} &\omega^2\left(\tau+b\right)+\omega\left(\omega_*-2\omega_\textrm{D}-6\omega_\textrm{D}\left(\tau+b\right)\right) \nonumber\\ & \quad +6\omega_\textrm{D}^2\left(1+\tau+b\right)+2\omega_\textrm{D}\omega_*\left(\eta-2\right) = 0, \end{align} \tag{ B5 }$$

where $\tau = T_i/T_e$ is the ion to electron temperature ratio, $b = k_\perp^2\rho_i^2$ (which is assumed to be small), ρ_i is the ion Larmor radius, $k_\perp$ is the perpendicular wavenumber, $\eta = L_{ni}/L_{Ti}$, where $L_{ni} = -\textrm{d}(\ln{n_i})/\textrm{d}x$ is the ion density logarithmic gradient, $L_{Ti} = -\textrm{d}(\ln{T_i})/\textrm{d}x$ is the ion temperature logarithmic gradient, $\omega_\textrm{D} = \textbf{v}_\textrm{D}\cdot\textbf{k}_\perp$ is the magnetic drift frequency and $\omega_* = \textbf{v}_*\cdot\textbf{k}_\perp$ is the diamagnetic drift frequency. Equation (B5) can be solved to find

$$\begin{align} \omega = \frac{2\omega_\textrm{D}+6\omega_\textrm{D}\left(\tau+b\right)-\omega_*}{2\left(\tau+b\right)}\pm\frac{\sqrt{\left(\omega_*-2\omega_\textrm{D}-6\omega_\textrm{D}\left(\tau+b\right)\right)^2-4\left(\tau+b\right)\left(6\omega_\textrm{D}^2\left(1+\tau+b\right)+2\omega_\textrm{D}\omega_*\left(\eta-2\right)\right)}}{2\left(\tau+b\right)}.\end{align} \tag{ B6 }$$

We can focus on the second term, as it corresponds to the growth rate γ of the mode. Moreover, to further simplify the expression, we will consider the flat density limit, which corresponds to $\eta\rightarrow\infty$ and $\omega_*\rightarrow0$. Therefore, we can write

$$\begin{align} \gamma = \pm\frac{\sqrt{2\left(\tau+b\right)\omega_\textrm{D}\omega_*\eta-3\omega_\textrm{D}^2\left(\tau+b\right)^2-\omega_\textrm{D}^2}}{\left(\tau+b\right)}.\end{align} \tag{ B7 }$$

Therefore, the condition for which the growth rate of the mode γ is maximized is given by

$$\begin{align} \omega_\textrm{D} = \frac{\left(\tau+b\right)\eta}{1+3\left(\tau+b\right)^2}\omega_*.\end{align} \tag{ B8 }$$

This is the optimal condition for the onset of an unstable toroidal ITG.

GENE, like many other GK codes, uses a field-aligned coordinate system to discretize the domain. This choice is made to exploit the large anisotropy of turbulence in order to reduce the computational cost. The field-aligned coordinate system is made of three coordinates: x selects a flux surface, y selects a field line on the flux surface and z selects a point on the field line. The magnetic field can be written in Clebsch form as

$$\begin{align} \textbf{B} = \boldsymbol{\nabla}x\wedge\boldsymbol{\nabla}y.\end{align} \tag{ C1 }$$

It is clear that this is not an orthogonal coordinate system, but a curvilinear one. Hence, the standard algebra for curvilinear systems can be applied [34]. We can define the metric coefficients g^ij and g_ij , which define the differential arc length along a curve in the curvilinear system

$$\begin{align} g_{ij} = \frac{\partial \textbf{R}}{\partial u^i}\cdot\frac{\partial \textbf{R}}{\partial u^j}\end{align} \tag{ C2 }$$

$$\begin{align} g^{ij} = \boldsymbol{\nabla}u^i\cdot\boldsymbol{\nabla}u^{\,j}\end{align} \tag{ C3 }$$

where R is the position in a Cartesian coordinate system and the uⁱ are the curvilinear coordinates. Hence, using GENE notation, $u^1 = x$, $u^2 = y$, $u^3 = z$.

We can also find the Jacobian of the transformation from Cartesian to curvilinear coordinates

$$\begin{align} \mathcal{J} = \frac{\partial \textbf{R}}{\partial u^1}\cdot\left[\frac{\partial \textbf{R}}{\partial u^2}\wedge\frac{\partial \textbf{R}}{\partial u^3}\right],\end{align} \tag{ C4 }$$

which is related to the metric coefficients through [34]

$$\begin{align} \sqrt{\det{\left[g^{ij}\right]}} = \mathcal{J}.\end{align} \tag{ C5 }$$

In this appendix we show how equation (10) can be derived.The magnetic drift frequency ω_D can be written as

$$\begin{align} \omega_\textrm{D}& = \textbf{v}_\textrm{D}\cdot\textbf{k}_\perp \nonumber\\ & = \left[\frac{1}{B\Omega}\left(\frac{v_\perp^2}{2}+v_\parallel^2\right)\textbf{b}\wedge\boldsymbol{\nabla}B\right]\cdot\left[ k_x\boldsymbol{\nabla}x+k_y\boldsymbol{\nabla}y\right] \nonumber\\ & = V\left[\textbf{b}\wedge\boldsymbol{\nabla}B\right]\cdot\left[ k_x\boldsymbol{\nabla}x+k_y\boldsymbol{\nabla}y\right], \end{align} \tag{ D1 }$$

where we used the definitions of $\textbf{v}_\textrm{D}$ and $\textbf{k}_\perp$, and we have hidden all the constants in V. Now, to make all the terms explicit, we can use the mixed product rule in curvilinear coordinates

$$\begin{align} \left(\textbf{A}\wedge\textbf{B}\right)\cdot\textbf{C} = \mathcal{J}\epsilon_{ijk}A^jB^kC^i,\end{align} \tag{ D2 }$$

where $\mathcal{J}$ is the Jacobian and we made use of Einstein's notation for repeated indices. As will be clear in the following, it is convenient to use the covariant version of B and C. Hence, using the contravariant metric tensor g^ij , we can write equation (D2) as

$$\begin{align} \left(\textbf{A}\wedge\textbf{B}\right)\cdot\textbf{C} = \mathcal{J}\epsilon_{ijk}g^{km}g^{il}A^jB_mC_l.\end{align} \tag{ D3 }$$

Then, we can set $\textbf{A} = \textbf{b}$ and use the fact that $\textbf{b}\parallel\textbf{z}$ (the coordinate parallel to the magnetic field lines). Hence, $b^j = b\delta^{j3}$, where δ^ij is Kronecker's delta. Therefore, we can simplify equation (D3)

$$\begin{align} \left(\textbf{b}\wedge\textbf{B}\right)\cdot\textbf{C} = \mathcal{J}b\epsilon_{i3k}g^{km}g^{il}B_mC_l,\end{align} \tag{ D4 }$$

which gives

$$\begin{align} \left(\textbf{b}\wedge\textbf{B}\right)\cdot\textbf{C} = \mathcal{J}b\left(\epsilon_{132}g^{2m}g^{1l}B_mC_l+\epsilon_{231}g^{1m}g^{2l}B_mC_l\right).\end{align} \tag{ D5 }$$

Hence, by applying the Levi-Civita tensor rule

$$\begin{align} \left(\textbf{b}\wedge\textbf{B}\right)\cdot\textbf{C} = \mathcal{J}b\left(g^{1m}g^{2l}-g^{2m}g^{1l}\right)B_mC_l.\end{align} \tag{ D6 }$$

Now, since the metric tensor g is symmetric, all the terms with l = m will cancel each other. Therefore, by expanding equation (D6), we can write

$$\begin{align} \left(\textbf{b}\wedge\textbf{B}\right)\cdot\textbf{C}& = \mathcal{J}b\left[\left(g^{11}g^{22}-g^{21}g^{12}\right)B_1C_2 \right.\nonumber\\ & \left. \quad + \left(g^{11}g^{23}-g^{21}g^{13}\right)B_1C_2 \right.\nonumber\\ & \left. \quad + \left(g^{12}g^{21}-g^{22}g^{11}\right)B_2C_1 \right.\nonumber\\ & \left. \quad + \left(g^{12}g^{23}-g^{22}g^{13}\right)B_2C_3 \right.\nonumber\\ & \left. \quad + \left(g^{13}g^{21}-g^{23}g^{11}\right)B_3C_1 \right.\nonumber\\ & \left. \quad + \left(g^{13}g^{22}-g^{23}g^{12}\right)B_3C_2\right]. \end{align} \tag{ D7 }$$

We can define some useful quantities to simplify the previous equation

$$\begin{align} \gamma_1 = g^{11}g^{22}-g^{21}g^{12}\end{align} \tag{ D8 }$$

$$\begin{align} \gamma_2 = g^{11}g^{23}-g^{21}g^{13},\end{align} \tag{ D9 }$$

$$\begin{align} \gamma_3 = g^{12}g^{23}-g^{22}g^{13}\end{align} \tag{ D10 }$$

which allow us to write equation (D7) as

$$\begin{align} \left(\textbf{b}\wedge\textbf{B}\right)\cdot\textbf{C}& = \mathcal{J}b\left[\gamma_1\left(B_1C_2-B_2C_1\right)+\left(\gamma_2B_1+\gamma_3B_2\right)C_3 \right.\nonumber\\ & \left. \quad -B_3\left(\gamma_2C_1+\gamma_3C_2\right)\right]. \end{align} \tag{ D11 }$$

Finally, if we replace $\textbf{B} = \boldsymbol{\nabla}B$ and $\textbf{C} = k_x\boldsymbol{\nabla}x+k_y\boldsymbol{\nabla}y$ in the previous equation, we can write equation (D1) as

$$\begin{align} \omega_\textrm{D} = \Lambda\left(k_x\underbrace{\left(-\partial_y B-\frac{\gamma_2}{\gamma_1}\partial_z B\right)}_{v_{Dx}}+k_y\underbrace{\left(\partial_x B -\frac{\gamma_3}{\gamma_1}\partial_z B\right)}_{v_{Dy}}\right)\end{align} \tag{ D12 }$$

where the scalar Λ comprises all the constants, v_Dx is the radial component of the magnetic drift velocity and v_Dy the binormal one.

In this appendix, we apply our NT ITG physical picture to elongation and magnetic shear scans performed for the large A pITG-2 scenario with adiabatic electrons. While in general it is necessary to use kinetic electrons to accurately model NT, previous simulations showed that adiabatic electrons are sufficient in this particular case.

E.1. Elongation scan

We took the NT and PT large A pITG-2 equilibria and performed NL simulations where the elongation was changed in a range from 0.5 to 2.0. Figure 25 shows the ratio between the nonlinear heat fluxes from the NT cases and the PT cases as a function of the elongation. We observe a monotonic decrease in the ratio with elongation. Therefore, larger values of κ increase the beneficial effect of NT on confinement and very low values of elongation make NT detrimental.

Zoom In Zoom Out Reset image size

To understand this behaviour, we can apply the NT ITG physical picture and explain the trend in terms of magnetic drifts and FLR effects. Figure 26 shows how the binormal component of the magnetic drifts v_Dy and the g^yy element of the metric tensor change as the elongation is varied. We observe that the behaviour of these two factors correlates well with the beneficial effect of NT. Indeed, as the elongation is decreased, the drift velocity of NT geometry varies less with poloidal angle. For κ = 0.5 it is actually closer to the resonance condition than PT around the outboard midplane. Similar behaviour can be observed for the FLR effect. At large elongation, the g^yy metric of the NT geometry is above the one of the PT case. However, when elongation is decreased below κ = 1, NT's g^yy is below PT's around the outboard midplane. Therefore, applying our NT ITG physical picture, we can say that the situation is the opposite when at κ = 0.5 compared to the rest of the paper: it is PT that has faster magnetic drifts and stronger FLR effects.

Zoom In Zoom Out Reset image size

E.2. Magnetic shear scan

We performed the same study as in the previous section, but with the magnetic shear $\hat{s}$ instead of elongation. Magnetic shear $\hat{s}$ was changed from large positive to negative values of shear, in the range $\hat{s} = [-0.5,2.0]$. Figure 27 shows the ratio between the nonlinear heat flux from the NT cases and the PT cases as a function of the magnetic shear. We observe a monotonic decrease of the ratio with $\hat{s}$. We could not find a value of magnetic shear for which NT becomes detrimental, but NT becomes less beneficial at lower $\hat{s}$. These observations are consistent with [11]. However, here we can try to apply our NT ITG physical picture to explain this trend.

Zoom In Zoom Out Reset image size

In figure 28 we show the binormal component of the magnetic drift v_Dy and the g^yy metric element as a function of the poloidal angle, for different values of magnetic shear. We observe that as the magnetic shear is decreased towards negative values, v_Dy becomes more similar between the two geometries, though NT remains slightly farther away from the resonance than PT. If we look at g^yy , we observe that the same is true for FLR effects—they become more similar as $\hat{s}$ is reduced, though FLR effects remain slightly stronger in NT at the lowest value of $\hat{s}$. Therefore, we can conclude that the beneficial effect of NT on ITG decreases at low magnetic shear because the drift velocity of ions gets closer to the resonance condition, and the FLR effects are reduced.

Zoom In Zoom Out Reset image size