2 The system analysed

This strongly magnetised plasma system, an idealised tokamak, is described via a gyrokinetic (GK) framework (Hahm Reference Hahm1988). The GK equations describes particle motion in magnetised plasmas in a self-consistent electromagnetic field by evolving the gyrotropic particle distribution function $f(x,y,\unicode[STIX]{x1D703},\unicode[STIX]{x1D707},v_{\Vert })$ , where $x,\unicode[STIX]{x1D703}$ and $y$ are spatial coordinates and $\unicode[STIX]{x1D707}$ and $v_{\Vert }$ are velocity-space coordinates. The $x$ coordinate parameterises the radial direction, $\unicode[STIX]{x1D703}$ is the poloidal (straight field line) angle and the $y$ coordinate is a field-line label, with $y\propto (\unicode[STIX]{x1D701}-q(x)\unicode[STIX]{x1D703})$ , with $\unicode[STIX]{x1D701}$ the toroidal direction and $q$ the safety factor (Beer, Cowley & Hammett Reference Beer, Cowley and Hammett1995). Note that changing $y$ and keeping other spatial variables fixed corresponds to a displacement in the toroidal direction of symmetry. The $x$ and $y$ coordinate are scaled such that $|\unicode[STIX]{x1D735}x|,|\unicode[STIX]{x1D735}y|\sim 1$ .

Figure 1. Geometry of the plasma system near the outboard mid-plane. Black circles indicate convective vortices generated by the drift instability. The velocity arrows indicate the sheared background flow.

The basic instability driving turbulence in this system is the ion temperature gradient instability, driven by pressure gradients aligned with magnetic curvature. A sketch of the geometry is provided in figure 1. Simulating tokamak turbulence in a kinetic rather than fluid model allows certain details such as parallel Landau damping, perpendicular particle resonances and finite-Larmor radius effects to be retained. These features are essential to recover a good quantitative match against experimental data: qualitatively, however, many of the dynamical features match a highly simplified one-dimensional model of turbulence (McMillan et al. Reference McMillan, Jolliet, Tran, Bottino, P. and Villard2009; Pringle et al. Reference Pringle, McMillan and Teaca2017). For the GK system of interest, quantities peak near $\unicode[STIX]{x1D703}=0$ and vary slowly along the field line (with field-line spatial dependence resembling linear eigenmodes), and we believe much of the important nonlinear dynamics can be understood by examining the $(x,y)$ dependence on the $\unicode[STIX]{x1D703}=0$ mid-plane in our analysis.

The usual control parameter of interest in plasma microturbulence is the gradient driving the instability, and a normalised measure is the drive rate compared to the parallel streaming time. Due to all spatial scales being subject to Landau damping by parallel streaming, but only some being favourable for instability drive, a somewhat narrow range of wavelengths is strongly activated in a typical simulation (especially in the simplified system here with only one active particle species). Finer scale structures are generated through spatial and velocity-space mixing, but their influence on collective behaviour is reduced due to gyro-averaging and velocity-space integration. Physical or numerical dissipation provides a means of saturation for these fine structures, but quantities like the spectrum of excited modes are often not very sensitive to the value of this dissipation. This should be contrasted with the situation in fluid turbulence where the typical control parameter is the Reynolds number, which directly determines the size of the range of wavenumbers that are dynamically relevant.

Compared to fluids, the overall dissipative nature of kinetic systems is much more complex. Collisions (Balescu Reference Balescu1960; Lenard Reference Lenard1960; Landau Reference Landau1981) are ultimately responsible for setting the dissipation scale. However, in typical under-resolved ‘collisionless’ simulations, anomalous dissipation (Eyink Reference Eyink2018) due to the truncation of the physical nonlinear couplings can be seen as the cause for the effective removal of energy from fine scales. In addition, hyperviscosity terms, such as the ones employed in our work, can also play this role. They can be seen as the simplest form of large eddy simulation models, applied to a gyrokinetic system (Navarro et al. Reference Navarro, Teaca, Jenko, Hammett and Happel2014). Moreover, as the route to dissipation occurs in phase space (Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Plunk, Quataert and Tatsuno2008), is linked to the phase-space mixing (Watanabe & Sugama Reference Watanabe and Sugama2006; Tatsuno et al. Reference Tatsuno, Dorland, Schekochihin, Plunk, Barnes, Cowley and Howes2009) and the energy cascade (Navarro et al. Reference Navarro, Morel, Albrecht-Marc, Carati, Merz, Görler and Jenko2011; Watanabe et al. Reference Watanabe, Sugama, Nunami, Tanaka and Nakata2012; Teaca, Navarro & Jenko Reference Teaca, Navarro and Jenko2014), non-intuitive behaviours can emerge (Schekochihin et al. Reference Schekochihin, Parker, Highcock, Dellar, Dorland and Hammett2016). The dissipation for the GK problem occurs over a wide range of spatial scales, including relatively large perpendicular spatial scales (Hatch et al. Reference Hatch, Terry, Jenko, Merz and Nevins2011, Reference Hatch, Jenko, Bratanov and Navarro2014), so in conjunction with the injection of energy and subsequent weak energy cascade (Howes et al. Reference Howes, Dorland, Cowley, Hammett, Quataert, Schekochihin and Tatsuno2008), the overall turbulence dynamics may be considered strongly dissipative. For example, in the fluid model of McMillan et al. (Reference McMillan, Jolliet, Tran, Bottino, P. and Villard2009) and Pringle et al. (Reference Pringle, McMillan and Teaca2017), an explicit dissipation of order 1 in terms of the characteristic scales is applied, which is modelling the effect of kinetic processes; this dissipative system nonetheless qualitatively captures many of the aspects of the nominally collisionless gyrokinetic system.

3 Formulation details and symmetry

We use the GKW code (Peeters et al. Reference Peeters, Camenen, Casson, Hornsby, Snodin, Strintzi and Szepesi2009) to evolve the electrostatic local (flux-tube) gyrokinetic equations, with adiabatic electrons, in the presence of a background poloidal flow shear with shearing rate $S$ . The aim is to focus on a somewhat simplified system, in which analysis is relatively straightforward, rather than to perform a detailed device-relevant full-physics model. Our analysis treats this fairly standard simulation setup in some ways as a ‘black box’, so much of the detail presented in this section will not be referred to in later discussion. Despite this, we present the equations of this GK system for completeness, in direct space (rather than spectral) form; more details are provided in the code reference paper (Peeters et al. Reference Peeters, Camenen, Casson, Hornsby, Snodin, Strintzi and Szepesi2009). The dynamics is found by solving a Vlasov equation for the perturbation $\unicode[STIX]{x1D6FF}f$ to the distribution function,

(3.1) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FF}f}{\unicode[STIX]{x2202}t}=-{\dot{Z}}_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{Z}\unicode[STIX]{x1D6FF}f-{\dot{R}}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{R}(f_{0}+\unicode[STIX]{x1D6FF}f)-\dot{v}_{\Vert 1}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}v_{\Vert }}f_{0}+C(\unicode[STIX]{x1D6FF}f),\end{eqnarray}$$

where $f_{0}$ is the (Maxwellian) background distribution function, $Z$ represents the five-dimensional phase space $(x,y,\unicode[STIX]{x1D703},\unicode[STIX]{x1D707},v_{\Vert })$ , $R$ the spatial coordinates $(x,y,\unicode[STIX]{x1D703})$ , ${\dot{Z}}_{0}$ are the drift trajectories in the absence of the perturbing electrostatic field, ${\dot{R}}_{1}$ is the $E\times B$ drift and $\dot{v}_{\Vert 1}$ represents acceleration due to $E_{\Vert }$ . $C$ stands in for the collision operator(which is not used here), or for the numerical hyperviscosity.

For the zero- $\unicode[STIX]{x1D6FD}$ , weak-flow, electrostatic case of interest, for ions of charge $e$ , and equal ion and electron temperatures, the equations of motion may be written as

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle {\dot{Z}}_{0}=({\dot{R}}_{0},\dot{v}_{\Vert 0},\dot{\unicode[STIX]{x1D707}_{0}})=\left(v_{\Vert }\boldsymbol{b}+\boldsymbol{v}_{D}+\boldsymbol{v}_{E0},-\frac{\unicode[STIX]{x1D707}B}{m}\frac{\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B}{B^{2}},0\right) & \displaystyle\end{eqnarray}$$

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle {\dot{R}}_{1}=\boldsymbol{v}_{E} & \displaystyle\end{eqnarray}$$

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{v}_{\Vert 1}=\frac{1}{mv_{\Vert }}(-e[v_{\Vert }\boldsymbol{b}+\boldsymbol{v}_{D}]\boldsymbol{\cdot }\unicode[STIX]{x1D735}\langle \unicode[STIX]{x1D719}\rangle _{\unicode[STIX]{x1D6FC}}-\unicode[STIX]{x1D707}\boldsymbol{v}_{E}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B), & \displaystyle\end{eqnarray}$$

with

(3.5) $$\begin{eqnarray}\boldsymbol{v}_{E}=\frac{\boldsymbol{b}\times \unicode[STIX]{x1D735}\langle \unicode[STIX]{x1D719}\rangle _{\unicode[STIX]{x1D6FC}}}{B},\end{eqnarray}$$

where angle brackets with subscript $\unicode[STIX]{x1D6FC}$ denote a gyro-average, $\boldsymbol{v}_{E0}=\boldsymbol{E}_{0}\times \boldsymbol{b}/B$ and

(3.6) $$\begin{eqnarray}\boldsymbol{v}_{D}=\frac{1}{e}\left[\frac{mv_{\Vert }^{2}}{B}+\unicode[STIX]{x1D707}\right]\frac{\boldsymbol{B}\times \unicode[STIX]{x1D735}B}{B^{2}}.\end{eqnarray}$$

For the derivatives of the background distribution function we use

(3.7) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}v_{\Vert }}=f_{0}\frac{mv_{\Vert }}{T}\end{eqnarray}$$

and

(3.8) $$\begin{eqnarray}\unicode[STIX]{x1D735}f_{0}=-f_{0}\left(\frac{mv_{\Vert }^{2}+\unicode[STIX]{x1D707}B-3T/2}{T}\frac{1}{L_{T}}+\frac{1}{L_{n}}\right)\unicode[STIX]{x1D735}x,\end{eqnarray}$$

with

(3.9) $$\begin{eqnarray}f_{0}=\frac{n}{(2\unicode[STIX]{x03C0}T)^{3/2}}\exp -\frac{mv_{\Vert }^{2}+\unicode[STIX]{x1D707}B}{T},\end{eqnarray}$$

where $L_{n}$ and $L_{T}$ are the density and temperature gradient scale lengths. The local limit allows taking $n$ , $T$ , $L_{n}$ and $L_{T}$ constant. The background electric field $\boldsymbol{E}_{0}=xB_{0}S\unicode[STIX]{x1D735}x$ so that the shearing rate $\text{d}(\unicode[STIX]{x1D735}y\boldsymbol{\cdot }\boldsymbol{v}_{\boldsymbol{E}\mathbf{0}})/\text{d}x\sim S$ .

The quasineutrality equation relates the gyro-averaged charge density associated with $f$ to the perturbed electric field $\unicode[STIX]{x1D719}$ ,

(3.10) $$\begin{eqnarray}\int \text{d}Z\frac{e}{T}f_{0}\langle \langle \unicode[STIX]{x1D719}\rangle _{\unicode[STIX]{x1D6FC}}-\unicode[STIX]{x1D719}\rangle _{\unicode[STIX]{x1D6FC}}-\frac{ne}{T}\left(\unicode[STIX]{x1D719}-\langle \unicode[STIX]{x1D719}\rangle \right)+\int \text{d}Ze\langle \unicode[STIX]{x1D6FF}f\rangle _{\unicode[STIX]{x1D6FC}}=0.\end{eqnarray}$$

In the long-wavelength limit, explicit calculation of the first term on the right-hand side yields $(mn/eB^{2})\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}$ , which is the charge associated with the polarisation response of the ions. The second term represents the adiabatic electron response, with the angle bracket (without a subscript) indicating a zonal average (volume-weighted integration in $y$ and $\unicode[STIX]{x1D703}$ direction) as the electrons are bound to the flux surfaces and do not respond to zonal charge fluctuations.

The symmetries of these gyrokinetic equations may be found by inspecting the form of these terms. The electrostatic field and hence the $E\times B$ drift can be expressed as a linear function of $\unicode[STIX]{x1D6FF}f$ , so we have ${\dot{Z}}_{1}=L(\unicode[STIX]{x1D6FF}f)$ . Due to axisymmetry and the flux-tube limit, both $L$ and ${\dot{Z}}_{0}$ are spatially invariant to translations in the $x$ and $y$ directions. The boundary conditions in the $x$ and $y$ directions are doubly periodic, but in the $\unicode[STIX]{x1D703}$ direction, there is a twisted periodicity of the form $f(x,y,\unicode[STIX]{x03C0})=f(x,y+sx,-\unicode[STIX]{x03C0})$ , with $s$ dependent on magnetic shear. As a consequence, the overall system has a continuous symmetry in the $y$ (toroidal) direction, but only a discrete (not continuous) translation symmetry in the $x$ (radial) direction. There is also an inversion symmetry (Parra, Barnes & Peeters Reference Parra, Barnes and Peeters2011), which involves changing the sign of one of the parameters (the flow shear), and allows us, as usual, to consider only cases with $S\geqslant 0$ , since results for $S<0$ may be found using the symmetry. For example, propagating structures exist with the same radial velocity but opposite sign for $S<0$ .

In the following, units for amplitudes use the local gyrokinetic convention, so the electrostatic potential $\unicode[STIX]{x1D719}$ is in units of $\unicode[STIX]{x1D719}_{0}=\unicode[STIX]{x1D70C}^{\ast }e\unicode[STIX]{x1D719}/T_{e}$ , with $\unicode[STIX]{x1D70C}^{\ast }=\unicode[STIX]{x1D70C}_{0}/R$ , so fluctuations, although order $1$ in these plots, are small in terms of relative density (here, $T_{e}$ is the electron temperature $\unicode[STIX]{x1D70C}_{0}=(mT_{e})^{1/2}/qB$ is the ‘ion sound gyroradius’, $R$ is the tokamak major radius).

The simulations use the standard set of CYCLONE parameters, with $L_{n}=R/2.2$ , $q=1.4$ , $(\text{d}q/\text{d}r)r/q=0.8$ , but with a slightly reduced temperature gradient, $L_{T}=R/6.0$ , and a concentric circular equilibrium, with local aspect ratio $0.18$ . The size of the simulation box is $[L_{x},L_{y}]=[157,84]\unicode[STIX]{x1D70C}_{0}$ , with 20 toroidal modes used, $321x$ grid points and $16,16$ and $32$ grid points used in the $\unicode[STIX]{x1D703}$ , $\unicode[STIX]{x1D707}$ and $v_{\Vert }$ directions. A normalised forth-order numerical hyperviscosity parameter of $0.1$ is chosen in the parallel, $v_{\Vert }$ and $x$ directions (in addition to inevitable numerical diffusion); this corresponds to damping oscillations at the grid scales in these directions on a time scale of $10t_{0}$ , and helps to avoid numerical problems of spectral pile up at fine scales without unduly influencing the longer scales that will be of interest here.

Internally, GKW represents the distribution function using a Fourier series. Except where specified otherwise, the simulations use an initialisation of the form

(3.11) $$\begin{eqnarray}f(k_{x},k_{y},\unicode[STIX]{x1D703},v,\unicode[STIX]{x1D707})=Af_{0}(v,\unicode[STIX]{x1D707})\exp \left\{\!-\frac{(k_{x}-k_{x0})^{2}}{C_{x}^{2}}-\frac{(k_{y}-k_{y0})^{2}}{C_{y}^{2}}\right\}\frac{1}{2}\left(\cos (\unicode[STIX]{x1D703})+1\right)\end{eqnarray}$$

representing a field-aligned density fluctuation with typical wavenumber $(k_{y0},k_{x0})$ modulated by a Gaussian envelope in the $x$ and $y$ directions with width parameterised by $C_{x}$ and $C_{y}$ , with an overall amplitude $A$ . Unless noted otherwise, values $k_{x0}\unicode[STIX]{x1D70C}_{0}=0.24$ , $C_{x}\unicode[STIX]{x1D70C}_{0}=0.1601$ , $k_{y0}\unicode[STIX]{x1D70C}_{0}=0.37$ , $C_{y}\unicode[STIX]{x1D70C}_{0}=0.074$ will be used.

4 The edge state

The black traces of figure 2 show the evolution of the heat flux in the system for initialisations of varying amplitude versus time (in units of transit frequency $t_{0}=c_{s}/R$ ). The linear system is stable despite periods where the flux increases in time, and simulations with sufficiently small amplitude initialisations decay. Given a large enough initialisation, however, sustained turbulence is triggered.

The amplitude of the initial perturbation was systematically varied, by a bisection technique (Pringle et al. Reference Pringle, McMillan and Teaca2017), to find the threshold amplitude below which the system decays to a laminar state, but above which it remains in a turbulent regime at late time. The simulations very close to threshold remain for some time near the separator between the stable and unstable manifold in the system, i.e. the edge of chaos, before diverging away. In figure 2 the near-stationary flux ( $\log _{10}$ [flux] ${\approx}-1$ ) of the edge state indicates that the edge dynamics is considerably simpler than the turbulent dynamics. The ‘steps’ that appear in decaying simulations ( $\log _{10}$ [flux]  ${\approx}-4$ ) are associated with a time dependent (‘Floquet mode’) behaviour (Waltz, Dewar & Garbet Reference Waltz, Dewar and Garbet1998); in this case this dynamics is too slow to play a role in the transition to turbulence.

Figure 2. Heat flux versus time (gyroBohm units) for simulations with $S=0.12t_{0}^{-1}$ and successive initial condition amplitudes chosen using a bisection method to approach the critical amplitude. Red traces are restarted from $t=120$ , with the distribution function rescaled to track the edge state.

The edge of chaos represents the separatrix between the attractors for the laminar and turbulent dynamical states, and is an unstable manifold for the system. When the dynamics is restricted to the edge (by careful choice of initial trajectory), however, we find a local attractor within the edge, which we refer to as the edge state. To analyse this state, in addition to standard simulations, we use a series with a very small $y$ domain (narrow simulations), one fifth the size of the standard domain (in units of the thermal gyroradius $\unicode[STIX]{x1D70C}_{0}$ ). In the narrow simulations, the non-zonal component is dominated by the longest-wavelength mode that fits in the $y$ direction (at late time more than 90 % of the vorticity is in this mode). We use the narrow simulations to focus more directly on the relevant dynamics in a simpler system with fewer degrees of freedom. The properties of the edge state are qualitatively and quantitatively very similar for standard and narrow simulations.

We also considered simulations initialised from a white noise perturbation, where independent normally distributed pseudorandom numbers are loaded into the numerical grid of the distribution function as an initial condition, and multiplied by the background distribution function $f_{0}$ . Performing the same bisection method to search for the critical amplitude yields an effectively identical edge state (shifted in the $x$ direction). The insensitivity of the result to the initial perturbation is an indication that the edge state is in fact an attractor within the entire edge manifold.

Figure 3. (a) Mean $\unicode[STIX]{x1D719}^{2}$ (averaged over $y$ ) at the mid-plane versus time and position in the travelling-wave frame $x-vt$ for the edge state in narrow simulations with $S=0.12t_{0}^{-1}$ . A periodicity over $3.2t_{0}$ is visible. (b) Non-zonal potential $\unicode[STIX]{x1D719}$ at outboard mid-plane versus $x$ and $y$ for edge state at $t=120t_{0}$ for $S=0.12t_{0}^{-1}$ for (top) narrow and (bottom) standard simulations.

For narrow simulations, the edge state is found to be very close to a travelling wave. We determined the radial velocity $v$ of this translating structure using a linear fit of the $x$ position of the peak root-mean-square (r.m.s.) amplitude of the non-zonal potential $\unicode[STIX]{x1D719}^{2}$ versus time. Detailed inspection (figure 3 a), reveals a small time oscillation, with period ( $3.2t_{0}$ ) equal to the distance between lowest-order rational surfaces in the system (here there are 60 of these surfaces in the domain) divided by the travelling-wave velocity. This is a consequence of the fact that for finite magnetic shear, local gyrokinetics has a discrete, rather than continuous, spatial translation symmetry. The edge state for narrow simulations is thus a relative periodic orbit rather than an exact traveling wave. The r.m.s. variation from exact periodicity (when sampled once every $3.2t_{0}$ ) is found to be $0.5$  % over a period $64t_{0}$ . For the standard simulations, plots of zonal quantities solutions also suggest a near-periodic orbits on the 12-fold discrete radial symmetry of the larger system, with zonal averages similar to the narrow simulations.

Figure 4. (a) Mean (time averaged from 40 to $120t_{0}$ ) temperature gradient (blue) and zonal shear flow $\text{d}\langle E\rangle /\text{d}r$ (red), both normalised to background gradients, versus position $x-vt$ (in the travelling-wave frame) of the edge state, for $S=0.12$ , and both narrow (dashed) and standard (solid) simulations. (b) Mean (time averaged from 90 to $150t_{0}$ ) temperature gradient (blue) and zonal shear flow $\text{d}\langle E\rangle /\text{d}r$ (red), both normalised to background gradients, versus position $x-vt$ (in the approximate burst frame) for a turbulent simulation state, for $S=0.16$ .

Figure 5. Mean of squared non-zonal potential $\unicode[STIX]{x1D719}^{2}$ (averaged over $y$ ) at the mid-plane versus time and $x$ for (a) $S=0.15t_{0}^{-1}$ , (b) $0.16t_{0}^{-1}$ and (c) for an edge state with $S=0.12t_{0}^{-1}$ . In (a), long-lived turbulence is seen in a slowly expanding region centred around the excitation front. In (b), turbulence is excited transiently over a period of $100t_{0}$ , remaining localised near the travelling excitation front but then decays. The edge state (c) is much simpler and smoother, but also of considerably lower amplitude.

The quasi-travelling-wave solution in both narrow and standard simulations (figure 3 b) consists of a tilted finite $k_{y}$ travelling-wave mode, fed by the gradient drive, that produces a travelling zonal shear flow ahead (leftwards) of the pulse, that opposes the background shear flow (figure 4 a). The travelling perturbation strengthens the temperature gradient ahead of the pulse, and weakens the gradient behind, as expected from the energy transport equation, given the localised heat flux associated with the burst (the change in gradient in figure 4(a) is of comparable size to the background gradient). Those two mechanisms would be compatible with a travelling wave in either direction (McMillan et al. Reference McMillan, Jolliet, Tran, Bottino, P. and Villard2009; Pringle et al. Reference Pringle, McMillan and Teaca2017), but when the nonlinearity in the simulation is turned off, the $k_{y}\neq 0$ mode amplitude continues to propagate (not shown) in the same direction for $10a/c_{s}$ due to the group velocity, which depends on the mean $k_{x}$ value and thus the sign of $S$ (McMillan et al. Reference McMillan, Jolliet, Tran, Bottino, P. and Villard2009). Note that narrow simulations do not permit a vortex pair (Horton & Hasegawa Reference Horton and Hasegawa1994) advection mechanism, where a spatially localised ‘blob’ self-advects across the domain, as they are dominated by one $k_{y}$ mode (unlike $y$ -localised features seen elsewhere (van Wyk et al. Reference van Wyk, Highcock, Field, Roach, Schekochihin, Parra and Dorland2017)). Time snapshots (figure 3 b) of the mid-plane potential for narrow and standard simulations show similar tilting and localisation in $x$ for the two simulations, but despite close similarities in $y$ -averaged diagnostics, the standard simulation does not decay towards the narrow edge state. The spatial scale of the edge structure is of order $10\unicode[STIX]{x1D70C}_{0}$ in the $x$ direction, and the typical wavelength in the $y$ direction is approximately $16\unicode[STIX]{x1D70C}_{0}$ ; this is of similar scale to the wavelength of the most unstable mode, at $k\unicode[STIX]{x1D70C}_{0}\sim 0.5$ ; this is also comparable to typical scales in the fully turbulent simulations here and elsewhere (Dimits et al. Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000). The combination of linear physics and nonlinear interaction with zonal flows that set the relevant length scales in the turbulence physics (Plunk, Navarro & Jenko Reference Plunk, Navarro and Jenko2015) also appears to be responsible for sustaining the edge state. So the radial ( $x$ direction) scales of the edge state might also be estimated by considering the wavenumber of the secondary instability that drives zonal flows (Rogers, Dorland & Kotschenreuther Reference Rogers, Dorland and Kotschenreuther2000); this also gives a scale in the $x$ direction of approximately $10\unicode[STIX]{x1D70C}_{i}$ .

Even though the numerical resolution chosen is known to be sufficient to obtain converged results for turbulence simulations of this case, it is possible that the slightly different quantities of interest related to the edge state require higher numerical resolution. We therefore found the edge state (using the bisection method) in simulations with doubled resolution in each of the five phase-space directions, with the hyper-diffusivity in code units reduced by a factor of $16$ , for $S=0.12$ . Qualitatively, there are no striking differences observed (figure 6); the propagating edge state obtained for the high-resolution simulation has a mean squared amplitude, width, and propagation velocity within $9\,\%$ , $2\,\%$ and $1\,\%$ of that found in standard simulations (with these quantities averaged over the time period $45-75t_{0}$ ). We therefore conclude that these phenomena are very insensitive to the value of the numerical diffusivity in the simulation and numerical resolution in general.

Figure 6. Intensity of the non-zonal field $\unicode[STIX]{x1D719}$ (averaged over the $y$ -direction) versus $x$ for $t=60t_{0}$ at the mid-plane for the edge state found in a standard (blue trace) and doubled-resolution simulation (green trace).

5 Transition to a turbulent state

In typical simulations with a uniform shear flow, an isolated perturbation of sufficient amplitude produces a spreading region of turbulence. For sufficient shear flow (figure 5), the propagation of turbulence is entirely in one direction, and isolated propagating disturbances are seen in the simulation, described variously as avalanches and bursts in previous work (Candy & Waltz Reference Candy and Waltz2003; McMillan et al. Reference McMillan, Jolliet, Tran, Bottino, P. and Villard2009; van Wyk et al. Reference van Wyk, Highcock, Schekochihin, Roach, Field and Dorland2016). Propagating phenomena have been frequently observed in tokamak turbulence simulations, especially when a background shear flow is imposed. Although these features are also present when no overall background flow shear is imposed (McMillan et al. Reference McMillan, Jolliet, Tran, Bottino, P. and Villard2009), we see very clear propagating features for large shears, and, as in other works (van Wyk et al. Reference van Wyk, Highcock, Schekochihin, Roach, Field and Dorland2016), these become more isolated as the shearing rate approaches the critical value. The propagating edge state (figure 4 a) has some features in common with the late-time nonlinear state (figures 5 a and 4 b), such as similar propagation velocities (to within 10 %) and spatial scales (a comparison of typical amplitudes can be seen in figure 10). Both are associated with a moving turbulence front supporting a moving zonal electric field that destabilises the system in front of the turbulence front, so can both be seen as a travelling excitation wave. Near the critical flow shear beyond which turbulence is quenched, there structures are not particularly localised in the $y$ -direction (figure 7 a), unlike those seen in certain more detailed simulations (van Wyk et al. Reference van Wyk, Highcock, Schekochihin, Roach, Field and Dorland2016); the detailed structure of the propagating features is somewhat parameter dependent, and we present a second example of a propagating state in a lower aspect ratio, zero magnetic shear simulation (figure 7 b). The phenomenological commonalities between all these observations of propagating structures in gyrokinetic simulations with background shear flows are so striking that a common origin seems like the simplest explanation. Curved and elongated eddies, associated with a self-driven radially propagating zonal shear flow appear to be a ubiquitous phenomenon in tokamak turbulence simulations. Understanding what physical and numerical parameters control the details of these structures, however, still requires further study.

Figure 7. Non-zonal potential $\unicode[STIX]{x1D719}$ at outboard mid-plane versus $x$ and $y$ for (a) the standard simulation at near-critical shear value $0.15$ and (b) a low aspect ratio simulation with zero magnetic shear at a near-critical shear value.

As in neutral fluid simulation, if the shear is increased beyond a certain point (here, $S\sim 0.15t_{0}^{-1}$ ) we observe relatively long-lived turbulence that unpredictably decays to the quiescent state. It is clear from figures such as figure 5(a) that for large shear ( $S\gtrsim 0.1t_{0}^{-1}$ ), the excited region of turbulence has an overall drift, so that ‘puffs’ of excited turbulence travel through the system, returning to a locally quiescent state after the puff has passed. Unlike in, say, pipe flow turbulence (Wygnanski & Champagne Reference Wygnanski and Champagne1973), where these puffs travel in the direction of fluid flow, the bursts here travel either aligned or anti-aligned with the direction of the temperature gradient. In these simulations an unphysical periodic boundary condition is applied in the $x$ direction, so that the turbulence gradually fills the domain. We consider a simulation at $S=0.12t_{0}^{-1}$ using an ‘open’ boundary condition (in this case applying Dirichlet conditions to $f$ and the electrostatic potential), with the standard initial condition displaced in $x$ so that it peaks at $x=80x_{\text{MAX}}$ . Here, the system becomes quiescent after the puff travels to the boundary (figure 8 c,d). On the other hand, at lower flow shear the boundary conditions have less influence on the interior of the domain, and late-time behaviour is similar (figure 8 a–c). The sensitivity to boundary conditions is surprising in some sense because turbulent structures are very much smaller (in the $x$ direction) than the system size, and one might expect the local turbulence properties to mostly depend on local gradients, rather that the conditions at the $x$ boundaries. Nonetheless, the propagating bursts allow for patterns of activation to be set up that transmit information over longer length scales in the $x$ direction. We performed a scan in $S$ , and find that turbulence can be sustained over a wider range of flow shear values in a periodic simulation than a bounded simulation (figure 9) which may explain some of the differences in earlier benchmarks (Casson Reference Casson2011). The more effective quenching of the flux by background flow shear in the Dirichlet simulations does not appear to be a consequence of the specific initial conditions chosen; a simulation started with $S=0$ , and restarted after turbulence has attained a near steady state with $S=0.12$ also decays to the laminar state.

Figure 8. Mean $\unicode[STIX]{x1D719}^{2}$ (averaged over $y$ ) at the mid-plane versus time and position in a simulation with $S=0.06t_{0}^{-1}$ with (a) Dirichlet and (b) periodic boundary conditions and $S=0.12t_{0}^{-1}$ with (c) Dirichlet and (d) periodic boundary conditions.

Figure 9. Volume-averaged heat flux versus $S$ for Dirichlet versus periodic simulations.

6 Scaling of the transition threshold with background flow shear

The minimal seed is the initial simulation state of minimal amplitude that allows transition to a turbulent state in a subcritical system. The amplitude of the minimal seed may be seen as a ‘safety threshold’ below which all perturbations will eventually decay to the laminar regime. The amplitude of the minimal seed can also be used to quantify the degree to which linear and nonlinear processes can allow amplification of small fluctuations up to turbulent levels. Examining the minimal seed amplitude compared to the edge-state amplitude and the turbulence amplitude thus provides a quantification of important pieces of subcritical physics. In general the amplification factor from the minimal seed level to the edge-state amplitude is not equal to the transient linear amplification factor, and in many fluids the nonlinear processes allow overall amplification factors orders of magnitude larger than linear transient growth alone.

For a general initial perturbation, the value of the transition threshold depends on the functional form of the initialisation. We varied the parameters of a wavepacket-like initialisation to find the ‘most dangerous’ state with a minimal nonlinear instability threshold as an approximation to the ‘minimal seed’. It is in principle be possible to perform a complete minimisation (Pringle, Willis & Kerswell Reference Pringle, Willis and Kerswell2012) which optimises over all possible initial states, and we were able to do this for a drift-wave model (Pringle et al. Reference Pringle, McMillan and Teaca2017). However, in the gyrokinetic context, in would require writing an adjoint gyrokinetic solver and performing subsequent computationally demanding simulations. The choice of a wavepacket type initialisation (that is, equation (3.11)) is motivated by earlier study (Pringle et al. Reference Pringle, McMillan and Teaca2017) that found this to approximately capture the true minimal seed in a drift-wave model. We performed scans of initialisation parameters at a fixed shearing rate value $0.04$ to determine the values that allowed transition at the lowest $A$ value for the standard simulation, finding $k_{x0}\unicode[STIX]{x1D70C}_{0}=0.24$ , $C_{x}\unicode[STIX]{x1D70C}_{0}=0.1601$ , $k_{y0}\unicode[STIX]{x1D70C}_{0}=0.37$ , $C_{y}\unicode[STIX]{x1D70C}_{0}=0.074$ for the multi-mode simulations. In the narrow simulations other parameters are the same but we take $C_{y}\rightarrow 0$ .

To compare these state amplitudes here, we use two different measures. Because the minimal seed and edge states are radially localised, to compare amplitudes to the typical turbulent state, we use the maximum squared potential in $x$ (the ratio between these values should be relatively independent of the system size in $x$ unlike a spatial average measure). The other measure used is a global average vorticity. The transition threshold with shearing rate (figure 10) scales roughly like $\exp (-1/S)$ , except that for standard simulations at small $S$ the transition threshold drops more rapidly. Very small initial amplitudes produce instability in the small $S$ limit. The linear transient amplification also scales with $\exp (1/S)$ in these systems (Waltz et al. Reference Waltz, Dewar and Garbet1998).

For large enough shear ( $S\gtrsim 0.1t_{0}^{-1}$ ), the transition threshold found based on the wavepacket initialisation is actually slightly higher than the edge-state amplitude, when measured using the maximum measure (figure 10 a) (rather than the r.m.s. measure (figure 10 b)): since the true minimum seed would have a lower transition threshold than any other state, this demonstrates that the wavepacket initialisation is not the minimum seed in this norm. This is not surprising since the reasoning used to suggest a wavepacket-like minimum seed (Pringle et al. Reference Pringle, McMillan and Teaca2017) is clearly not valid except in the regime of large transient growth at low shear. Also, the amount of nonlinear transient growth depends quite strongly on which norm is used; this is expected even in linear problems, where the amount of transient amplification is a direct consequence of the choice of norm.

Figure 10. (a) Logarithm of maximum value of squared non-zonal potential $\max _{x}\langle \unicode[STIX]{x1D719}(x,y)^{2}\rangle _{y}$ and (b) of mean vorticity versus initial flow shear for several simulation phases, for narrow (red trace) and standard (blue trace) simulations. The amplitude of the critical state for the baseline initial perturbation shape is shown as solid traces, and the amplitude of the edge state is shown as dashed traces. At larger shearing rates $S>0.06t_{0}^{-1}$ these results for both simulation types are similar.

The edge-state amplitude gives an estimate of the amplitude for which the linear and nonlinear terms balance; this reduces with small flow shear. On the basis that the scaling of the transition threshold can be explained based on linear transient growth, the overall pathway for a near-critical mode to become unstable is hypothesised to be transient linear growth amplifying an initial perturbation, pushing it slightly beyond the edge-state amplitude, after which the unstable trajectory departing from the edge state allows access to the turbulent regime. The typical situation in neutral fluid experiments, is that the transition involves several stages of linear growth chained together as a result of nonlinear effects (Pringle et al. Reference Pringle, Willis and Kerswell2012). This more complex situation appears to arise for small flow shear in the gyrokinetic simulations, where the additional toroidal modes in the standard simulation allows transition to turbulence at lower initial amplitudes (and lower edge-state amplitudes) through coupling between non-zonal modes. The idea that scaling of subcritical thresholds in gyrokinetic systems (in that case for the maximum shearing rate at which turbulence can be sustained) can be found by considering linear transient amplification was suggested by the results of van Wyk et al. (Reference van Wyk, Highcock, Field, Roach, Schekochihin, Parra and Dorland2017). This also appears to be the case in our simulations, except at low shearing rates, where the details of the nonlinear dynamics becomes more important (as in neutral fluids).

A travelling-wave-type edge state is found for all shear rates for the narrow simulations and for $S>0.04t_{0}^{-1}$ for the standard simulations. The amplitudes of the edge state and the critical perturbation amplitude are not affected strongly by increasing the simulation box width for $S>0.06t_{0}^{-1}$ where the edge states are qualitatively similar, and the relevant nonlinearity in the critical transition to turbulence is the drift-mode/zonal-flow (and zonal gradient) interaction.

7 Conclusions and discussions

The behaviour of the edge of chaos is qualitatively similar to simple plasma-interchange (PI) model (Pringle et al. Reference Pringle, McMillan and Teaca2017), strengthening the thesis (McMillan et al. Reference McMillan, Jolliet, Tran, Bottino, P. and Villard2009) that qualitative features in the dynamics are the consequence of fluid-like behaviour, rather than details of tokamak geometry, or subtleties in the kinetic physics. Despite the complexity of GK model compared to neutral fluid models, the edge manifold contains a quasi-travelling-wave attractor, which is also seen in the PI model (note that this is attractive only within the edge manifold, not globally). The increasing amplitude of the edge state to levels comparable to the turbulent fluctuation level (figure 10 a) near the maximum background flow shear for which turbulence can be sustained, in conjunction with the relative simplicity of the edge state hints that, as in fluid turbulence theory, analysis of periodic orbits in plasma turbulence could be a powerful tool for understanding how and where turbulent states exist. The resemblances between the quasi-travelling wave in the edge of chaos, and the bursts seen in the turbulent state are notable, but as in the PI system, the nature of the relationship between these two phenomena is still unclear. We have found a way to estimate the transition threshold in this system, to quantify how robust the laminar state is against external forcings. The scaling of the transition threshold matches the PI model, except at very low shear, and the gyrokinetic system follows the same pathway to turbulence in this parameter space.

Propagating features seen in the turbulence (avalanches/bursts) have qualitative features that echo the travelling-wave edge state, with similar propagation velocities, but have stronger amplitudes and are more disordered. A simple state was seen in the limit where the flow shear was increased to just below the threshold for sustained turbulence in van Wyk et al. (Reference van Wyk, Highcock, Field, Roach, Schekochihin, Parra and Dorland2017): these investigations of the critical behaviour in these systems hint at the importance of periodic orbits in the critical dynamics of such systems. Because of the simplicity of the edge state, the mechanisms that allows the edge state to propagate could be illustrated in detail; the travelling wave destabilises the region in front of itself by removing the background shear flow and increasing the temperature gradient, and the tilting of the drift waves leads to a finite group velocity of the wavepacket-like finite $k_{y}$ modes. These propagation mechanisms appear to carry over to the avalanche/burst features in the fully turbulent state (McMillan et al. Reference McMillan, Jolliet, Tran, Bottino, P. and Villard2009). Long range propagation of features allows powerful nonlocality in these systems: at large flow shearing rate, the system is only convectively unstable, so at a fixed spatial location, the system will eventually return to a zero-flux state. On the other hand, there are a broad range of shearing rates (figure 10 a) for which a local perturbation 10 % as large as the typical turbulence level is required to initiate turbulence.

Ideas around the edge of chaos and exact solutions are well established in subcritical, neutral fluid problems. Some progress has also been made in applying them in astrophysical plasmas (Rincon et al. Reference Rincon, Ogilvie and Proctor2007). Here we have taken the first steps in applying them to study tokamaks. The results show that these methods can reveal intriguing aspects of this problem, but pose as many questions as they answer. Is the edge always dominated by a simple quasi-travelling wave for all parameter regimes? Do other such states exist? Can these states be extended into the turbulent attractor? How densely is state space packed with such solutions, and how are they connected?

The quasi-travelling wave presented could only be isolated as it was linearly stable within the edge. Even with this advantage, the bisection technique required is time consuming. To pursue these problems further more advanced techniques are required. Such techniques (primarily matrix-free Newton–Krylov solvers) have been widely applied in classical fluids to find, track and analyse steady states, travelling waves and other, more complex, classes of exact solution. Implementing these techniques within existing plasma codes is an ambitious but feasible problem which this paper motivates and begins to open the door to.