1. Introduction
Turbulence is a common feature of magnetised plasmas, appearing in systems as varied as the solar wind, astrophysical accretion disks and laboratory plasmas. According to the most common paradigm for such systems, a unique stationary turbulent state can be identified given a certain stirring mechanism and a set of equilibrium plasma parameters. Multistable solutions – for which identical parameters admit distinct turbulent states – are known to occur in neutral fluids (Snedeker & Donaldson Reference Snedeker and Donaldson1966; Burggraf & Foster Reference Burggraf and Foster1977; Schmucker & Gersten Reference Schmucker and Gersten1988), where they are associated with bifurcations and hysteretic behaviour (Shtern & Hussain Reference Shtern and Hussain1993; Ravelet et al. Reference Ravelet, Marié, Chiffaudel and Daviaud2004). Multistability has also been reported in weakly magnetised systems of charged fluids (Simitev & Busse Reference Simitev and Busse2009; Latter & Papaloizou Reference Latter and Papaloizou2012). In this work, we report the discovery of bistable turbulence in a strongly magnetised plasma, using direct numerical simulations. We find that bistability arises in such a plasma through the interplay of two crucial mechanisms: an externally imposed mean flow shear and self-generated ‘zonal’ flows. Our observations are made in a toroidal geometry typically encountered in magnetic-confinement-fusion experiments, although the results may be generalisable to other systems.
Previous studies have already established that sheared flows play an important role in regulating turbulence. In the absence of an externally imposed mean flow shear, the plasma is known to generate sheared ‘zonal’ flows spontaneously, contributing to the saturation of turbulence (Biglari, Diamond & Terry Reference Biglari, Diamond and Terry1990; Dimits et al. Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000; Rogers, Dorland & Kotschenreuther Reference Rogers, Dorland and Kotschenreuther2000; Diamond et al. Reference Diamond, Itoh, Itoh and Hahm2005; Colyer et al. Reference Colyer, Schekochihin, Parra, Roach, Barnes, Ghim and Dorland2017; van Wyk et al. Reference van Wyk, Highcock, Field, Roach, Schekochihin, Parra and Dorland2017; Ivanov et al. Reference Ivanov, Schekochihin, Dorland, Field and Parra2020). So far, it is therefore usually assumed that the effect of zonal flows is to suppress turbulence. When a mean flow shear is imposed, it provides an additional mechanism for suppressing turbulence. Specifically, the shear in the mean flow perpendicular to the magnetic field has been found to reduce turbulent fluctuations (Waelbroeck & Chen Reference Waelbroeck and Chen1991; Artun & Tang Reference Artun and Tang1992; Dimits et al. Reference Dimits, Williams, Byers and Cohen1996; Synakowski et al. Reference Synakowski, Batha, Beer, Bell, Bell, Budny, Bush, Efthimion, Hammett and Hahm1997; Casson et al. Reference Casson, Peeters, Camenen, Hornsby, Snodin, Strintzi and Szepesi2009; Mantica et al. Reference Mantica, Strintzi, Tala, Giroud, Johnson, Leggate, Lerche, Loarer, Peeters and Salmi2009; Highcock et al. Reference Highcock, Barnes, Schekochihin, Parra, Roach and Cowley2010; Barnes et al. Reference Barnes, Parra, Highcock, Schekochihin, Cowley and Roach2011; Highcock et al. Reference Highcock, Schekochihin, Cowley, Barnes, Parra, Roach and Dorland2012; Schekochihin, Highcock & Cowley Reference Schekochihin, Highcock and Cowley2012; van Wyk et al. Reference van Wyk, Highcock, Schekochihin, Roach, Field and Dorland2016, Reference van Wyk, Highcock, Field, Roach, Schekochihin, Parra and Dorland2017). It has also been shown that the effect of the mean shear weakens away from marginal stability (Fox et al. Reference Fox, van Wyk, Field, Ghim, Parra and Schekochihin2017; van Wyk et al. Reference van Wyk, Highcock, Field, Roach, Schekochihin, Parra and Dorland2017).
Here, we find that the transition from turbulent states where the mean shear plays an essential role to states where it appears to matter only marginally is characterised by a discontinuous jump in the level of turbulent transport. Most importantly, we show that, in a region of parameter space near this transition, two distinct turbulent states exist with identical equilibrium parameters but dramatically different levels of transport. We find that the presence of strong zonal flows is a feature of higher-transport states – the opposite of what is usually assumed. The main result is presented in figure 1.
Figure 1. Dependence of the turbulent ion heat flux on the inverse ion-temperature-gradient scale length. In the simulations labelled by green ‘$+$
’ signs, the externally imposed mean flow shear was set to zero. For all other simulations, $\gamma _E=-0.079$
. Zonal modes are artificially zeroed out in simulations labelled by black crosses. The black circle denotes a simulation where amplitudes decay with time and no saturated turbulent state is observed. The dashed line marks the temperature gradient below which there is no effective linear instability ($\langle \gamma \rangle _t < 0$
) in the presence of mean flow shear.
This discovery has important implications for research in nuclear fusion. In experiments, turbulent fluxes are set by the external injection of particles, heat and momentum into the plasma, with profile gradients evolving until a stationary state is reached. However, because of computational cost, direct numerical simulations only consider a small fraction of the device's volume, in which they solve the inverse problem: for given local equilibrium quantities (e.g. profile gradients), the simulations determine the associated turbulent fluxes. The flux-to-gradient problem and its inverse can be considered equivalent if a one-to-one correspondence exists between the turbulent transport and the equilibrium parameters. Our work shows that this correspondence is not always one-to-one, which poses a challenge for modelling transport – and thus for designing future fusion devices. Finally, bistability has some remarkable consequences, such as the possibility for bifurcations of turbulent transport and gradient-relaxation cycles to develop (these have previously been considered in the absence of mean flow shear Peeters et al. Reference Peeters, Rath, Buchholz, Camenen, Candy, Casson, Grosshauser, Hornsby, Strintzi and Weikl2016).
2. Modelling plasmas with a sheared mean flow
We consider equilibrium parameters obtained from a fusion experiment conducted at the Joint European Torus facility (discharge no. 68448 Siren et al. Reference Siren, Varje, Weisen and Giacomelli2019). The plasma is confined by magnetic fields that trace out nested toroidal surfaces, with the equilibrium density and temperature staying constant along the field lines. External heating sources sustain an ion temperature gradient between the hotter core and the colder edge of the plasma, which then drives the dominant linear instability (Romanelli Reference Romanelli1989; Cowley, Kulsrud & Sudan Reference Cowley, Kulsrud and Sudan1991). A sheared mean toroidal flow is generated by injecting beams of neutral hydrogenic atoms into the plasma. The ratio of thermal to magnetic pressure is small, so the turbulent fluctuations can be assumed electrostatic. We also neglect any trace impurities in the plasma and only consider two kinetic species – the electrons and the main deuterium ions. The simulations include collisions, as well as a small amount of numerical hyperviscosity (Belli Reference Belli2006). The numerical parameters used for this work are provided as supplementary material available at https://doi.org/10.1017/S0022377822000691 alongside the present publication.
The model used for this work is presented in Appendix A. The time evolution of turbulent fluctuations is described by following the approach of local $\delta f$
-gyrokinetics (Catto Reference Catto1978; Frieman & Chen Reference Frieman and Chen1982; Sugama & Horton Reference Sugama and Horton1998; Abel et al. Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013). This approximation relies on the scale separations present in the system by defining an asymptotic-expansion parameter $\rho _{*s} = \rho _s/a \ll 1$
for a species $s$
, where $\rho _s$
is the particle's gyroradius around a magnetic field line and $a$
is the minor radius of the device. As a result, the rapid gyromotion of particles can be averaged out. In this approach, the kinetic equation and the quasineutrality condition form a closed system of equations for the fluctuating probability distribution function of charged rings and for the electrostatic potential $\varphi$
. The system is solved numerically with the code GS2 (Kotschenreuther, Rewoldt & Tang Reference Kotschenreuther, Rewoldt and Tang1995; Barnes et al. Reference Barnes, Abel, Dorland, Ernst, Hammett, Ricci, Rogers, Schekochihin and Tatsuno2009; Highcock Reference Highcock2012) in a filament-like simulation domain (Beer, Cowley & Hammett Reference Beer, Cowley and Hammett1995) that follows a magnetic field line as it wraps around the torus. The code then computes the turbulent contributions to the heat and momentum fluxes exiting the core of the plasma, which we denote by $Q_s$
and $\varPi _s$
, respectively. In the following figures, we normalise $Q_s$
to the so-called gyro-Bohm value $Q_{\mathrm {gB}} = \langle \left \lvert \boldsymbol {\nabla }\psi \right \rvert \rangle _{\psi }n_iT_iv_{\mathrm {th},i}\rho _{*s}^2$
, where $\psi$
is the poloidal magnetic flux, $\langle \cdot \rangle _{\psi }$
the average over a magnetic-flux surface, $T_i$
the ion temperature multiplied by the Boltzmann constant $k_{B}$
, $v_{\mathrm {th},i}=\sqrt {2 T_i/m_i}$
the ion thermal speed and $m_i$
the ion mass. The ion temperature gradient is specified through $a/L_{T_i}$
, where $L_{T_i}$
is the local $e$
-folding length scale of $T_i$
. Finally, we denote by $\gamma _E$
the rate at which the mean flow is sheared across magnetic surfaces, normalised by $v_{\mathrm {th},i}/a$
.
When $\gamma _E\neq 0$
, linear modes (known as Floquet modes) are advected along the magnetic field lines, passing through regions of the plasma that are alternately stable (inboard of the torus) and unstable (outboard) to the ion temperature gradient (Waelbroeck & Chen Reference Waelbroeck and Chen1991). As a consequence, their linear growth rate is time dependent with a Floquet period $T_{{F}} = 2{\rm \pi} \hat {s}/\gamma _E$
. Here, $\hat {s}$
(defined in Appendix A) measures how the twisting of magnetic field lines around the torus changes with the minor radius. In the following, we denote the time-averaged growth rate by $\langle \gamma \rangle _t$
, and the maximum instantaneous growth rate by $\gamma _{\mathrm {max}}$
. In this work, the code GS2 was used with an improved algorithm for background flow shear described in McMillan, Ball & Brunner (Reference McMillan, Ball and Brunner2019) and Christen, Barnes & Parra (Reference Christen, Barnes and Parra2021), although it has been verified that the same conclusions are reached when using the original algorithm devised by Hammett et al. (Reference Hammett, Dorland, Loureiro and Tatsuno2006).
3. Two distinct turbulent states
Near marginal stability, we find that two distinct turbulent states can be obtained at identical equilibrium parameters. This is shown in figure 1, where saturated values of $Q_i/Q_{\mathrm {gB}}$
are plotted against $a/L_{T_i}$
for a particular value of $\gamma _E$
. We find that the fluxes computed in the low-transport state match the levels of transport observed in the experiment, while the fluxes computed in the high-transport state differ from it by an order of magnitude. For equilibrium parameters where the two states exist, it is the initial size of the fluctuation amplitudes that determines which state is observed in a simulation. While the impact of initial conditions on gyrokinetic simulations was explored in previous work, such as Pueschel, Kammerer & Jenko (Reference Pueschel, Kammerer and Jenko2008), our work is the first to obtain two distinct, saturated and finite-amplitude turbulent states with identical equilibrium parameters. Interestingly, we note that both low-transport and high-transport states can exist above and below the threshold for linear instability. Previous work had already established this for a single, finite-amplitude turbulent state sustained either by a linear instability (known as supercritical turbulence) or by transient linear growth (known as subcritical turbulence Highcock et al. Reference Highcock, Barnes, Schekochihin, Parra, Roach and Cowley2010; Barnes et al. Reference Barnes, Parra, Highcock, Schekochihin, Cowley and Roach2011; Schekochihin et al. Reference Schekochihin, Highcock and Cowley2012).
The two states observed here are distinguished by significant differences in the amplitudes of their turbulent fluctuations and by the spatial structure of the turbulence. In figure 2, we show a typical snapshot of turbulence in a low-transport state. The contours of the fluctuating electrostatic potential are plotted in the plane perpendicular to the magnetic field at the outboard of the torus. The $x$
coordinate measures the distance along the normal to a magnetic-flux surface and the $y$
coordinate labels the magnetic field lines within the surface. The simulation is done in the frame moving with the mean flow at $x=0$
: the $y$
-component of the mean flow thus has the opposite sign to $x$
. The turbulent eddies feature a clear tilt as they are being sheared by the mean flow, similarly to Shafer et al. (Reference Shafer, Fonck, McKee, Holland, White and Schlossberg2012), van Wyk et al. (Reference van Wyk, Highcock, Schekochihin, Roach, Field and Dorland2016, Reference van Wyk, Highcock, Field, Roach, Schekochihin, Parra and Dorland2017) and Fox et al. (Reference Fox, van Wyk, Field, Ghim, Parra and Schekochihin2017). In figure 3, we show consecutive snapshots of turbulence with the same equilibrium parameters as in figure 2, but in the high-transport state: bands of high-amplitude eddies propagate radially across the simulation domain, and eddies do not feature any clear tilt. This intermittent high-transport state is reminiscent of advecting structures reported in McMillan et al. (Reference McMillan, Jolliet, Tran, Villard, Bottino and Angelino2009), McMillan, Pringle & Teaca (Reference McMillan, Pringle and Teaca2018) and Chandrarajan Jayalekshmi (Reference Chandrarajan Jayalekshmi2020).
Figure 2. Consecutive snapshots of the turbulence in real space for the low-transport state where $a/L_{T_i} = 1.76$
and $\gamma _E=-0.079$
. In the top panels, the fluctuating electrostatic potential is plotted at three successive times at the outboard of the torus, in the plane perpendicular to $\boldsymbol {B}$
. In the bottom panels, the zonal flow is plotted at the same times.
Figure 3. Same as figure 2 but for the high-transport state. In the bottom panels, the zonal shear averaged between the two vertical dashed lines is compared with the externally imposed mean flow shear.
4. The role of zonal modes
We find that the presence of ‘zonal’ modes is a crucial distinguishing feature between the two states. Modes are called zonal when they have no spatial variation other than in the radial ($x$
) direction. They are linearly stable and cannot be sheared by a toroidal mean flow, but they can exchange energy with non-zonal modes via nonlinear interactions. Zonal modes include zonal flows with a shearing rate $\gamma _{{Z}}$
, which can affect the rest of the turbulence in a manner analogous to the mean flow shear $\gamma _E$
. The zonal flows are known to develop through a secondary instability of the modes driven unstable by the temperature gradient (Rogers et al. Reference Rogers, Dorland and Kotschenreuther2000; Diamond et al. Reference Diamond, Itoh, Itoh and Hahm2005; Ivanov et al. Reference Ivanov, Schekochihin, Dorland, Field and Parra2020).
When the amplitudes of zonal modes become large enough, we find that the zonal shear can compete with, and indeed obviate, the mean flow shear. Such a negation of the equilibrium shear by a zonal shear was already explored in previous work, e.g. McMillan et al. (Reference McMillan, Jolliet, Tran, Villard, Bottino and Angelino2009, Reference McMillan, Pringle and Teaca2018). In the lower panels of figure 3, we plot the zonal flow $V_{{Z}}$
($\propto \partial \varphi _{{Z}} / \partial x$
where $\varphi _{{Z}}$
is the zonal part of the potential). We observe radially propagating bands within which the zonal shear $\gamma _{{Z}}$
($\propto \partial V _{{Z}} / \partial x$
) is of the same order of magnitude as the background shear $\gamma _E$
, but carries the opposite sign. In such bands, where the zonal and mean shears oppose each other, non-zonal fluctuations grow faster and feed the zonal modes nonlinearly, until the system settles in the high-transport state. We also show in figure 1 that the transport obtained in the complete absence of mean flow shear (green ‘$+$
’ symbols) is much closer to the high-transport states than to the low-transport states.
In low-transport states, zonal modes do not seem to play a crucial role for the turbulent dynamics: unlike in the high-transport states, no long-lived structures with $\gamma _{{Z}}$
opposing $\gamma _E$
are observed. Further evidence of the weaker impact of zonal modes in low-transport states can be seen in simulations where we artificially set zonal modes to zero at every time step, indicated by black crosses in figure 1. Despite this unphysical truncation introduced in the system, and independently of the initial condition, a saturated state is obtained that is – apart from a slight change in the flux – indistinguishable from the low-transport state.
The occurrences of low-transport and high-transport states for a range of mean flow shear rates are shown in figure 4, where we plot the ratio $\textrm {rms}[\gamma _{{Z}}]/\gamma _E$
. Here, we define the root mean square of the zonal shear as $\textrm {rms}[\gamma _{{Z}}] = \sqrt {\langle \varphi _{{Z}}^2\rangle _{t,x}}/\ell _{x,Z}^2$
and we denote by $\ell _{x,Z}$
the radial correlation length of the zonal modes, which we define in Appendix B. As a result of the interplay between the zonal modes and the mean flow, high-transport states are only obtained when the initial fluctuation amplitudes are sufficiently large, or when the fluctuation amplitudes become large enough for the zonal shear to start competing with $\gamma _E$
(e.g. when $a/L_{T_i}$
is increased or $\gamma _E$
is decreased past a certain threshold). Figure 4(a) confirms that the magnitude of the zonal shear in the high-transport states is comparable to that of the mean shear and panel (b) indicates that low-transport states only survive when the zonal shear is much smaller than $\gamma _E$
(roughly by an order of magnitude). This last result could be due to a feedback mechanism whereby even a weak $\gamma _{{Z}}$
can partially oppose $\gamma _E$
, allowing for larger turbulent amplitudes – and therefore larger zonal modes – to develop.
Figure 4. Root-mean-square zonal shear versus $\gamma _E$
for the two turbulent states. High-transport states are shown in (a), and low-transport states in (b). Black-bordered markers indicate parameters at which either a high-transport or a low-transport state can be obtained, depending on the initial size of fluctuation amplitudes. The parameters of the experiment considered here are shown by a black star in (b). The dashed line marks the temperature gradient as a function of $\gamma _E$
below which turbulence is subcritical ($\langle \gamma \rangle _t < 0$
).
5. Two correlation time scales
In a saturated turbulent state, the correlation time of the turbulence can be estimated from the gyrokinetic equation by $\tau = [c \left [ \varphi _{\mathrm {NZ}} \right ]_{\mathrm {rms}} / (B \ell _x \ell _y)]^{-1}$
, where $\left [ \varphi _{\mathrm {NZ}} \right ]_{\mathrm {rms}}$
is the root-mean-square value of the non-zonal electrostatic potential, $\ell _y$
is the eddy correlation length in the $y$
direction (defined in Appendix B), $c$
is the speed of light and $B$
is the magnetic field strength. For high-transport states, figure 5 shows that $\tau \sim 1/\gamma _{\mathrm {max}}$
, where $\gamma _{\mathrm {max}}$
is the maximum instantaneous linear growth rate in the presence of flow shear (close to the linear growth rate in the absence of flow shear). For low-transport states, $1/\gamma _{\mathrm {max}} < \tau < 1/\langle \gamma \rangle _t$
, which may suggest that the average linear growth rate plays a role in setting the saturated turbulent amplitudes in those states. In order for $\langle \gamma \rangle _t$
to be relevant in a turbulent state, eddies must be able to survive longer than a Floquet period, i.e. $\tau \gtrsim T_{{F}}$
, as is approximately the case for the low-transport states in figure 5. While the competition between zonal and mean flow shear is likely a generic feature of magnetised plasma turbulence, the existence of $\langle \gamma \rangle _t \neq \gamma _{\mathrm {max}}$
requires toroidicity. Further studies are needed to determine if these two distinct growth rates are a necessary feature of the bistable turbulence reported here – and thus if similar bistable states are likely to be found beyond toroidal plasmas.
Figure 5. Correlation time of turbulent eddies in the low-transport state, high-transport state and in the absence of a mean flow shear. The vertical dashed line marks the temperature gradient below which there is no effective instability in the presence of a mean flow shear, i.e. $\langle \gamma \rangle _t \leq 0$
. In the simulations labelled by green ‘$+$
’ signs, the externally imposed mean flow shear was set to zero. For all other simulations, $\gamma _E=-0.079$
was used. Considering other $\gamma _E$
values leads to similar results.
6. Consequences of bistability
The bistability reported in this work may lead to the existence of bifurcations. As we have argued, low-transport states cease to exist when the fluctuation amplitudes increase past a certain threshold value. If we now consider a plasma in which, instead of being fixed, the temperature gradient is slowly increasing in time, a discontinuous jump will be triggered from a low-transport state to the high-transport branch. The same jump can be triggered by decreasing the mean flow shear in a low-transport state. As shown in figure 4, we observe that the subcritical low-transport states exist closer to marginal stability in the $(\gamma _E,a/L_{T_i})$
plane than the subcritical high-transport states. We attribute this to the intermittent nature of the high-transport state associated with the radially propagating bands shown in figure 3. Previous work with neutral fluid flows (Faisst & Eckhardt Reference Faisst and Eckhardt2004), accretion disks (Rempel, Lesur & Proctor Reference Rempel, Lesur and Proctor2010) and fusion plasmas (Barnes et al. Reference Barnes, Parra, Highcock, Schekochihin, Cowley and Roach2011; Highcock et al. Reference Highcock, Barnes, Parra, Schekochihin, Roach and Cowley2011; van Wyk et al. Reference van Wyk, Highcock, Field, Roach, Schekochihin, Parra and Dorland2017) has indeed shown that the survival of subcritical turbulence over long times is compromised by rare, large fluctuations. Similarly to the transition from low to high transport, we argue that subcritical high-transport states can drop to the low-transport branch if the temperature gradient slowly decreases in time. From figure 4, we expect that the same transition could be achieved by increasing the mean flow shear in a subcritical high-transport state.
Existence of bifurcations opens up the possibility for relaxation cycles of the mean gradients to develop. This hinges on two findings that we show in figure 6. First, we observe a significant gap between the highest heat flux obtained in low-transport states and the lowest heat flux obtained in high-transport states. Second, we observe that the ratio of the turbulent momentum flux to heat flux, $\varPi _i/Q_i$
, is almost identical in the low-transport and high-transport states.
Figure 6. Dependence of the ion heat flux (a,b) and the momentum-to-heat-flux ratio (c,d) on the imposed flow shear and the inverse ion-temperature-gradient scale length. The top panels show results for the high-transport states, the bottom panels for the low-transport states. Dotted areas in the upper (respectively lower) panels indicate areas where no high-transport (respectively low-transport) state could be obtained. The grey areas indicate parameter ranges where no simulations were run. There is a gap between the values of the heat flux obtained in (a) and those obtained in (b). The path defined by points A, B, C and D gives an example of the successive stages of a gradient-relaxation cycle, when the heat injected into the plasma corresponds to a flux within the aforementioned gap.
We now consider a thought experiment in which an external power $P$
is injected into the volume bounded by a given magnetic-flux surface of area $S$
, via a beam of neutral atoms with energy $E$
. As is argued in Parra et al. (Reference Parra, Barnes, Highcock, Schekochihin and Cowley2011), the turbulent heat flux exiting the magnetic-flux surface is $Q_i \sim P/S$
, and $\varPi _i/Q_i \sim E^{-1/2} v_{\mathrm {th},i}/a$
. Thus, a given beam configuration corresponds to a unique pair $(Q_i, \varPi _i/Q_i)$
. We consider an initial situation where the input power is such that $(P,E)$
corresponds to the levels of turbulent fluxes of a low-transport state (point A in figure 6). From this initial stationary state, we increase $P$
by small successive increments, keeping $E$
fixed. In response, the plasma equilibrium will evolve through a succession of low-transport stationary states with ever larger $Q_i$
, but with $\varPi _i/Q_i$
staying constant (along the solid arrow up to point B in figure 6).
Above a certain threshold, we find that there is a range of values of $P$
with no corresponding solutions for the turbulent fluxes: in figure 6, these are the powers too high to match the low-transport state at point B, but too low to match the high-transport state at point C. It is then interesting to ask what will happen to a plasma where the input power falls into this gap. One (unexciting) possibility is that an actual solution may exist outside of the region of parameter space explored here, and that the plasma migrates to that solution. Another (more interesting) possibility would be for the temperature gradient to continue increasing until the plasma transitions to a high-transport state (jumping from B to C in figure 6). In this state, the outgoing heat flux is larger than what can be sustained by the external power input, and the temperature gradient starts to flatten. As $a/L_{T_i}$
decreases (from C to D), the turbulent fluctuations remain too large to allow a transition back to a low-transport state. Eventually, $a/L_{T_i}$
becomes too small for the high-transport state to survive, and the system transitions back to the lower state (from D to A). The flux is now too low compared with the power input, so the temperature gradient builds up again, and the cycle repeats itself. In this scenario, no proper steady state is reached when the input power falls within a ‘forbidden’ gap, and the temperature gradient and mean flow shear would experience periodic relaxation cycles.
7. Discussion
We have found that near-marginal turbulence in fusion devices is bistable, and regulated by the competition between external shear and zonal modes. The existence of bistability suggests a new approach to long-standing questions around bifurcations and gradient-relaxation cycles observed in fusion devices (von Goeler, Stodiek & Sauthoff Reference von Goeler, Stodiek and Sauthoff1974; Wagner et al. Reference Wagner, Becker, Behringer, Campbell, Eberhagen, Engelhardt, Fussmann, Gehre, Gernhardt and Gierke1982; Hastie Reference Hastie1997; Connor Reference Connor1998). This work also presents a new challenge for a research area where the prevailing assumption has been a one-to-one correspondence between plasma parameters and turbulent transport. Further work could focus on how the extent of the bistable region might be modified, for example by exploring the effect of collisions on the saturation of zonal modes (Colyer et al. Reference Colyer, Schekochihin, Parra, Roach, Barnes, Ghim and Dorland2017; Weikl et al. Reference Weikl, Peeters, Rath, Grosshauser, Buchholz, Hornsby, Seiferling and Strintzi2017). Another avenue of interest may be to determine how bistability manifests itself in flux-driven gyrokinetic simulations. Experiments could test the existence of bistability in fusion devices, following a scenario similar to the one described in figure 6. The idea of relaxation cycles discussed in § 6 could be considered in the context of subcritical turbulence, where transitions might occur from a situation with no turbulent transport to a state with a finite level of turbulent transport. Lastly, we note that the details of the plasma analysed here, such as the exact nature of the drive for turbulence and perhaps toroidicity, do not appear to be crucial to our understanding of bistable states: the only requirements we have identified so far are an applied flow shear and the ability of the plasma to generate zonal flows. Therefore, we expect that similar effects may be observed in a variety of systems.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/S0022377822000691.
Acknowledgements
The authors are grateful to H. Weisen and P. Sirén for providing the experimental data used in this work. They especially thank F. Parra for his insightful comments on the manuscript. They would also like to thank O. Beeke, J. Parisi and J. Ruiz Ruiz for very fruitful discussions. N.C. was supported by a Berrow Foundation Scholarship, the Steppes Fund for Change and the Fondation Hélène et Victor Barbour. The work of M.H. was funded by EPSRC grant EP/R034737/1, as was, in part, the work of M.B. and A.A.S. Computing resources were provided on the ARCHER High Performance Computer through the Plasma HEC Consortium, EPSRC grant EP/L000237/1 under project e607, on the EUROfusion High Performance Computer (Marconi-Fusion) under the projects FUA34_MULTEI and FUA35_OXGK and on the JFRS-1 supercomputer system at the International Fusion Energy Research Centre's Computational Simulation Centre (IFERC-CSC) at the Rokkasho Fusion Institute of QST (Aomori, Japan) under the project MULTEI.
Editor Paolo Ricci thanks the referees for their advice in evaluating this article.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Gyrokinetic system
In this work, we follow the $\delta f$
gyrokinetic approach (Catto Reference Catto1978; Frieman & Chen Reference Frieman and Chen1982; Sugama & Horton Reference Sugama and Horton1998; Abel et al. Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013), which relies on the scale separations present in the plasma to describe the time evolution of turbulent fluctuations. The ratio of gyroradius to machine size $\rho _{*s} = \rho _s/a \ll 1$
for species $s$
is used as the asymptotic-expansion parameter. The minor radius of the device is denoted by $a$
and the gyroradius is given by $\rho _s = \vert \hat {\boldsymbol {b}}\times \boldsymbol {v}/\varOmega _s\vert$
, where $\hat {\boldsymbol {b}}$
is the unit vector in the direction of the magnetic field $\boldsymbol {B}$
and $\boldsymbol {v}$
is the velocity of the particle. The gyrofrequency of the particle is $\varOmega _s=eZ_s B/m_s c$
, where $Z_s$
and $m_s$
are, respectively, the charge number and mass of the particle, $e$
is the elementary charge, and $c$
is the speed of light. The amplitudes of the fluctuations are ordered to be $O(\rho _{*s})$
smaller than the corresponding mean quantities. The turbulent time scale is ordered to be $O(\rho _{*s}^2)$
shorter than the time scale of the evolution of mean plasma parameters, but $O(\rho _{*s}^{-1})$
longer than the Larmor periods of the particles. Moreover, it is assumed that fluctuations can stretch far along magnetic field lines, but that they only span a few gyroradii across field lines. The orderings in time and space can be summarised as
where $\delta f_s$
is the fluctuating part of the distribution function of particles and $F_s$
is their mean distribution function (averaged over the turbulent time scales and over the turbulent length scales across $\boldsymbol {B}$
). Here, ${\rm d}/{\rm d} t = \partial /\partial t + \boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla }$
is the convective time derivative with respect to the mean flow $\boldsymbol {u}$
.
The geometry of the system considered here is typical of magnetic-confinement-fusion experiments. The plasma is confined in a toroidally shaped magnetic cage. The magnetic field lines of this cage wind around the torus, tracing out nested toroidal surfaces, known as magnetic-flux surfaces. The rapid gyromotion about magnetic field lines limits the ability of charged particles to move across these flux surfaces.
We focus on plasmas with a mean flow such that $\rho _{*s} \ll \left \lvert \boldsymbol {u}\right \rvert\!/ v_{\mathrm {th},i} \ll 1$
, where $v_{\mathrm {th},i} = \sqrt {2T_i/m_i}$
is the ion thermal speed and $T_i$
is the ion temperature multiplied by the Boltzmann constant $k_{B}$
. In this ‘intermediate-flow’ ordering, we can neglect the centrifugal force, and the mean flow is purely toroidal (Catto, Bernstein & Tessarotto Reference Catto, Bernstein and Tessarotto1987; Abel et al. Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013): $\boldsymbol {u}=\varOmega _{\phi } R^2 \boldsymbol {\nabla }\phi$
with $\varOmega _{\phi }$
the angular rotation frequency, $R$
the major radius of the torus and $\phi$
the toroidal angle. It follows that the perpendicular and parallel flow shear rates are related by a geometrical factor, and both can be expressed in terms of the shearing rate $\gamma _E = (r_{\psi,0}/q_0) \partial \varOmega _{\phi }/\partial r _{\psi } \rvert _{r_{\psi,0}}$
, where the subscript ‘$0$
’ denotes quantities evaluated on the flux surface of interest, $r_{\psi }$
is the half-width of the flux surface at the height of the magnetic axis and the safety factor $q=(2{\rm \pi} )^{-1}\int _0^{2{\rm \pi} }{\rm d}\theta (\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }\phi )/(\boldsymbol {B}\boldsymbol {\cdot } \boldsymbol {\nabla }\theta )\vert _\psi$
is the number of toroidal turns required by a field line to wind once around the torus poloidally. The magnetic shear appearing in § 2 is defined as $\hat {s} = (r_{\psi,0}/q_0) \partial q/\partial r _{\psi } \vert _{r_{\psi,0}}$
. We further restrict consideration to cases with low thermal-to-magnetic-pressure ratio (plasma beta) and hence only retain electrostatic fluctuations. We neglect all effects associated with impurities in the plasma, and only consider electrons and the main hydrogenic ion species.
After averaging over the rapid gyromotion of particles, the gyrokinetic equation can be written as
in $(\boldsymbol {R}_s, \varepsilon _{s}, \mu _{s}, \vartheta )$
coordinates, where $\boldsymbol {R}_s = \boldsymbol {r} - \boldsymbol {\rho }_s$
is the particle's gyrocentre, $\boldsymbol {r}$
is its position, $\varepsilon _{s}=m_s w^2/2$
is its kinetic energy, $\mu _{s}=m_s w_\perp ^2/2B$
its magnetic moment and $\vartheta$
its gyrophase. Here, $\langle \cdot \rangle _{\boldsymbol {R}_s}$
denotes an average over $\vartheta$
at fixed $\boldsymbol {R}_s$
, $\varphi$
is the fluctuating electrostatic potential, $\boldsymbol {w}$
is the particle velocity relative to $\boldsymbol {u}$
, subscripts $\parallel$
and $\perp$
indicate components along and across $\boldsymbol {B}$
respectively, $F_{0,s}$
is a local Maxwellian velocity distribution, $C$
is the collision operator and $B_\phi$
is the toroidal component of $\boldsymbol {B}$
. The drift velocity due to magnetic curvature and $\boldsymbol {\nabla } B$
is $\boldsymbol {V}_{B,s} = \hat {\boldsymbol {b}}/\varOmega _s\times [ w_\perp ^2\boldsymbol {\nabla }\ln (B)/2 + w_\parallel ^2 \hat {\boldsymbol {b}}\boldsymbol {\cdot }\boldsymbol {\nabla }\hat {\boldsymbol {b}} ]$
, and the Coriolis drift velocity is $\boldsymbol {V}_{C,s} = (2w_\parallel \varOmega _{\phi }/\varOmega _s)\hat {\boldsymbol {b}}\times (\hat {\boldsymbol {z}}\times \hat {\boldsymbol {b}})$
with $\hat {\boldsymbol {z}}$
the unit vector in the vertical direction. The nonlinearity in (A3) stems from the fluctuating $\boldsymbol {E}\times \boldsymbol {B}$
drift $\boldsymbol {V}_E = c\hat {\boldsymbol {b}}/B\times \boldsymbol {\nabla }\varphi$
advecting $\delta f_s$
on the left-hand side. Perpendicular flow shear enters via the convective time derivative, while the drives from the shear in the parallel flow and the temperature gradient, respectively, enter via the $\boldsymbol {\nabla }\varOmega _{\phi }$
and $\boldsymbol {\nabla }F_{0,s}$
terms on the right-hand side. The temperature gradient is specified by the normalised inverse gradient length $a/L_{T_s}=-a\,{\rm d}(\ln T_s)/{\rm d}r_{\psi }$
. The set of equations is closed by the quasineutrality condition
where $n_s$
is the particle density and $\langle \cdot \rangle _{\boldsymbol {r}}$
denotes an average over $\vartheta$
at fixed particle position.
Fluctuations with no spatial variation other than in the radial direction are known as ‘zonal’ fluctuations, and produce sheared $\boldsymbol {E}\times \boldsymbol {B}$
drifts in the $y$
direction. The zonal flow is given by
where $\varphi _{{Z}}$
is the zonal part of the electrostatic potential. The shear of the zonal flow is $\gamma _{{Z}} = \partial V _{{Z}} / \partial r _{\psi }$
.
The system of (A3) and (A4) is solved for $\delta f_s$
and $\varphi$
using the local gyrokinetic code GS2 (Kotschenreuther et al. Reference Kotschenreuther, Rewoldt and Tang1995; Barnes et al. Reference Barnes, Abel, Dorland, Ernst, Hammett, Ricci, Rogers, Schekochihin and Tatsuno2009; Highcock Reference Highcock2012; Christen et al. Reference Christen, Barnes and Parra2021) in a filament-like simulation domain (known as a flux tube Beer et al. Reference Beer, Cowley and Hammett1995) that follows a magnetic field line around the flux surface of interest. The flux-surface label $x = (q_0/r_{\psi,0} B_r)(\psi -\psi _0)$
and the field-line label $y = (1/B_r)(\partial \psi /\partial r _{\psi })\rvert _{r_{\psi,0}}(\alpha -\alpha _0)$
are used as coordinates across $\boldsymbol {B}$
. The poloidal angle $\theta$
serves as the coordinate along $\boldsymbol {B}$
. Here, $B_r$
is a reference magnetic field strength, $\psi = \int _{0}^{r}{\rm d} r' r'R\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }\theta$
is the poloidal magnetic flux, $r$
is the minor radius of the torus and $\alpha =\phi -\int _0^\theta {\rm d}\theta ^\prime (\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }\phi )/(\boldsymbol {B} \boldsymbol {\cdot }\boldsymbol {\nabla }\theta )\vert _\psi$
labels field lines on a given flux surface. The code computes the turbulent contribution to the heat and momentum fluxes given by
respectively, with $\langle \cdot \rangle _{\psi }$
denoting the volume average over the flux tube.
Appendix B. Correlation time and correlation lengths
Given a saturated turbulent state, we estimate the eddy correlation time as being
where $\ell _x$
and $\ell _y$
denote the eddy correlation length in the $x$
and $y$
directions, respectively, and where $\left [ \varphi _{\mathrm {NZ}} \right ]_{\mathrm {rms}} = \sqrt {\langle \varphi _{\mathrm {NZ}}^2\rangle _{t,x,y}}$
is the root mean square of the non-zonal part of the electrostatic potential $\varphi _{\mathrm {NZ}}$
, averaged over time, $x$
and $y$
. The expression (B1) is obtained from the nonlinear term in the gyrokinetic equation (A3). We then define the two-point spatial correlation function
The correlation lengths $\ell _x$
and $\ell _y$
are chosen to correspond to the $e$
-folding of $\mathrm {Cor}[\varphi _{\mathrm {NZ}}]$
along the $\delta _x$
direction (adjusted to match the tilt induced by the flow shear) and the $\delta _y$
direction, respectively. Typical examples of $\mathrm {Cor}[\varphi _{\mathrm {NZ}}]$
are shown in figure 7. The zonal correlation length $\ell _{x,Z}$
corresponds to the $e$
-folding of $\mathrm {Cor}[\varphi _{{Z}}]$
along the $\delta _x$
direction. Note that the exact definition of the correlation lengths is somewhat arbitrary. Another choice, that is commonly found in the literature, is to define $\ell _x$
and $\ell _y$
as the integral of $\mathrm {Cor}[\varphi _{\mathrm {NZ}}]$
along the $\delta _x$
and $\delta _y$
axes – which, in our case, yields similar results to the $e$
-folding lengths.
Figure 7. Two-point spatial correlation function for a low-transport state (a), a high-transport state (b) and a state with no mean flow shear (c). For the three states, $a/L_{T_i}=1.76$
. In (a,b), $\gamma _E = -0.079$
.
Appendix C. Source code
The version of the GS2 code used for this work is available at https://bitbucket.org/gyrokinetics/gs2/branch/ndc_branch, with the newest commit at the time of writing being 0abdcda. The associated version of ‘Makefiles’ necessary for compilation is available at https://bitbucket.org/gyrokinetics/makefiles/branch/ndc_branch under the commit ba24979, and the additional ‘utils’ files required to run the code are available at https://bitbucket.org/gyrokinetics/utils/branch/ndc_branch under the commit 8e41f9a.
