1 Introduction
Controlling turbulence in magnetically confined plasmas is the key to achieving sustained nuclear fusion as an energy source (Krushelnick & Cowley Reference Krushelnick and Cowley2005). Typically, unstable perturbations driven by the pressure gradient and other sources of free energy grow exponentially and eventually saturate nonlinearly, leading to turbulence. Recent work (Newton, Cowley & Loureiro Reference Newton, Cowley and Loureiro2010; Highcock et al. Reference Highcock, Barnes, Schekochihin, Parra, Roach and Cowley2010, Reference Highcock, Barnes, Parra, Schekochihin, Roach and Cowley2011; Barnes et al. Reference Barnes, Parra, Highcock, Schekochihin, Cowley and Roach2011a ; Schekochihin, Highcock & Cowley Reference Schekochihin, Highcock and Cowley2012) has shown that in the presence of sheared flows, such systems can be subcritical. This means that all perturbations are linearly stable and a transition to a turbulent state only occurs if large enough initial perturbations undergo sufficient transient growth to allow nonlinear interaction. Understanding the transition to a turbulent state is a long-standing challenge in fluids (Barkley et al. Reference Barkley, Song, Mukund, Lemoult, Avila and Hof2015) and, more recently, in fusion plasmas, where a quiescent state leads to dramatically improved confinement. Experimental studies in simple devices (Weixing et al. Reference Weixing, Huang, Wang and Yu1993; Klinger et al. Reference Klinger, Latten, Piel, Bonhomme and Pierre1997; Riccardi et al. Reference Riccardi, Xuantong, Salierno, Gamberale and Fontanesi1997; Burin et al. Reference Burin, Tynan, Antar, Crocker and Holland2005) have proposed that this transition occurs through an increasing number of unstable frequencies leading to a turbulent state with a broadband spectrum. There is, however, currently very little known about a subcritical transition to turbulence in fusion-relevant plasmas. Here we use first-principles gyrokinetic simulations of a turbulent plasma in the outer core of the Mega-Ampere Spherical Tokamak (MAST) to demonstrate that the experimentally observed state is near the transition threshold, that the turbulence in this state is subcritical, and that transition to turbulence occurs via accumulation of long-lived, intense, finite-amplitude coherent structures, which dominate the near-threshold state. This represents a conceptually new and thus far unexplored scenario for transition to turbulence in magnetised plasmas.
The seemingly abrupt transition to turbulence from a quiescent state is a fluid-dynamical phenomenon that has fascinated scientists since the seminal experiments by Reynolds (Reference Reynolds1883). In fusion plasmas, transition to turbulence is a similarly tantalising challenge, both as a matter of fundamental physics and because understanding and controlling turbulence remains the greatest challenge to technologically and commercially effective fusion-power generation (Krushelnick & Cowley Reference Krushelnick and Cowley2005). This challenge arises from the turbulent eddies’ propensity for transporting heat and particles out of the core of the device, often leading to dramatically degraded plasma confinement (Kotschenreuther et al. Reference Kotschenreuther, Dorland, Beer and Hammett1995).
Extensive experimental (Baker et al.
Reference Baker, Greenfield, Burrell, DeBoo, Doyle, Groebner, Luce, Petty, Stallard and Thomas2001; Mantica et al.
Reference Mantica, Strintzi, Tala, Giroud, Johnson, Leggate, Lerche, Loarer, Peeters and Salmi2009; Field et al.
Reference Field, Michael, Akers, Candy, Colyer, Guttenfelder, Ghim, Roach and Saarelma2011; Ghim et al.
Reference Ghim, Field, Schekochihin, Highcock and Michael2014) and numerical (Kotschenreuther et al.
Reference Kotschenreuther, Dorland, Beer and Hammett1995; Dimits et al.
Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000; Roach et al.
Reference Roach, Abel, Akers, Arter, Barnes, Camenen, Casson, Colyer, Connor and Cowley2009; Barnes et al.
Reference Barnes, Parra, Highcock, Schekochihin, Cowley and Roach2011a
; Field et al.
Reference Field, Michael, Akers, Candy, Colyer, Guttenfelder, Ghim, Roach and Saarelma2011; Highcock et al.
Reference Highcock, Schekochihin, Cowley, Barnes, Parra, Roach and Dorland2012) work has identified the key parameters that trigger the transition to a turbulent regime at certain critical values. In particular, the ion-temperature gradient (ITG),
$\unicode[STIX]{x1D705}_{T}=-\text{d}\ln T_{i}/\text{d}r$
(
$T_{i}$
is the ion temperature and
$r$
is an appropriate dimensionless radial coordinate defined later), acts as a source of free energy (Coppi, Rosenbluth & Sagdeev Reference Coppi, Rosenbluth and Sagdeev1967), driving turbulent fluctuations, whereas differential toroidal rotation of the plasma, quantified by the (non-dimensionalised) flow shear perpendicular to the confining magnetic field,
$\unicode[STIX]{x1D6FE}_{E}=(r/q)\,\text{d}\unicode[STIX]{x1D714}/\text{d}r\,(a/v_{thi})$
(
$q$
is the ‘magnetic safety factor’,
$\unicode[STIX]{x1D714}$
is the angular frequency of toroidal rotation,
$a$
is the minor radius of the toroidal device and
$v_{thi}=\sqrt{2T_{i}/m_{i}}$
is the ion thermal velocity) can suppress turbulence (Burrell Reference Burrell1997).
In a steady state, a given amount of power injected (or, in a fusion-power plant, fusion generated) in the core of the device must be transported out through the plasma. Since the heat flux typically increases with the temperature gradient, more power requires maintaining a larger ITG. In the most straightforward scenario, above a certain threshold value of ITG, the plasma becomes linearly unstable, perturbations to the equilibrium are amplified and saturate, giving rise to a turbulent state. In this state, heat is transported by turbulent fluctuations, whose amplitude becomes larger at larger ITG (Barnes, Parra & Schekochihin Reference Barnes, Parra and Schekochihin2011b ). As a result, heat flux typically scales very strongly with the ITG, leading to ‘stiff transport’, with the practical consequence that the system cannot stray too far above the threshold value of ITG (Kotschenreuther et al. Reference Kotschenreuther, Dorland, Beer and Hammett1995).
The situation becomes more complicated in the presence of differential rotation. Perpendicular flow shear has been shown to have a suppressing effect on the linear instabilities and can even render the plasma completely linearly stable, i.e. all modes decay exponentially at large times. However, this may still entail substantial transient growth of perturbations (Roach et al. Reference Roach, Abel, Akers, Arter, Barnes, Camenen, Casson, Colyer, Connor and Cowley2009; Newton et al. Reference Newton, Cowley and Loureiro2010; Schekochihin et al. Reference Schekochihin, Highcock and Cowley2012) and, given finite initial perturbations, can lead to a saturated nonlinear state – a phenomenon known as ‘subcritical’ turbulence (Highcock et al. Reference Highcock, Barnes, Schekochihin, Parra, Roach and Cowley2010, Reference Highcock, Barnes, Parra, Schekochihin, Roach and Cowley2011; Barnes et al. Reference Barnes, Parra, Highcock, Schekochihin, Cowley and Roach2011a ). It is then an intriguing question how the heat flux due to such turbulence can change continuously from zero below the transition threshold to small but finite ITG-dependent values just above it (as it indeed does, in our simulations). In the case of supercritical turbulence, the saturated fluctuation amplitude everywhere increases continuously from small values near the threshold to finite ones far above it. In contrast, in the regime leading to subcritical turbulence, we find that small-amplitude perturbations decay and only finite-amplitude ones can survive and saturate. Therefore, the turbulent heat flux must increase with increasing ITG by some mechanism other than a continuous increase in fluctuation amplitude. Here we will identify this mechanism (which will lead us to a very different transition scenario than the conventional one outlined above) and ascertain that it is relevant to real experimental situations.
2 Gyrokinetic simulations
As an example of such a real experimental situation, we consider turbulence in the MAST tokamak, which is a major current experimental machine that is well diagnosed and actively used to test novel fusion concepts. We pick a magnetic configuration and plasma parameters describing the outer core of MAST for a particular discharge (no. 27268; see Field et al. Reference Field, Dunai, Ghim, Hill, McMillan, Roach, Saarelma, Schekochihin and Zoletnik2014 for its detailed description). We then solve numerically for the turbulent fluctuations in a local ‘flux tube’, by means of gyrokinetic simulations (Dimits et al. Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000; Fasoli et al. Reference Fasoli, Brunner, Cooper, Graves, Ricci, Sauter and Villard2016) with the widely used code GS2Footnote 1 . Our simulations are electrostatic, restricted to ion scales, include both kinetic ions and kinetic electrons and model collisions using a linearised Fokker–Planck collision operator (see appendix A for the specific equilibrium and resolution parameters used).
A GS2 flux tube is a twisted box of finite radial and poloidal width that follows a magnetic field line once around the tokamak in the poloidal direction (figure 1
a). We assume that the underlying equilibrium is axisymmetric, so a single flux tube can be used to simulate the entire flux surface, allowing a large saving in computational cost. In this local approach, all equilibrium gradients (of density, velocity and temperature) are assumed constant across the radial extent of the flux tube. We pick the time
$t=0.25$
s from the beginning of the discharge and the radial location
$r=D/2a=0.8$
, where
$D$
is the diameter of the flux surface hosting our flux tube, at the height of the magnetic axis and
$2a$
is the diameter of the last closed flux surface (the ‘edge’ of the plasma). This definition of the radial location is convenient because it coincides with the Miller et al. (Reference Miller, Chu, Greene, Lin-Liu and Waltz1998) parametrisation of the flux-tube geometry used by GS2. Figure 1(b) shows a poloidal projection of the MAST flux surfaces for discharge no. 27268 with the flux surface at
$r=0.8$
highlighted.
(a) The flux surface at
$r=0.8$
traced by the field lines that lie in it. The field line marked in red is the centre line of a GS2 ‘flux tube’. The actual GS2 flux tube is rectangular in the outboard midplane but, as the magnetic field is sheared, the flux tube twists as it follows the field line along the flux surface. (b) Poloidal projection of the MAST flux surfaces. The flux surface at
$r=0.8$
, the last closed flux surface and the separatrix are indicated. The Miller et al. (Reference Miller, Chu, Greene, Lin-Liu and Waltz1998) radial coordinate
$r$
of a flux surface is the ratio of the diameter of this flux surface at the elevation of the magnetic axis (
$Z=0$
in this plot) to that of the last closed flux surface.
The gyrokinetic equation (Frieman & Chen Reference Frieman and Chen1982; see appendix A) solved by GS2 gives us the time evolution of the perturbed distribution function
$\unicode[STIX]{x1D6FF}f_{s}(t,\boldsymbol{r},\boldsymbol{v})$
of particles of species
$s$
(ions or electrons), where
$\boldsymbol{r}$
is the spatial position and
$\boldsymbol{v}$
the velocity. Gyrokinetic theory makes use of the fact that turbulent fluctuations in a tokamak plasma occur at much longer time scales than the Larmor motion of the particles (although still much shorter than the evolution of the background thermal and magnetic equilibrium state). This allows the Larmor motion to be averaged over analytically, leading to a kinetics of Larmor rings whose distribution depends on two, rather than three, velocity variables (parallel
$v_{\Vert }$
and perpendicular
$v_{\bot }$
velocities but not the gyroangle), reducing the phase space from six to five dimensions – a more tractable problem than the Vlasov–Boltzmann equation for the evolution of the full distribution function.
3 Turbulent heat flux
With the knowledge of the distribution function, one can calculate any characteristics of the turbulence. In particular, the turbulent heat flux as a function of the temperature gradient, flow shear or any other equilibrium parameters is a key quantity of interest. Focusing on the heat flux will allow us to diagnose the transition to turbulence, as values greater than zero indicate a turbulent state. We focus here on the heat flux due to the ions (deuterium in our simulations). The radially outward, time-averaged ion heat (energy) flux through a volume
$V$
enveloping a given flux surface is
$$\begin{eqnarray}Q_{i}=\left\langle \frac{1}{V}\int \text{d}^{3}\boldsymbol{r}\int \text{d}^{3}\boldsymbol{v}\frac{m_{i}v^{2}}{2}\unicode[STIX]{x1D6FF}f_{i}\boldsymbol{V}_{E}\boldsymbol{\cdot }\unicode[STIX]{x1D735}r\right\rangle ,\end{eqnarray}$$
where
$\boldsymbol{V}_{E}$
is the
$\boldsymbol{E}\times \boldsymbol{B}$
drift velocity due to the perturbed electric field, calculated from
$\unicode[STIX]{x1D6FF}f_{s}$
via the plasma quasineutrality constraint (see appendix A) and
$\left<\cdots \,\right>$
indicates an average in time.
$Q_{i}$
is typically normalised by the so-called gyro-Bohm value,
$Q_{gB}=n_{i}T_{i}v_{thi}\unicode[STIX]{x1D70C}_{i}^{2}/a^{2}$
, where
$n_{i}$
,
$T_{i}$
and
$\unicode[STIX]{x1D70C}_{i}$
are the ion density, temperature and Larmor radius, respectively (it is a feature of the asymptotic ordering on which the gyrokinetic theory is based that
$Q_{i}/Q_{gB}$
is a finite number; see Abel et al.
Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013).
To map out the transition to turbulence in our system, we vary the flow shear
$\unicode[STIX]{x1D6FE}_{E}$
and the temperature gradient
$\unicode[STIX]{x1D705}_{T}$
around their experimental values (
$\unicode[STIX]{x1D6FE}_{E}=0.16\pm 0.02$
and
$\unicode[STIX]{x1D705}_{T}=5.1\pm 1$
) and covering the range
$\unicode[STIX]{x1D6FE}_{E}\in [0,0.19]$
and
$\unicode[STIX]{x1D705}_{T}\in [4.3,8.0]$
. Figure 2 shows the turbulent heat flux
$Q_{i}/Q_{gB}$
in a part of this range, close to the threshold. All simulations were run until they reached a statistically steady state, i.e. until the running time-average value became independent of time. An average was then taken over a period of
${\sim}100(a/v_{thi})$
(which represents
${\sim}400~\unicode[STIX]{x03BC}$
s) during this steady state. Examining both the range of values of
$(\unicode[STIX]{x1D6FE}_{E},\unicode[STIX]{x1D705}_{T})$
compatible with experiment and the experimentally determined value
$Q_{i,\mathit{exp}}/Q_{gB}\approx 2.0\equiv Q_{\mathit{exp}}$
(Field et al.
Reference Field, Dunai, Ghim, Hill, McMillan, Roach, Saarelma, Schekochihin and Zoletnik2014), we see that the turbulent state found in a real device is close to the threshold (perhaps not a surprising conclusion, but an important one to be able to make quantitatively). Experimental investigations (e.g. Mantica et al.
Reference Mantica, Strintzi, Tala, Giroud, Johnson, Leggate, Lerche, Loarer, Peeters and Salmi2009) corroborate this observation for other tokamaks and show that the proximity to threshold is enforced by the rapid increase in
$Q_{i}$
as the stability parameter
$\unicode[STIX]{x1D705}_{T}$
is increased from its critical value (‘stiff transport’). A similar conclusion can be drawn from figure 2: small increases in
$\unicode[STIX]{x1D705}_{T}$
lead to order of magnitude changes in
$Q_{i}$
. We find that very small departures of the flow shear
$\unicode[STIX]{x1D6FE}_{E}$
from the threshold also lead to large increases of
$Q_{i}$
, showing that flow shear matters at experimentally relevant values and heat transport is highly sensitive to it.
Normalised turbulent heat flux as a function of flow shear
$\unicode[STIX]{x1D6FE}_{E}$
and temperature gradient
$\unicode[STIX]{x1D705}_{T}$
. The marked region indicates the range of
$(\unicode[STIX]{x1D6FE}_{E},\unicode[STIX]{x1D705}_{T})$
compatible with the experimentally measured values on our chosen flux surface and the dashed line is the contour corresponding to the experimental value
$Q_{i}/Q_{gB}\approx 2$
(Field et al.
Reference Field, Dunai, Ghim, Hill, McMillan, Roach, Saarelma, Schekochihin and Zoletnik2014). Thus, the experiment is close to the turbulence threshold. The two black dots mark the cases for which the density-fluctuation field is shown in figure 4.
(a) Transient growth of initial perturbations of the electrostatic potential
$\unicode[STIX]{x1D711}$
for
$\unicode[STIX]{x1D6FE}_{E}=0.16$
(experimentally measured value) and several values of the ion-temperature gradient
$\unicode[STIX]{x1D705}_{T}$
. These time-evolution curves are obtained in purely linear simulations, for a single mode poloidal wavenumber
$k_{y}\unicode[STIX]{x1D70C}_{i}=0.2$
, which is approximately the wavenumber of the most vigorous transient growth. The radial wavenumber,
$k_{r}$
, is allowed to vary with time to resolve the effect of flow shear that causes growing modes to be shifted to neighbouring wavenumbers. As the stability parameter
$\unicode[STIX]{x1D705}_{T}$
increases away from the (nonlinear) threshold (cf. figure 2), initial perturbations are amplified by an ever larger factor before decaying. (b) Illustration of subcritical turbulent state (as measured by the normalised heat flux,
$Q_{i}/Q_{gB}$
) being reached starting from a finite perturbation, whereas small perturbations decay. The three nonlinear simulations shown here differ only in their initial perturbation level, all other parameters being the same and lying within the window of experimentally consistent values (see figure 2),
$\unicode[STIX]{x1D6FE}_{E}=0.16$
and
$\unicode[STIX]{x1D705}_{T}=5.5$
(so the linear transient amplification in this case is very low; see (a)). The time histories start at
$t=10$
because we excluded some initial evolution involving various numerical adjustments.
4 Subcritical turbulence
Usually one would also carry out a sequence of linear simulations to ascertain whether the turbulence threshold found nonlinearly coincides with the linear stability boundary. Doing these simulations showed that all modes in our system at all parameter values that we have investigated (except
$\unicode[STIX]{x1D6FE}_{E}=0$
) were formally linearly stable. Initial perturbations did, however, exhibit transient growth, typically for a longer period in the cases far from the nonlinear threshold, as illustrated in figure 3(a) (cf. Highcock et al.
Reference Highcock, Barnes, Schekochihin, Parra, Roach and Cowley2010, Newton et al.
Reference Newton, Cowley and Loureiro2010, Barnes et al.
Reference Barnes, Parra, Highcock, Schekochihin, Cowley and Roach2011a
, Highcock et al.
Reference Highcock, Schekochihin, Cowley, Barnes, Parra, Roach and Dorland2012, Schekochihin et al.
Reference Schekochihin, Highcock and Cowley2012). Nonlinearly, this means that, beyond a certain threshold in
$\unicode[STIX]{x1D6FE}_{E}$
and
$\unicode[STIX]{x1D705}_{T}$
, and given a large enough initial perturbation, subcritical turbulence can be sustained.Footnote
2
This is illustrated in figure 3(b), showing the effect of changing the initial perturbation amplitude. There is clearly a critical amplitude above which the nonlinearity can pick up the transiently amplified perturbations (very weakly amplified, when close to threshold) and give rise to a non-zero saturated state. Importantly, the saturation level does not depend on the size of the initial perturbation (as long as the latter is large enough). Thus, in the experimental instance that we have considered, ion-scale turbulence in MAST in the presence of flow shear is subcritical and so tokamak plasmas join a plethora of neutral fluid systems where the transition to turbulence depends strongly on the (size of) initial perturbation (Trefethen et al.
Reference Trefethen, Trefethen, Reddy and Driscoll1993; Darbyshire & Mullin Reference Darbyshire and Mullin1995): e.g. both Poiseuille and Couette flows are formally stable (Salwen, Cotton & Grosch Reference Salwen, Cotton and Grosch1980; Trefethen et al.
Reference Trefethen, Trefethen, Reddy and Driscoll1993), but still able to transition to a nonlinear, turbulent state; a similar situation arises in Keplerian shear flows believed to exist in accretion disks (Riols et al.
Reference Riols, Rincon, Cossu, Lesur, Longaretti, Ogilvie and Herault2013). Recent theoretical work, involving very simple models, suggested that this may also be possible in plasmas (Newton et al.
Reference Newton, Cowley and Loureiro2010; Schekochihin et al.
Reference Schekochihin, Highcock and Cowley2012; Landreman, Plunk & Dorland Reference Landreman, Plunk and Dorland2015), as did simulations of simplified tokamak equilibria (Barnes et al.
Reference Barnes, Parra, Highcock, Schekochihin, Cowley and Roach2011a
; Highcock et al.
Reference Highcock, Barnes, Schekochihin, Parra, Roach and Cowley2010, Reference Highcock, Barnes, Parra, Schekochihin, Roach and Cowley2011, Reference Highcock, Schekochihin, Cowley, Barnes, Parra, Roach and Dorland2012), but ours appears to be the first demonstration of subcritical ion-scale turbulence in a specific, experimentally diagnosed tokamak plasma.
5 Scenario for transition to turbulence
How does the transition to subcritical turbulence occur, i.e. what sequence of turbulent states does the system go through as either
$\unicode[STIX]{x1D6FE}_{E}$
or
$\unicode[STIX]{x1D705}_{T}$
crosses the critical threshold and moves away from it into ever more strongly driven regimes? It is clear that the transition cannot occur via near-threshold states featuring arbitrarily small fluctuations everywhere, because sustaining subcritical turbulence requires finite initial perturbations. These initial perturbations must be larger near the threshold than far from it because the amount of amplification expected during transient growth tends to decrease close to the threshold (Highcock et al.
Reference Highcock, Schekochihin, Cowley, Barnes, Parra, Roach and Dorland2012; Schekochihin et al.
Reference Schekochihin, Highcock and Cowley2012). If the typical maximum amplitude of the fluctuations in the saturated state must remain finite, one way to reduce the turbulent heat flux to low values near the threshold is by reducing the fraction of the system’s volume taken up by turbulence, i.e. by concentrating intense fluctuations in a shrinking part of space. This is precisely what happens, as we will now demonstrate.
In figure 4, we show real-space snapshots of the turbulent density-fluctuation field in the poloidal cross-section of our flux tube at
$\unicode[STIX]{x1D6FE}_{E}=0.16$
and for two different temperature gradients:
$\unicode[STIX]{x1D705}_{T}=4.8$
, which is very close to the threshold, and
$\unicode[STIX]{x1D705}_{T}=5.2$
, a case that represents more strongly driven turbulence away from the threshold (both points are within the experimentally consistent range; see figure 2). We find that the near-threshold turbulent state is dominated by long-lived, intense coherent structures, which travel across the domain both radially and poloidally, whereas far from the threshold, we observe a more conventional chaotic turbulent state characterised by interacting eddies. These two cases are representative of the relevant regions of our parameter space. We always find that long-lived, large-amplitude structures form in the near-threshold cases and survive against the background of very weak ambient fluctuations. As the system is taken away from the threshold by increasing
$\unicode[STIX]{x1D705}_{T}$
or decreasing
$\unicode[STIX]{x1D6FE}_{E}$
, these structures become more numerous (i.e. more volume filling) while retaining comparable amplitude, eventually start interacting with each other, and break up. Finally, far from the threshold, we observe no discernible long-lived structures, but rather strong time-variable fluctuations everywhere with amplitudes that increase with
$\unicode[STIX]{x1D705}_{T}$
or decreasing
$\unicode[STIX]{x1D6FE}_{E}$
.
Density-fluctuation field
$\unicode[STIX]{x1D6FF}n_{i}/n_{i}$
in the poloidal plane for simulations close to and far from the turbulence threshold (the two cases marked by points in figure 2): (a)
$\unicode[STIX]{x1D6FE}_{E}=0.16$
and
$\unicode[STIX]{x1D705}_{T}=4.8$
(near threshold), (b)
$\unicode[STIX]{x1D6FE}_{E}=0.16$
and
$\unicode[STIX]{x1D705}_{T}=5.2$
(strongly driven). Half of the full simulation domain in the vertical direction is shown.
Let us make the statements above more quantitative. Consider how the maximum amplitude
$A_{max}$
of the density perturbations found across the domain (and averaged over time) changes as our stability parameters
$\unicode[STIX]{x1D705}_{T}$
and
$\unicode[STIX]{x1D6FE}_{E}$
change. Since
$Q_{i}$
is a strong function of both
$\unicode[STIX]{x1D705}_{T}$
and
$\unicode[STIX]{x1D6FE}_{E}$
, we can measure the ‘distance from threshold’ by just using
$Q_{i}/Q_{gB}$
as a stability parameter. A naive estimate based on (3.1) is
$$\begin{eqnarray}\frac{Q_{i}}{Q_{gB}}\sim \frac{a^{2}}{\unicode[STIX]{x1D70C}_{i}^{2}}\frac{\unicode[STIX]{x1D6FF}n_{i}}{n_{i}}\frac{V_{Er}}{v_{thi}}\sim k_{y}\unicode[STIX]{x1D70C}_{i}\frac{T_{e}}{T_{i}}\left(\frac{a}{\unicode[STIX]{x1D70C}_{i}}\frac{\unicode[STIX]{x1D6FF}n_{i}}{n_{i}}\right)^{2}\sim A^{2},\end{eqnarray}$$
where
$A=(a/\unicode[STIX]{x1D70C}_{i})\unicode[STIX]{x1D6FF}n_{i}/n_{i}$
(which in the gyrokinetic theory is an order-unity quantity; see Abel et al.
Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013). We have estimated the radial
$\boldsymbol{E}\times \boldsymbol{B}$
velocity as
$V_{Er}\sim (c/B)k_{y}\unicode[STIX]{x1D711}\sim k_{y}\unicode[STIX]{x1D70C}_{i}v_{thi}e\unicode[STIX]{x1D711}/T_{i}$
, where
$k_{y}$
is the typical poloidal wavenumber (
${\sim}\unicode[STIX]{x1D70C}_{i}^{-1}$
in this regime) of the fluctuations of the electrostatic potential
$\unicode[STIX]{x1D711}$
, which are related (by order of magnitude) to the electron (and, therefore, ion) density via the Boltzmann response
$e\unicode[STIX]{x1D711}/T_{e}\sim \unicode[STIX]{x1D6FF}n_{e}/n_{e}$
. Figure 5(a) shows the relationship between
$A_{max}$
and
$Q_{i}/Q_{gB}$
for a number of simulations with different values of
$\unicode[STIX]{x1D705}_{T}$
and
$\unicode[STIX]{x1D6FE}_{E}$
. While the naive scaling (5.1) is indeed manifest far from the threshold, it is a striking feature of figure 5(a) that the maximum fluctuation amplitude hits a finite ‘floor’ as
$Q_{i}/Q_{gB}$
decreases to and below its experimental value – this coincides with the appearance of long-lived structures illustrated in figure 4(a). Thus, while in the conventional supercritical turbulence, we might have observed smaller fluctuation amplitudes corresponding to lower heat fluxes all the way to the threshold, in the present subcritical turbulent system, we see the heat flux decrease while the maximum fluctuation amplitude remains constant.
(a) Maximum amplitude of the density fluctuations versus the time-averaged ion heat flux. The naive scaling (5.1) is shown and holds far from threshold, whereas for small values of
$Q_{i}/Q_{gB}$
(around and below the experimental value
$Q_{\mathit{exp}}$
), the amplitude flattens. (b) Number of turbulent structures (amplitude within 75 % of the maximum) versus the time-averaged ion heat flux. It grows up to and slightly beyond the experimental value
$Q_{\mathit{exp}}$
, eventually the volume is filled with turbulent structures and their number tends to a constant. The scaling
$Q\propto N$
is shown for reference: the heat flux near threshold is controlled by the accumulation of distinct structures; far from the threshold, the number of structures is simply set by their scale and the size of the domain, whereas the heat flux is controlled by the overall fluctuation amplitude (see (a)).
As we argued above, it does this via reduction of the volume taken up by large fluctuations. We will now show this by measuring the typical number of the turbulent structures as a function of distance to threshold. While two-dimensional structures are easily discerned by a human eye (e.g. in figure 4(a) there are two), counting them systematically is a not entirely trivial problem, which is often encountered in computer-vision and pattern-recognition applications. It has been considered before in the context of experimental measurements of turbulence (Kauschke Reference Kauschke1999, Müller et al. Reference Müller, Fasoli, Labit, McGrath, Pisaturo, Plyushchev, Podestà and Poli2005; see review of various techniques by Love & Kamath Reference Love and Kamath2007). Structure counting can be reduced to an image-labelling, or segmentation, problem by applying a threshold function to our density fluctuations: setting values below a certain percentile (here 75 %) of the maximum amplitude to 0 and above it to 1. We are then left with an array of 1’s representing our structures against a background of 0’s. We employ a general-purpose image processing package scikit-image (van der Walt et al. Reference van der Walt, Schönberger, Nunez-Iglesias, Boulogne, Warner, Yager, Gouillart and Yu2014), which implements an efficient labelling algorithm (Fiorio & Gustedt Reference Fiorio and Gustedt1996), to label connected regions, i.e. turbulent structures. In order to improve the reliability of the labelling algorithm, we applied a Gaussian image filter (with a standard deviation of the order of the grid scale) as a pre-processing step and also removed structures below 10 % of the mean structure size as a post-processing step. These steps are justified because we are hunting intense, relatively large-scale structures.
Figure 5(b) shows the results of the above analysis: the number
$N$
of turbulent structures with amplitudes above the 75
$\text{th}$
percentile versus the heat flux. As in figure 5(a), there are two distinct regimes:
$N$
grows with
$Q_{i}$
until the structures have filled the simulation domain (which happens just beyond the experimental value of the flux), whereupon
$N$
tends to a constant. Taking figures 5(a) and 5(b) in combination, we have, approximately,
$Q_{i}/Q_{gB}\sim NA^{2}$
, i.e. near the threshold, turbulent heat flux increases because coherent structures become more numerous (but not more intense:
$N$
grows,
$A$
stays constant), whereas far from the threshold, it does so because the fluctuation amplitude
$A$
increases (while
$N$
stays constant).
6 Conclusion
In conclusion, we have discovered, using numerical simulations of an experimentally relevant fusion plasma, a novel scenario for transition to turbulent state – which applies when the turbulence is subcritical. Above a certain critical value of
$\unicode[STIX]{x1D705}_{T}$
and below a critical value of
$\unicode[STIX]{x1D6FE}_{E}$
, a large enough initial perturbation will ignite turbulence. Near the threshold, the density and temperature fluctuations (and hence heat flux) are concentrated in long-lived, intense coherent structures (interestingly, this is reminiscent of the transition via localised patches in pipe flows; see Barkley et al.
Reference Barkley, Song, Mukund, Lemoult, Avila and Hof2015, or transition to Magnetorotational instability turbulence in Keplerian accretion flows; see Riols et al.
Reference Riols, Rincon, Cossu, Lesur, Longaretti, Ogilvie and Herault2013, Reference Riols, Rincon, Cossu, Lesur, Ogilvie and Longaretti2016). As the stability parameters
$(\unicode[STIX]{x1D705}_{T},\unicode[STIX]{x1D6FE}_{E})$
depart slightly from their critical values into the more strongly driven regime, the number of these structures increases rapidly while their amplitude stays roughly constant (in contrast to the conventional supercritical turbulence, where the amplitude increases with
$\unicode[STIX]{x1D705}_{T}$
because arbitrarily low-amplitude turbulence can be supported close to threshold). Increasing the turbulent drive further leads to the turbulent structures filling the simulation domain and any further increase in the heat flux is caused by an increase in turbulent amplitude. The latter regime is similar to the conventional plasma turbulence.
It is the presence of flow shear that appears to be the key feature that allows tokamak turbulence to exist in two distinct regimes – more strictly speaking, these are two regimes of anomalous transport, rather than turbulence: it is not obvious that the structure-dominated near-threshold state can be viewed as properly turbulent, representing perhaps a persistent nonlinear solution rather than full-scale chaos (cf. Riols et al. Reference Riols, Rincon, Cossu, Lesur, Ogilvie and Longaretti2016). It will be very interesting to see if a structure-dominated regime can be detected in MAST or in other tokamaks where significant flow shear is present: so far, there are some tentative but encouraging indications that such a regime might manifest itself in experimentally observed skewed probability distributions of density fluctuations (Fox et al. Reference Fox, van Wyk, Field, Ghim, Parra and Schekochihin2016).
The new regime of tokamak turbulence described above, besides raising interesting questions of theory of the subcritical transition and its degree of universality, also raises potentially useful ones regarding ways (e.g. optimal combinations of momentum and power input) in which such a turbulence could be controlled and the associated heat flux further reduced, leading to better confined plasmas.
Acknowledgements
We would like to thank M. Barnes, J. Ball, G. Colyer, M. Fox and F. Parra for many useful discussions. The simulations were carried out in part using the HELIOS supercomputer system at International Fusion Energy Research Centre, Aomori, Japan, under the Broader Approach collaboration between Euratom and Japan, implemented by Fusion for Energy and JAEA. Simulations carried out using the ARCHER UK National Supercomputing Service (http://www.archer.ac.uk) were provided by the Plasma HEC Consortium (EP/L000237/1) and the Collaborative Computational Project in Plasma Physics funded by UK EPSRC (EP/M022463/1). This work has been carried out within the framework of the EUROfusion Consortium and received funding from the European Union’s Horizon 2020 research and innovation programme and Euratom research and training programme 2014–2018 under grant agreement no. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. The work of A.A.S. was supported in part by grants from UK EPSRC and STFC.
Appendix A. Gyrokinetic model and numerical set-up
Here we outline the theoretical framework that we have used for modelling plasma turbulence and give all the information necessary to reproduce our GS2 simulations.
The gyrokinetic equation (Frieman & Chen Reference Frieman and Chen1982; Abel et al.
Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013) describes the time evolution of the perturbed (from a background Maxwellian
$F_{s}$
) distribution function
$\unicode[STIX]{x1D6FF}f_{s}$
of particles of species
$s$
. Here we solve an approximate form of this equation arising from assuming, formally, that the Mach number
$M$
of the plasma rotational flow is small, but that the flow shear is large enough to matter: namely, if
$\unicode[STIX]{x1D714}$
is the angular frequency of toroidal rotation, then
$$\begin{eqnarray}\frac{R\unicode[STIX]{x1D714}}{v_{thi}}=M\ll 1,\quad |a\unicode[STIX]{x1D735}\ln \unicode[STIX]{x1D714}|\sim \frac{1}{M},\end{eqnarray}$$
where
$R$
is the major radius of the toroidal device at the location of our flux tube. This ordering allows one to formulate local gyrokinetics in a rotating flux surface, neglecting such effects as Coriolis and centrifugal forces, but retaining flow shear. With this ordering and assuming also purely electrostatic perturbations (no fluctuating magnetic fields), the gyrokinetic system of equations is written as follows (see § 11 of Abel et al.
Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013 or appendix A of Schekochihin et al.
Reference Schekochihin, Highcock and Cowley2012). The perturbed distribution function is split into a part corresponding to Boltzmann density response and the distribution of the gyrocentres:
$$\begin{eqnarray}\unicode[STIX]{x1D6FF}f_{s}=-\frac{Z_{s}e\unicode[STIX]{x1D711}(\boldsymbol{r})}{T_{s}}F_{s}+h_{s}(\boldsymbol{R},\unicode[STIX]{x1D700},\unicode[STIX]{x1D707},\unicode[STIX]{x1D70E}),\quad F_{s}=n_{s}(\unicode[STIX]{x1D713})\left[\frac{m_{s}}{2\unicode[STIX]{x03C0}T_{s}(\unicode[STIX]{x1D713})}\right]^{3/2}\text{e}^{-\unicode[STIX]{x1D700}/T_{s}(\unicode[STIX]{x1D713})},\end{eqnarray}$$
where
$Z_{s}e$
,
$m_{s}$
,
$n_{s}$
,
$T_{s}$
are the charge, mass, density, and temperature of particles of species
$s$
,
$\unicode[STIX]{x1D711}$
is the electrostatic potential perturbation,
$\boldsymbol{R}=\boldsymbol{r}-\boldsymbol{b}\times \boldsymbol{v}_{\bot }/\unicode[STIX]{x1D6FA}_{s}$
is the position of the centre of a particle’s Larmor orbit,
$\boldsymbol{r}$
is the position of the particle,
$\boldsymbol{b}=\boldsymbol{B}/B$
is a unit vector in the direction of the magnetic field
$B$
,
$\boldsymbol{v}_{\bot }$
is the velocity of the particle perpendicular to the magnetic field,
$\unicode[STIX]{x1D6FA}_{s}=Z_{s}eB/m_{s}c$
is the cyclotron frequency,
$c$
is the speed of light,
$\unicode[STIX]{x1D700}=m_{s}v^{2}/2$
the particle energy,
$\unicode[STIX]{x1D707}=m_{s}v_{\bot }^{2}/2B$
its magnetic moment and
$\unicode[STIX]{x1D70E}$
the sign of its parallel velocity
$v_{\Vert }=\pm [2(\unicode[STIX]{x1D700}-\unicode[STIX]{x1D707}B)/m_{s}]^{1/2}$
(
$\unicode[STIX]{x1D700}$
,
$\unicode[STIX]{x1D707}$
and
$\unicode[STIX]{x1D70E}$
are the velocity-space variables used by the GS2 code); the velocities are taken in the frame rotating with the angular frequency
$\unicode[STIX]{x1D714}(\unicode[STIX]{x1D713})$
, which, like
$n_{s}(\unicode[STIX]{x1D713})$
and
$T_{s}(\unicode[STIX]{x1D713})$
, is a function of the poloidal flux
$\unicode[STIX]{x1D713}$
only. The evolution equation for
$h_{s}$
is then
$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\right)\left(h_{s}-\frac{Z_{s}e\langle \unicode[STIX]{x1D711}\rangle _{R}}{T_{s}}F_{s}\right)+\left(v_{\Vert }\boldsymbol{b}+\boldsymbol{V}_{Ds}+\langle \boldsymbol{V}_{E}\rangle _{R}\right)\boldsymbol{\cdot }\unicode[STIX]{x1D735}h_{s}-\langle C[h_{s}]\rangle _{R}\nonumber\\ \displaystyle & & \displaystyle \displaystyle \qquad =-\langle \boldsymbol{V}_{E}\rangle _{R}\boldsymbol{\cdot }\unicode[STIX]{x1D735}r\left[\frac{\text{d}\ln n_{s}}{\text{d}r}+\left(\frac{\unicode[STIX]{x1D700}}{T_{s}}-\frac{3}{2}\right)\frac{\text{d}\ln T_{s}}{\text{d}r}+\frac{m_{s}v_{\Vert }}{T_{s}}\frac{RB_{\unicode[STIX]{x1D719}}}{B}\frac{\text{d}\unicode[STIX]{x1D714}}{\text{d}r}\right]F_{s},\end{eqnarray}$$
where
$\boldsymbol{u}=\unicode[STIX]{x1D714}(\unicode[STIX]{x1D713})R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}$
, with
$\unicode[STIX]{x1D719}$
being the toroidal angle, is the toroidal rotation velocity,
$$\begin{eqnarray}\boldsymbol{V}_{Ds}=\frac{c}{Z_{s}eB}\boldsymbol{b}\times \left[m_{s}v_{\Vert }^{2}\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{b}+\unicode[STIX]{x1D707}\unicode[STIX]{x1D735}B\right],\end{eqnarray}$$
is the magnetic drift velocity,
$$\begin{eqnarray}\boldsymbol{V}_{E}=\frac{c}{B}\boldsymbol{b}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D711}\end{eqnarray}$$
is the
$\boldsymbol{E}\times \boldsymbol{B}$
drift velocity,
$\langle \,\cdots \,\rangle _{R}$
is an average over a particle orbit at constant guiding-centre position
$\boldsymbol{R}$
,
$C[h_{s}]$
is the linearised collision operator,
$B_{\unicode[STIX]{x1D719}}$
is the toroidal component of the magnetic field. To calculate
$\unicode[STIX]{x1D711}$
in (A 3), the quasineutrality condition is used:
$$\begin{eqnarray}\mathop{\sum }_{s}Z_{s}\unicode[STIX]{x1D6FF}n_{s}=0\Rightarrow \mathop{\sum }_{s}\frac{Z_{s}^{2}e\unicode[STIX]{x1D711}}{T_{s}}n_{s}=\mathop{\sum }_{s}Z_{s}\int \text{d}^{3}\boldsymbol{v}\langle h_{s}\rangle _{r},\end{eqnarray}$$
where
$\langle \,\cdots \,\rangle _{r}$
means a gyroaverage at constant
$\boldsymbol{r}$
.
The last term in (A 3) represents the advection by the gyro-averaged
$\boldsymbol{E}\times \boldsymbol{B}$
velocity of the Maxwellian plasma equilibrium distribution characterised by
$n_{s}$
,
$T_{s}$
and
$\unicode[STIX]{x1D714}$
. These are functions only of the flux surface, conventionally labelled by the poloidal flux
$\unicode[STIX]{x1D713}$
, but here this dependence has been converted to the Miller et al. (Reference Miller, Chu, Greene, Lin-Liu and Waltz1998) radial coordinate
$r=D/2a$
, where
$D$
is the diameter of the flux surface of interest at the elevation of the magnetic axis and
$2a$
is the diameter of the last closed flux surface (see figure 1
b). Since
$r$
is also a flux-surface label,
$\unicode[STIX]{x1D735}r=(\text{d}\unicode[STIX]{x1D713}/\text{d}r)^{-1}\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$
. The gradients that appear in the right-hand side of (A 3) are sources of free energy in the plasma. In local flux-tube calculations, these are approximated as constant parameters and the following definitions for them are introduced:
$$\begin{eqnarray}-\frac{\text{d}\ln n_{s}}{\text{d}r}=\unicode[STIX]{x1D705}_{ns},\quad -\frac{\text{d}\ln T_{s}}{\text{d}r}=\unicode[STIX]{x1D705}_{Ts},\quad \frac{a}{v_{thi}}\frac{r}{q}\frac{\text{d}\unicode[STIX]{x1D714}}{\text{d}r}=\unicode[STIX]{x1D6FE}_{E},\end{eqnarray}$$
where
$q(\unicode[STIX]{x1D713})$
is the ‘safety factor’. Since this quantity in a tokamak is, approximately,
$q\sim (r/R)B/B_{p}$
, where
$B_{p}=|\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}|/R$
is the poloidal component of the magnetic field,
$\unicode[STIX]{x1D6FE}_{E}$
has the meaning of the (non-dimensionalised) part of the toroidal flow shear that is perpendicular to the local magnetic field. The special relevance of this quantity becomes obvious if we ‘unpack’ what is meant by
$\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}$
in the left-hand side of (A 3). Since GS2 solves the gyrokinetic equation locally in the vicinity of a particular flux surface
$\unicode[STIX]{x1D713}_{0}$
, we may expand, inside our flux tube,
$\unicode[STIX]{x1D714}\approx \unicode[STIX]{x1D714}_{0}+xB_{p}R\,\text{d}\unicode[STIX]{x1D714}/\text{d}\unicode[STIX]{x1D713}$
, where
$x$
is the distance from the flux surface. ThenFootnote
3
$$\begin{eqnarray}\boldsymbol{u}=\unicode[STIX]{x1D714}R^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}\approx \left(\unicode[STIX]{x1D714}_{0}R+xB_{p}R^{2}\frac{\text{d}\unicode[STIX]{x1D714}}{\text{d}\unicode[STIX]{x1D713}}\right)\left(\frac{B_{\unicode[STIX]{x1D719}}}{B}\boldsymbol{b}+\frac{B_{p}}{B}\boldsymbol{e}_{y}\right),\end{eqnarray}$$
where
$\boldsymbol{e}_{y}=\boldsymbol{b}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}/(B_{p}R)$
is the unit vector in the direction perpendicular to the field line but tangent to the flux surface. If we now go to the frame rotating with the flux surface at the rate
$\unicode[STIX]{x1D714}_{0}$
and also use the fact that, in gyrokinetics, gradients parallel to
$\boldsymbol{b}$
are always small compared to those perpendicular to it, we find
$$\begin{eqnarray}\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\approx x\frac{B_{p}^{2}R^{2}}{B}\frac{\text{d}\unicode[STIX]{x1D714}}{\text{d}\unicode[STIX]{x1D713}}\boldsymbol{e}_{y}\boldsymbol{\cdot }\unicode[STIX]{x1D735}=\left(\frac{qRB_{p}}{rB}|\unicode[STIX]{x1D735}r|\right)x\unicode[STIX]{x1D6FE}_{E}\frac{v_{thi}}{a}\boldsymbol{e}_{y}\boldsymbol{\cdot }\unicode[STIX]{x1D735},\end{eqnarray}$$
with
$\unicode[STIX]{x1D6FE}_{E}$
as defined in (A 7). The prefactor enclosed in the parentheses is close to unity and so
$\unicode[STIX]{x1D6FE}_{E}$
is the non-dimensionalised flow shear that operates on the distribution function. There is also free-energy injection associated with the presence of the flow shear, as is manifest in the presence of a term proportional to
$\unicode[STIX]{x1D6FE}_{E}$
on the right-hand side of (A 3), but, at the values of
$\unicode[STIX]{x1D6FE}_{E}$
considered here, the (destabilising) effect of this term, while included in our simulations, is irrelevant in comparison with that of the ion-temperature gradient (see Schekochihin et al.
Reference Schekochihin, Highcock and Cowley2012 for further details on this).
GS2 simulation parameters obtained from diagnostic measurements of the MAST discharge no. 27268 and appropriately normalised. Here
$\unicode[STIX]{x1D713}_{tor}=\left(1/2\unicode[STIX]{x03C0}\right)^{2}\int _{0}^{V}\,\text{d}V\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}$
is the toroidal magnetic flux and
$\unicode[STIX]{x1D713}_{pol}=\left(1/2\unicode[STIX]{x03C0}\right)^{2}\int _{0}^{V}\,\text{d}V\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}$
is the poloidal magnetic flux. See http://gyrokinetics.sourceforge.net for instructions on how the code is run with these parameters.
GS2 solves the GK equations (A 3) and (A 6) in a flux tube following the flux surface once around the torus poloidally (see figure 1
a and discussion in the main text). Figure 1(b) shows the poloidal projection of the MAST flux surfaces with the flux surface at
$r=0.8$
marked. The marked flux surface is the location at which we solve the GK equation. It can be described as located in the outer core of the device.
Each GS2 simulation requires input of a number of constant parameters that define the magnetic-field geometry (e.g. elongation of the flux surface, its triangularity, magnetic shear, etc.), the properties of the mix of the participating particle species (their masses, charges, densities, temperatures, collisionalities, etc.) and the local thermal equilibrium properties of the plasma – in particular, the gradients defined in (A 7). The use of the local formulation of gyrokinetics requires that
$\unicode[STIX]{x1D70C}^{\ast }\equiv \unicode[STIX]{x1D70C}_{i}/a\ll 1$
. At the radial location considered in this paper,
$\unicode[STIX]{x1D70C}_{i}\approx 6\times 10^{-3}$
m and the minor radius
$a\approx 0.6$
m, which implies
$\unicode[STIX]{x1D70C}^{\ast }\sim 10^{-2}$
and justifies the local approximation of gyrokinetics that we have used. table 1 gives a list of the equilibrium parameters, which have been determined via diagnostic measurements of the MAST discharge no. 27268 (at
$t=0.25$
s). In particular,
$T_{i}$
and
$\unicode[STIX]{x1D714}$
were obtained from charge-exchange-recombination spectroscopy measurements (Conway et al.
Reference Conway, Carolan, McCone, Walsh and Wisse2006) and
$T_{e}$
and
$n_{e}$
were obtained from a Thomson scattering diagnostic (Scannell et al.
Reference Scannell, Walsh, Dunstan, Figueiredo, Naylor, O’Gorman, Shibaev, Gibson and Wilson2010). The magnetic geometry in our simulations is described by the Miller et al. (Reference Miller, Chu, Greene, Lin-Liu and Waltz1998) parametrisation and the geometric parameters were obtained from a reconstruction of the equilibrium using the EFIT code (Lao et al.
Reference Lao, St. John, Stambaugh, Kellman and Pfeiffer1985). All these parameters were fixed at the same values in all our simulations, except the ion-temperature gradient
$\unicode[STIX]{x1D705}_{T}$
and the flow shear
$\unicode[STIX]{x1D6FE}_{E}$
(their experimental values are
$\unicode[STIX]{x1D705}_{T}=5.1$
and
$\unicode[STIX]{x1D6FE}_{E}=0.16$
; see figure 2 for error bars on these values).
The resolution of our simulations (with corresponding GS2 input parameters) was as follows: 128 radial modes (nx), 96 binormal modes (ny), 20 parallel grid points (ntheta), 16 energy grid points (negrid) and 27 pitch-angle grid points (
$\mathtt{ngauss}=8$
). The box sizes were approximately
$200\unicode[STIX]{x1D70C}_{i}$
(
$\mathtt{x0}=10$
and
$\mathtt{jtwist}=80$
) and
$62\unicode[STIX]{x1D70C}_{i}$
(
$\mathtt{y0}=10$
) in the radial and binormal directions, respectively, and
$2\unicode[STIX]{x03C0}$
in the parallel direction, given that GS2 uses the poloidal angle as the parallel coordinate.
