Transitioning to tokamak turbulence via states of increasing complexity

Abstract

We provide a fundamentally new perspective on subcritical turbulence in plasmas, based on coherent structures, which are obtained and characterised via direct numerical solution. The domains where these coherent states exist appear to be closely connected to the those where related turbulent states can exist, so there may be a deep connection between the stability of these coherent structures and the domain where sustained turbulence is possible. In contrast to previous descriptions of turbulence in terms of a stochastic collection of linear waves, we present a fundamentally nonlinear representation based on more general classes of translating oscillatory nonlinear solutions. In turbulent tokamak plasmas, the transport can often be completely suppressed by introducing a background shear flow, whose amplitude is an important control parameter. As this parameter is decreased below a critical value, radially localised structures appear, becoming larger and more complex, in both gyrokinetic simulations and a simpler fluid model of the plasma. For the fluid model, we directly solve for a particular class of nonlinear solutions, relative periodic orbits, and determine their stability, thus explaining why these isolated structures appear in initial-value simulations. The increase of complexity as the flow shear is reduced is explained by a series of Hopf bifurcations of these nonlinear solutions, which we quantify via stability analysis. In gyrokinetic simulations, we are able to indirectly determine the underlying relative periodic orbits by imposing symmetry conditions on the simulations.

1. Introduction

1.1. Background

Transport of heat in tokamak core plasmas is dominated by turbulent transport arising from microinstabilities. Here, we consider the ion-temperature gradient (ITG) instability (Conner & Wilson 1994; Horton 1999), which is pervasive in tokamaks (and other devices), and often responsible for significant heat transport.

Unlike in large Reynolds number fluid turbulence, which displays self-similar turbulence over several orders of magnitude in scale, observations and simulations of ITG (Candy & Waltz 2003) typically reveal a spatial dynamics dominated by a relatively narrow range of length scales (Barnes et al. 2011), suggesting an analogy with the low Reynolds number transition regime. Geometrical considerations lead to preferential excitation of zonal modes perpendicular to the gradients in the plasma, at similar scales to the ITG modes. The resulting dynamics can show long-lived propagating structures (Benkadda et al. 2001; McMillan et al. 2009; Görler et al. 2011; Jolliet & Idomura 2012; Van Wyk et al. 2016; Ivanov et al. Ivanov et al., 2020; Zhou et al. 2020), oscillations (Villard et al. 2019) and multiple attractors (Christen et al. 2022). When shear flow is present, subcritical non-steady dynamical states exist (Newton et al. 2010; Schekochihin et al. 2012), just as in classical fluid flows; these states are precursor states that occur during the transition to fully developed scale-free turbulence in fluid systems, and we follow the laboratory-plasma convention of calling them turbulence.

In the classic paradigm of tokamak turbulence, the actors are drift instability eigenmodes, which are plane waves in two dimensions (in the perpendicular directions of ballooning space), and the zonal flows, plus damped parasitic modes that serve as a free-energy sink (Bañón Navarro et al. 2011; Stoltzfus-Dueck et al. 2022; Connaughton et al. 2015). This amounts to considering linear modes as the main basis for understanding turbulence. But tokamak turbulence is often dominated by coherent and highly intermittent structures, where techniques to understand the dynamics in a wave-mode (Fourier) decomposition, such as the random-phase approximation, are invalid. We will instead examine whether a class of nonlinear solutions, RPOs (relative periodic orbits), can capture these structures and illuminate the transition to turbulence for this subcritical system. In particular, we will determine whether the domain of existence and stability of these RPOs is related to the region of parameter space where turbulence exists. Long-lived structures, which may be RPOs (or related to them), inherently balance free-energy generation and dissipation, and are thus naturally connected to the question of nonlinear saturation.

Periodic orbits are oscillations of the system state (distinct from the periodic orbits of individual particles) and play a key role in tokamak plasma physics, especially in the form of limit cycles in e.g. regular edge-localised mode production (Itoh et al. 1991), sawteeth (von Goeler et al. 1974) and chirping (Berk et al. 1997); linear oscillations, and their nonlinear extensions, are one way periodic orbits arise. Propagating features like blobs in the edge (Myra et al. 2006) are also pervasive: RPOs are a general class of periodic orbits that can also capture displacement in symmetry directions including propagation. These nonlinear structures are challenging to identify and isolate where they do not arise out of linear oscillations about a simple equilbrium, and may be unstable; there are a variety of techniques to find them, some of which will be explored in this manuscript.

This paper discusses the nature of the transition to turbulence in tokamak physics, defines these RPOs and characterises them in a simple reduced-model system (summarising aspects of Smith 2023); the gyrokinetic system is shown to exhibit analogous behaviour, with RPOs the key player in the transition regimes from simple behaviour (laminar/zero-turbulence or isolated propagating structures) to domain-filling turbulence.

Figure 1. (a) Density plot of heat flux (units of $c_s n_0 T_0 (\rho _s/R)^2$ versus radial position $x$ (in $\rho _s$ units) and time $t$ for a standard-box gyrokinetic simulation with $S=0.35$ initialised with a localised perturbation (details in § 3). (b) Density plot of turbulence intensity versus radial position $x$ and time $t$ in the plsma interchange model at $S=1.15$ (details in § 2.

1.2. The nonlinear nature of gyrokinetic (GK) turbulence

A key control parameter in tokamak turbulence is the level of background flow shear. For lower values, turbulence is domain filling but various authors (McMillan et al. 2009; Benkadda et al. 2001; Candy & Waltz 2003; Van Wyk et al. 2016) have observed localised propagating structures near the transition point of tokamak turbulence, which have been isolated by careful parameter scans (Van Wyk et al. 2016) near the transition point, and by finding the attractor in the edge of chaos (McMillan et al. 2018) (the ‘edge of chaos’ is the boundary in phase space separating the laminar and turbulent regimes). The radially localised nature of these structures indicates their nonlinear nature; in the assumed flux tube geometry (Beer et al. 1995), linear modes are space-filling plane waves in the perpendicular plane. The experimental signature of these may be ballistically propagating heat pulses found in certain devices (Zhu et al. 2015). Note also certain commonalities with the propagating structures (blobs) universally observed in the tokamak edge (Myra et al. 2006), which are also spatially isolated, with a fundamentally nonlinear propagation mechanism. These propagating structures split and evolve into space-filling turbulence at lower values of shear, in both detailed gyrokinetic simulations (figure 1 a) and in a simpler fluid plasma model (figure 1 b), suggesting the propagating features play a role as precursors from which a more complex dynamics arises (the set-up of these simulations is explained in the respective later sections).

1.3. A nonlinear description of turbulence

The presence of a background flow shear often has a strong stabilising effect on core (Waltz et al. 1998; Biglari et al. 1990; Waltz et al. 1997) and pedestal (Rogers et al. 1998) regions of tokamaks. Indeed, linear analysis of these configurations often finds absolute stability, so any infinitesimal perturbation decays to nothing over time (Waltz et al. 1998), although there can be significant transient growth of the perturbation amplitude before it undergoes damping. Possibly as a consequence of the transient growth, despite linearly stability, sufficiently large amplitude disturbances may initiate sustained turbulence (Newton et al. 2010). In this case there is a critical amplitude (dependent on the disturbance’s form) below which disturbances decay and above which turbulence ensues.

Working in configuration space, where a single experimental run is viewed as a trajectory, this critical amplitude defines the edge of chaos, a separatrix dividing initial conditions which are in the laminar state’s basin of attraction from those in that of the turbulent state. Continuing with ideas borrowed from dynamical systems, turbulence can be viewed as a strange attractor made from a tangle of exact nonlinear solutions of the underlying equations (steady states, travelling waves and (relative) periodic orbits) and their interconnected stable and unstable manifolds. The exact solutions provide a skeleton for the observed dynamics which weaves its way between them.

The dynamics of runs constrained within the edge of chaos is of lower energy but can still be chaotic and a similar description of the behaviour within it exists – of exact solutions structuring the behaviour. From this perspective we view the edge of chaos solutions as controlling the transition to turbulence, and the higher energy solutions as describing the turbulence itself. Intriguingly, the solutions embedded within these two parts of phase space are in fact connected. Typically, solutions emerge in pairs via saddle node bifurcations with one branch embedded in the edge of chaos and the other within turbulence. This drives the conclusion that it is the same nonlinear, self-sustaining processes at work within the edge of the chaos state, persistent turbulence and the isolated structures seen in tokamak simulation.

For gyrokinetic simulations of initial condition within the edge of chaos (found via bisection) the system has been seen to evolve into a travelling wave (McMillan et al. 2018) which acts as an attractor within the edge.

1.4. Relative periodic orbits

We hypothesise that propagating structures may be RPOs, or at least remain in the vicinity of such orbits. Relative periodic orbits are configurations which return to their initial state after some time, up to the application of a symmetry operator allowing, e.g. for propagating pulsations. This broad class of solution includes the simpler cases of steady states, travelling waves (steady states in a moving frame of reference) and stationary periodic orbits.

The configuration (system state) of the plasma will be represented symbolically by a quantity $\mathbf {u}(t)$ which varies in time $t$ (this state might, for example, be a velocity field $\mathbf {v}(\mathbf {R},t)$ ). A periodic orbit is a configuration for which $\mathbf {u}(\tau +t) = \mathbf {u}(t)$ for some $\tau$ .

Mathematically, RPOs are roots of

(1.1) \begin{eqnarray} G(\mathbf {u};\textbf { k},\textbf { K},T)\equiv Q(\textbf { k}, \textbf { K}) Y^T\mathbf {u}-\mathbf {u}, \end{eqnarray}

where $Q(\textbf { k},\textbf { K})$ is a symmetry operator parameterised by the continuous ( $\mathbf {k}$ ) and discrete ( $\mathbf {K}$ ) parameters, $Y^T$ is the time evolution operator that captures the effect of applying the equations of motion for time $T$ . For any $\mathbf {u}$ , we can evaluate $G$ using standard initial-value simulation codes (e.g. by time evolving the distribution function using the Vlasov equations). With the trivial (identity) symmetry operator $Q=1$ , this equation reduces to states that repeat themselves after time $T$ . Both $\mathbf {u}$ and symmetry parameters $\mathbf {k}$ are varied to find roots and we denote $\mathbf {w}=(\mathbf {u},\mathbf {k})$ .

Finding roots of high-dimensional, nonlinear functions such as $G$ is a notoriously difficult problem. We follow the approaches laid out by Viswanath (2007) and Duguet et al. (2008). The classical method for solving these problems is to iteratively apply Newton–Raphson to improve an initial guess

(1.2) \begin{equation} \mathbf {w}^{n+1}=\mathbf {w}^n+\delta \mathbf {w}^n, \end{equation}

where

(1.3) \begin{equation} \mathbf {J}(\mathbf {w}^{n})\delta \mathbf {w}^{n}=-\mathbf {G}(\mathbf {w}^k), \end{equation}

with $\mathbf {J}$ being the Jacobian matrix for $\mathbf {G}$ . The large dimension of such problems renders calculating the Jacobian explicitly at best time consuming and frequently impossible. A Newton–Krylov method is therefore utilised to avoid explicitly calculating it. To this end we use the Generalised minimum residual algorithm (Saad & Schultz 1986) which gives an approximation to the solution of (1.3) using only matrix-vector products which can be found directly.

The second difficulty with finding the roots of nonlinear equations is the necessity of finding a good initial guess. The tightness of this constraint is relaxed somewhat by using a so-called globally convergent method (here the double dogleg step method (Dennis & Schnabel 1996)) which expands the basin of attraction for the root-finding algorithm. Once an RPO has been found, it can then be tracked through parameter space by slowly altering one parameter and using the RPO from the previous value as an initial guess.

2. The plasma interchange model

There are a wide variety of stylised models of plasma turbulence, but one of the simplest is the plasma interchange (PI) model of plasma instability derived from resistive magnetohydro-dynamics (MHD) (Beyer & Spatschek 1996; Benkadda et al. 2001). The underlying two-dimensional (2-D) resistive MHD model in the $(x,y)$ -plane has a curvature-driven instability perpendicular to the magnetic field, giving rise to a Rayleigh–Taylor instability that transports thermal energy in the $x$ direction, with a dynamics identical to a 2-D buoyancy-driven fluid model (the Rayleigh–Benard–Couette equations). This 2-D model is not directly proposed as a model for drift-wave turbulence (i.e. unlike microturbulence models it is unstable at the system scale); instead we use a quasi-1-D reduction of the model to capture certain of the dynamics of gyrokinetic turbulence, which was shown to persist when only a single $k_y$ (and the zonal dynamics) were retained.

The 2-D model is reduced to a quasi-1-D model by representing only two Fourier components in the $y$ direction; the constant (zonal) component and a finite-wavenumber drift mode. This model is selected as a minimal model that avoids certain plasma-specific details. For example, the linear modes are zero real frequency, whereas slightly more complex models like Hasegawa–Wakatani have waves with finite group velocity, allowing propagation of structures in the linear limit (McMillan et al. 2009; Zhou et al. 2020). The hypotheses that motivate this model is that much of the phenomenology seen in complex shaped plasmas is just a neutral fluid dynamics, and that it is the interaction between zonal and non-zonal features that is the key nonlinearity.

Earlier investigations of the PI model (McMillan et al. 2009; Pringle et al. 2017) show that it exhibits dynamical features of drift-wave turbulence, such as zonal flow generation and the formation of avalanche-like structures (McMillan et al. 2009; Van Wyk et al. 2016). For high levels of flow shear any initial condition returns to a laminar (zero-turbulence) state. Below a critical shear value a sustained non-zero flux is possible and patches of spatially localised turbulence appear (figure 1 b). Reducing the shear still further leads to the turbulence spreading to fill more of the domain. For some intervals in shear the turbulent state is transient and back transitions to the laminar state are observed. As the shear becomes still smaller the turbulence ultimately becomes domain filling.

The edge of chaos in the PI model has previously been found to contain RPOs (in fact travelling waves (Pringle et al. 2017)), which are straightforward to isolate using a bisection method. This provides a starting point, and these RPOs may then then tracked through parameter space using the Newton–Krylov method.

The PI model evolves direct space functions ${\tilde {n}}(x,t)$ , $\tilde {\phi }(x,t)$ , $\bar {n}(x,t)$ and $E(x,t)$ , defined on a periodic radial domain $x \in [0, L)$ (which is the 1-D reduction of the 2-D shearing-box domain). Here, $\bar {n}$ and $E$ are real, describing the zonal density and zonal flow, respectively, with spatial boundary conditions $\bar {n}(L,t) = \bar {n}(0,t)$ and $E(L,t) = E(0,t)$ . The parameters $\tilde {\phi }$ and $\tilde {n}$ describe the drift-wave electrostatic potential and density, respectively (they represent a single $k_y$ plane wave in the $y$ direction, which is not explicitly modelled), with boundary conditions ${\tilde {n}}(L,t) = {\tilde {n}}(0,t) \exp (i S L t)$ and $\tilde {\phi }(L,t) = \tilde {\phi }(0,t) \exp (i S L t)$ ; these boundary conditions arise due to the large-scale sheared flow $S$ in the shearing box. Evolution equations are

(2.1) \begin{align} \partial _t {\tilde {n}} = - i {\tilde {n}} (E+Sx) + i \tilde {\phi } ( \partial _x \bar {n} -1 ) + \partial _{xx} {\tilde {n}}, \\[-24pt] \nonumber \end{align}

(2.2) \begin{align} \partial _t \tilde {\phi } = - i \tilde {\phi } (E+Sx) + i {\tilde {n}} + \partial _{xx} \tilde {\phi }, \\[-24pt] \nonumber \end{align}

(2.3) \begin{align} \partial _t \bar {n} = i \partial _x ({\tilde {n}}^* \tilde {\phi } - {\tilde {n}} \tilde {\phi }^* ), \\[-24pt] \nonumber \end{align}

(2.4) \begin{align} \partial _t E = i \partial _x ( \tilde {\phi } \partial _x \tilde {\phi }^* - \tilde {\phi }^* \partial _x \tilde {\phi } ) + \nu ' \partial _{xx} E, \end{align}

representing the instability (linear terms in $\tilde {n}$ and $\tilde {\phi }$ evolution equation), diffusion of the zonal potential and coupling between the zonal and non-zonal fields. There is a single control parameter $S$ representing the background flow shear amplitude.

The model is solved using a Fourier method in the $x$ direction and a second-order implicit solver in time. For zero $S$ , a Rayleigh–Taylor instability arises, with $\tilde {n}$ and $\tilde {\phi }$ growing exponentially in time, and maximum growth rate at zero $k_x$ (corresponding to motion aligned with the driving gradient).

Although these equations have no explicit time dependence, the boundary conditions mean that the simulation states can only be directly compared at time intervals $\delta t = 1/S L$ , because this is when they are in the same function space (this is also the case in local gyrokinetic simulations with background flow shear). Detail of the derivation of these equations and the qualitative behaviour of solutions is in earlier work (McMillan et al. 2009; Pringle et al. 2017; Smith 2023). The box length $L$ needs to be chosen larger than typical correlation sizes for the same reasons as in a radially periodic gyrokinetic system; the radial periodicity is artificial and the physical limit is $L\rightarrow \infty$ . However, smaller $L$ makes finding RPOs simpler. We use $L=50$ here, which appears to be a good compromise based on investigations in Smith (2023).

Figure 2. Top: The RPOs identified in the PI model, plotted as circles, with their stability labelled by colour, with the time average (over one orbit) of the spatially averaged flux for a series of turbulent simulations at different $S$ in the $L=50$ sized box. The maximum and minimum flux (over time) of these orbits are represented using whiskers (note that these are not error bars). The envelope of the time evolution of an initial-value simulation is plotted as a light-blue region, and a smoothed average as the intermediate blue trace. Red labels (a) and (b) identify the RPOs plotted in figures 3(a) and 3(b) respectively. Bottom: spectral density of the same turbulent simulations against frequency and $S$ .

There are a number of symmetries of this system. The system has discrete time translation invariance. It also possesses a continuous radial ( $x$ ) translation symmetry, complex phase invariance in joint rotation of $n$ and $x$ and invariance when a constant is added to $E$ . The phase invariance is inherited from the 2-D model, where it corresponds to invariance to translation in the $y$ (binormal) direction. It also possesses a (discrete) mirror symmetry under $x \rightarrow -x$ and $\phi \rightarrow \phi ^*, n \rightarrow n^*$ (which means burst propagation does not have a preferred direction (McMillan et al. 2009)). The symmetry operator $Q$ is then the composed operation of all these possible symmetries, with a specific choice of associated discrete and continuous parameters. For example this might correspond to applying a mirror symmetry as well as a displacement $x \rightarrow x-1$ .

Figure 2 (upper) shows x-domain-averaged fluxes ( $\langle i {\tilde {n}}^* \tilde {\phi } - i {\tilde {n}} \tilde {\phi }^* \rangle _x$ ) of the initial-value simulation and of RPOs; the stability of the RPOs is also presented. We first describe the value of initial-value simulations initialised near the critical shear value $S=1.6$ , with the shear parameter $S$ slowly reduced as the simulation is evolved forwards in time. The light blue line shows the long-time average simulation flux for a range of shear values. The blue-shaded region is the observed range of that flux over time. Above $S\sim 1.125$ average flux remains relatively invariant while the range of the flux gradually decreases. The range then contracts more rapidly around $S=1.46$ and then gradually decreases to zero at $S\sim 1.57$ . Above this, flux is time constant (and the simulation state is a simple travelling wave) until a sudden transition to laminar flow occurs at $S\sim 1.62$ .

We then search for RPOs in this system of equations using the Newton–Krylov solver, initialised with solutions derived either from the standard initial value or an edge-state solution (found using bisection (McMillan et al. 2018)). Once an RPO has been successfully identified, we may attempt to find a related RPOs for slightly different control parameter values $S$ , and given a pair of nearby solutions we can attempt to linearly extrapolate to provide an initial guess for a third solution; this allows us to track RPOs through phase space. The stability of these RPOs is of significant interest; this may be characterised in terms of the dimension of the space spanned by the unstable eigenvectors of the map $G$ linearised for small $\mathbf {u}$ (we can afford to construct the full linearised map as a matrix; Krylov methods would be needed for higher-dimensional models). This describes how stable the RPO is to small perturbations. If all eigenvalues are stable (bar those in the symmetry directions, which will have zero eigenvalue), the RPO is an attractor and we might expect to see it arise as a late-time solution to the initial-value problem. Higher numbers of unstable eigendirections (eigendirections are the linear spaces associated with the eigenvectors of the linear mapping) indicate that the RPO is not a stable steady-state solution and the initial-value simulations would approach this solution at most transiently. The nature of the eigenspace is of interest too; when complex-conjugate pairs of eigenvalues switch from positive to negative, that indicates a (sub- or super-critical) Hopf bifurcation. A Hopf bifurcation, for example, can give rise to a time-invariant steady state evolving into an oscillatory limit cycle of increasing amplitude as the control parameter is increased. On the other hand, in other classes of bifurcation, a stable solution can disappear without a nearby RPO appearing.

The behaviour of the initial-value simulation may be explained in terms of the family of RPOs identified. The simplest of these are travelling wave solutions (nonlinear states that appear steady in an appropriate co-moving frame of reference). Consider the RPOs depicted in orange, with one unstable eigendirection, on the bottom right of the upper frame of figure 2: these RPOs have been found for $S \in [0.2, 1.62]$ (the branch visible in figure 3 a extends to low $S$ (Smith 2023)) with increasing flux as $S$ increases: direct tracking of the RPOs allows us to follow the solution around the corner, where the amplitude of these states increases as $S$ decreases. The lower branch solution has earlier been identified using a bisection technique and is the attractor within the edge of chaos (Pringle et al. 2017); its one unstable direction pointing out of the edge of chaos – up to turbulence and down to the laminar flow – and within the edge of chaos it is a stable attractor. At $S=1.6$ the number of unstable directions associated with this travelling wave (TW) switches from one to zero, this indicates a saddle node bifurcation in the parameter $S$ . The upper branch (figure 3 a) is stable for a short range of values of $S$ , and it is this stable branch that turbulence was invariantly tracking before the sudden relaminarisation at the saddle node.

Figure 3. Density plots of zonal electric field versus time and spatial position in the PI model for the travelling wave RPO at (a) $S=1.579$ and full RPO at (b) $S=1.441$ (see labelled positions on figure 2). Only part of the full ( $L=50$ ) $x$ -domain is shown.

Tracking the upper branch back in $S$ , the TW undergoes two changes in stability in quick succession (just to the left of the point labelled “a” in figure 2), each corresponding to a pair of complex-conjugate eigenvalues becoming unstable - two Hopf bifurcations. These Hopf bifurcations give rise to new RPO solutions (in the range $1.47 \lt S \lt 1.57$ ) closely related to the TW solutions but now with a periodic variation whose period matches that of the eigenvalues that become unstable here (figure 3 b). In this range of $S$ , these RPOs have average flux corresponding to the average flux of the observed turbulence and a flux range also matching the observations. A series of Hopf bifurcations is a natural path to complex behaviour in frequency space (figure 3 b), because the combination of multiple frequencies approaches the broadband signal expected in a fully chaotic system.

Oscillatory RPOs are more difficult to track from a numerical perspective as the shearing-box boundary conditions in this model implies RPOs must have a period an exact multiple of the box period; there is a “base” period of 2 units seen in the spectrum, and RPOs near this period can only be found at a small number of isolated values of $S$ where it matches the boundary condition period. Potentially there are more solutions at multiples of the base period at intermediate $S$ values. Just below $S = 1.47$ the family of RPOs tracked no longer corresponds to the dynamics seen in the initial-value dynamics, there is a larger oscillation in the initial-value dynamics than the RPO; the oscillatory RPO, which is stable at higher shear, undergoes two Hopf bifurcations near $S=1.47$ , so the divergence is unsurprising. We were not able to identify additional nearby RPOs in the parameter regime near the point labelled “b” that corresponded to these orbits. Note that, at this point, multiple frequencies enter the spectrogram, and the dominant two frequencies do not appear to be divisors of a common larger frequency. This suggests that the initial-value simulation is not following an RPO, but possibly undergoing some more complex quasiperiodic behaviour with two different time periods, tracing out a 2-D torus in configuration space. Nonetheless, for $S \in [1.28,1.32]$ , we were able to identify a new family of RPOs matching the dynamics of the turbulent simulations (see Smith (2023) for details of the structures).

Physically, the onset of oscillation in the dynamics occurs as the TWs begin shedding a train of zonal flow oscillations. This is related to the secondary instability (Rogers et al. 2000) that is frequently cited as the origin of zonal flow formation; here, the zonal flows occur as a temporal modulational instability to an saturated nonlinear state. As the background flow shear decreases, nonlinearly generated zonal flows increase, maintaining the overall flow shear, and providing a route to saturated turbulence.

The increasing complexity of the initial-value state as shear decreases may be illustrated by plotting the frequency spectrum of the flux at each shear value (figure 2, lower). The Hopf bifurcation at $S=1.58$ gives rise to a narrow signal at frequency $\omega =3.5$ and harmonics (not plotted). Near $S=1.48$ , an additional low frequency signal appears, at $\omega =0.92$ , along with harmonics of this frequency, and, at slightly lower $S$ , combinations of these frequencies, creating a rich spectrum. The amplitude of the initial $\omega =3.5$ oscillation reduces significantly at lower shear, so that the $\omega =0.92$ spectral line and harmonics dominate. At $S=1.16$ , a sudden transition away from this state occurs and a switch to a state with dominant frequency $\omega =1.1$ , but with background of broadband activity possibly indicating the onset of chaos. Overall, this process is reminiscent of generic transitions to chaos that occur due to the successive onset of instabilities in, for example, the period-doubling route to chaos.

3. Turbulent dynamics in gyrokinetics

In order to examine the transition from simple propagating states to gradually fully developed domain-filling turbulence in gyrokinetics, and make the connection to the earlier PI model results, we perform simulations using the GENE (Dannert & Jenko 2005) code in local triply periodic geometry, where the x (radial), z (along the field line) and y (bi-normal to x and z, toroidal-like) coordinates are periodic. Usually, the $z$ coordinate has a more complex periodicity (Beer et al. 1995) but we run with zero magnetic shear which has the advantage for our purposes of allowing a continuous radial translation symmetry (McMillan et al. (2018) shows that the absence of this symmetry leads to rapid oscillations in the travelling disturbances). Apart from the choice of zero magnetic shear (and a slightly lower density gradient), these simulations use CYCLONE (Dimits et al. 2000) (a standard simplified tokamak test case) parameters, and are electrostatic adiabatic-electron collisionless runs in a concentric circular geometry. A sheared poloidal background zonal flow is introduced, with the poloidal shearing rate parameterised by $S = (R/c_s) d\lt v\gt /dr$ , with $c_s/R$ the ion sound transit time. Note that, formally, the simulation geometry is then a “shearing box”, like the PI model (McMillan et al. 2019).

Standard-width simulations have $y$ domain length $L_y =60 \rho _s$ in line with typical practice; for the purpose of identifying RPOs, a subspace of a standard simulation is selected by running narrow simulations with size $L_y = 20 \rho _s$ in the binormal direction, which permits only one significantly unstable toroidal mode at $k_y \rho _s =\,0.3$ to develop in the simulation (as well as zonal fluctuations). This choice is prompted by fairly strong peaking of the flux spectrum near $k_y \rho _s = \, 0.3$ , especially at large flow shear, and earlier simulations (McMillan & Sharma 2016), which found that narrow-box simulations with this choice of wavenumber gave similar qualitative and quantitative results. The use of a restricted toroidal domain probes whether drift-mode/drift-mode coupling is important; this coupling is absent in the single-mode simulations. Restricted domain size also allows us to find RPOs that may be unstable in the larger simulation box. Otherwise parameters are x domain length $L_x = 560$ , safety factor $q=1.4$ , zero magnetic shear, local inverse aspect ratio $r/R=0.18$ and background temperature and density gradients ${\rm d}T/{\rm d}r = 6.9 T/R$ and ${\rm d}n/{\rm d}r = 2.0 n/R$ , respectively (note that CYCLONE has ${\rm d}n/{\rm d}r = 2.2 n/R$ ). Twenty-four grid points were used in $v_{||}$ , $8$ in $\mu$ , $768$ in $x$ , 16 in $z$ and $1$ and $15$ (non-zonal) modes in $k_y$ treated for the narrow- and standard-box simulations repectively.

Figure 4. (a) Density plot of electrostatic potential (in units of $T_0 \rho _s/e R$ ) at the outboard mid-plane verus radial ( $x$ ) and binormal ( $y$ ) coordinate for standard-box gyrokinetic simulations for $S=0.35$ . This is the state at $t=250 a/c_s$ in this simulation (cf. Figure 1 a). (b) Density plot of electrostatic potential at the outboard mid-plane versus radial ( $x$ ) and binormal ( $y$ ) coordinate for standard-box gyrokinetic simulations for $S=0.75$ . In both cases the bursts are propagating in the positive $x$ direction.

For finite shear, the simulations are linearly stable, but may be nonlinearly unstable even for very small amplitudes. Initial conditions for the simulations are chosen so that a localised turbulent patch forms by slightly increasing the amplitude of the ’edge of chaos’ state found in earlier work (McMillan et al. 2018). The choice of perturbation is less important at lower shear values $S \lesssim 0.5$ , but as shear approaches the critical value, it becomes more difficult to initialise turbulence (Casson 2011). Whether real tokamaks will see finite flux in regimes where the turbulent attractor is small is unclear.

In narrow-box simulations the flow shear, $S$ , is slowly scanned downwards from $S=0.8$ (figure 5, black line), where the state is a single, time-invariant TW in the solution. These are isolated structures in the $x$ direction with typical extent in $10\rho _s$ in this direction, with a strong associated perturbed zonal shear flow of similar magnitude to the background shear flow. These TWs are qualitatively the same as those found in similar simulations (with the same physical parameters except magnetic shear $\hat {s}=0.76$ ) of McMillan et al. (2009) and McMillan & Sharma (2016), and the restriction to a single non-zero $k_y$ means that the TWs cannot localise in the $y$ direction unlike those seen in Van Wyk et al. (2016).

At around $S=0.725$ , this state splits into a pair of TWs. At $S=0.71$ , the system transitions to an oscillating state, which is seen in the spectrum as a coherent peak and a second harmonic; the oscillation has a roughly 30 % peak-to-peak oscillation, so the appearance of the second harmonic is not surprising. Below about $S=0.53$ , the simulation dithers between a state which is just below $1c_s/R$ and one which which has a fundamental oscillation frequency just above $2c_s/R$ (near, but not exactly at the harmonic), switching finally to the low frequency state below $S=0.46$ until an abrupt transition to domain-filling turbulence at $S=0.38$ . The details of this dynamics are unclear, but the overall time-averaged structure of the state does not dramatically change.

Figure 5. Top: evolution of the flux versus shear in narrow-box (black) and standard-box (red) gyrokinetic simulations where the shear is slowly varied. The narrow-box scan is an downwards scan in shear, and the standard-box scan is upwards in shear to the first back transition to laminar. The standard-box simulation is restarted at higher shears where there are narrow shear windows with propagating turbulent puffs, outside which only laminar solutions were found. Bottom: density plot of spectral amplitude of flux versus shear $S$ , and frequency, for the narrow-box gyrokinetic simulation with decreasing shear. Narrowband features in the spectrum are evident in the region $0.38 \lesssim S \lesssim 0.72$ where the dynamics consists of a single turbulent puff, but disappears above $S \sim 0.72$ where the propagating structure is a time-invariant TW.

At this point, the TW splits up and we see domain-filling turbulence arising, and an associated broadband frequency spectrum. An upwards scan from $S=0.8$ find that the simulation transitions back to laminar at $S=0.81$ . We find that the stable state above $S=0.75$ and the (unstable) “edge of chaos state”, found by the bisection method (McMillan et al. 2018), have structure and properties that converge as $S$ approaches the critical value where both solutions vanish. That is, the edge of chaos state and the stable state at high $S$ are two branches originating from a saddle node bifurcation, just as in the PI model system.

In simulations with standard binormal resolution (i.e. the physically relevant regime), the quasi-periodic states found in the narrow-box simulation are generally found to be unstable to modulational instabilities. This leads to a state that is not monochromatic in $k_y$ , but usually with only two significantly non-zero $k_y$ components, so it is weakly localised in the $y$ direction (figure 4 b). At $S\gt 0.62$ this state usually decays back to the laminar state. On the other hand, in the interval $S\in [0.52, 0.62]$ the state remains qualitatively similar to that found in the narrow box. Even at lower shear, where the burst has split into two, near the transition to space-filling turbulence, the state is dominantly monochromatic in $y$ (figure 4 a). This finding, along with qualitative similarities in behaviour (as found in earlier investigations on the edge of chaos state) suggests that a more complex dynamics in the large box is structured around the periodic states found in the small simulation box. That is, the turbulent attractor is structured around the stable RPO found in the narrow-box simulation. Simulations run initialised near the RPO found in the narrow-box simulation slowly diverge away, even though this is also (by symmetry) an RPO in the wide-box simulation, it is unstable to a modulational instability which causes an amplitude modulation in the $y$ direction.

The RPOs found in the small box thus play a conceptually similar role to plane waves; states with high degrees of symmetry, not generally found as a pure state in nonlinear simulation, but that nonetheless provide a vital description of the processes underlying the dynamics. Unlike linear theory, the RPOs also encode the nature of the nonlinear saturation process. In this case, as the RPOs involve a single binormal mode coupling with zonal behaviour, the nonlinear saturation process is interaction between the drift wave and zonal fluctuations.

4. Conclusion

In all the systems considered, the onset (in the sense of decreasing $S$ ) of turbulence is linked to the emergence of nonlinear solutions via a saddle node bifurcation. This is shown explicitly for the PI model and via careful initial-value simulations for GK turbulence. The high energy solution that emerges from this saddle node describes the initial stages of turbulence either exactly (PI model) or closely (GK simulation). As $S$ decreases further, the increasing complexity of the dynamics is linked directly to the decreasing stability of these solutions. This sequence of increasing numbers of nonlinear solutions underlying the onset of turbulence is analagous to that seen in classical shear flows (Gibson et al. 2009; Budanur et al. 2017), where the critical value for sustained turbulence is related to the increasing number of RPO orbits that provide a skeleton for a sustained turbulent state.

Determining where sustained turbulence is possible, and understanding transitions between various kinds of plasma states, is a longstanding problem of broad interest in tokamak physics. The Low-to-high confinement transition, for example, is often discussed in terms of bifurcations in systems of differential equations (Kim & Diamond 2003), with turbulence represented in some mean-field fashion, but here we are able to explain transitions between turbulent states as bifurcations while explicitly modelling the turbulent structures.

These tools also provide a promising new technique to examine nonlinear zonal flow oscillations in plasmas that are not necessarily related to the linear dynamics of the equilibrium state (Villard et al. 2019), but may instead be fundamentally nonlinear. More generally, there are a range of coherent, sometimes propagating, nonlinear structures observed in plasmas where RPOs may be useful as a conceptual framework, and directly obtaining these solutions via numerical root finding may be fruitful. Directly obtaining solutions is more systematic than attempting to isolate them “by eye” in post hoc analysis, and stability analysis explains links between related states.

We have presented a method to isolate classes of relative periodic solutions applied to tokamak turbulence, and isolated a series of these nonlinear states in a simple plasma model. These states are seen to form a skeleton about which the transition to domain-filling plasma turbulence occurs, and to explain the prevalence of propagating structures in this system. Increasingly complex time dynamics is allowed as a pure propagating state destabilises, gives rise to a propagating oscillating state, which then itself becomes unstable, reminiscent of the period-doubling transition to chaos.

Qualitatively similar behaviour is found in gyrokinetic simulations of turbulence where, under certain symmetry restrictions, RPOs are stable in a large parameter regime; these RPOs are less stable in a full-scale simulation, and so serve as organising centres for the more complex dynamics seen there. By identifying nonlinear solutions directly, a direct understanding of sustainment of turbulence in sub-critical turbulence may be obtained. Since the simplest RPOs here are directly connected to the edge-state solution seen in earlier work, this demonstrates that the same sustainment processes are in play, with the drift mode creating a self-destabilising sheared zonal flow.

Although these propagating structures in core tokamak turbulence have been a subject of considerable discussion in the literature, beyond the usual motivations of curiosity and verification they are not central to pragmatic efforts to achieve fusion energy. For example, the fluxes associated with these structures have gyro-Bohm scaling, and the energy releases associated with them are not a significant fraction of the reactor-scale device energy. On the other hand, carefully characterising these propagating structures was required to establish that they are events of a specific characteristic amplitude, not global avalanches causing occasional extreme flux events.

It is really the propagating structures in the edge of tokamaks that are of more obvious importance, and since these are caused by the intermittent dynamics inside the last closed flux surfaces, the dynamics of the pedestal is where attempts to characterise propagating and time-periodic structures are of obvious utility. Propagating fronts of turbulence in the pedestal, which may be related to RPOs, could be the root cause of blobs ejected into the edge. Some forms of periodic or quasi-periodic variations of the edge region might also be characterised directly as RPOs.