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.