Abstract

In tokamak plasmas, sheared flows perpendicular to the driving temperature gradients can strongly stabilize linear modes. While the system is linearly stable, regimes with persistent nonlinear turbulence may develop, i.e. the system is subcritical. A perturbation with small but finite amplitude may be sufficient to push the plasma into a regime where nonlinear effects are dominant and thus allow sustained turbulence. The minimum threshold for nonlinear instability to be triggered provides a criterion for assessing whether a tokamak is likely to stay in the quiescent (laminar) regime. At the critical amplitude, instead of transitioning to the turbulent regime or decaying to a laminar state, the trajectory will map out the edge of chaos. Surprisingly, a quasi-traveling-wave solution is found as an attractor on this edge manifold. This simple advecting solution is qualitatively similar to, but simpler than, the avalanche-like bursts seen in turbulent simulations and provides an insight into how turbulence is sustained in subcritical plasma systems. For large flow shearing rate, the system is only convectively unstable, and it will eventually return to a laminar state at a fixed spatial location.

pacs:

52.35.Ra, 52.30.-q, 47.27.Cn

Introduction.— The strong effect of sheared flows on linear plasma instabilities results in a broad range of subcritical configurations, which are linearly stable but allow long-lived turbulence to develop given a large enough displacement from equilibrium. Such subcritical configurations in plasmas Johnson and Gammie (2005); Friedman and Carter (2015) and tokamaks more specifically Casson et al. (2009); McMillan et al. (2009); Roach et al. (2009); van Wyk et al. (2016) have recently come under extensive study. Computationally, the late-time properties of the turbulent state in a subcritical plasma can be determined by giving the plasma a sufficiently large initial kick Casson et al. (2009), but whether or not an experimental plasma would enter this regime depends both on the threshold for nonlinear instability, and details of the experimental time-history.

In this letter, we investigate the threshold in phase space (not parameter space) between the quiescent (laminar) and turbulent state, the edge of chaosItano and Toh (2001); Skufca et al.(2006); Pringle et al.(2017), and the dynamics of the plasma on that threshold, to answer practical questions about how large an initial perturbation is required to induce long-lived turbulence in a tokamak configuration. It also provides insight into how nonlinear and linear terms combine to allow quasi-steady states in plasma turbulence. We find a simple propagating state on the edge of chaos, and this allows us to provide a relatively simple picture of how nonlinear effects can sustain the dynamics in a linearly stable regime. This propagating edge-state mirrors many of the features of the propagating bursts/avalanchesCandy and Waltz (2003); McMillan et al. (2009) seen in the turbulent regime.

The system analyzed.— This strongly magnetized plasma system, an idealised tokamak, is described via a gyrokinetic (GK) framework Hahm (1988). The GK equations describes particle motion in magnetized plasmas in a self-consistent electromagnetic field by evolving the gyrotropic particle distribution function f(x,y,θ,μ,v∥), where x,y and θ are spatial coordinates and μ and v|| are velocity space coordinates. We use the GKW codePeeters et al. (2009) to evolve the local gyrokinetic equations for the standard CYCLONE tokamak parametersDimits et al. (2000). In addition to the self-organization of the electrostatic potential ϕ and the emergence of zonal flows Diamond et al. (2005), a background poloidal flow shear of amplitude ωs is imposed. For zero shear, the system is linearly unstable to ion temperature gradient instabilities, resulting in a strong radial heat flux.

Finding the edge of chaos.— The black traces of Fig. 1 show the evolution of the heat flux in the system for initializations of varying amplitude versus time (in units of transit frequency t0). The linear system is stable despite periods where the flux increases in time, and simulation with sufficiently small amplitude initializations decay. Given a large enough initialization, however, sustained turbulence is triggered.

The amplitude of the initial perturbation was systematically varied, by a bisection technique Pringle et al. (2017), 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 Fig. 1 the near-stationary flux (log10[flux] ≈−1) of the edge state indicates that the edge dynamics are considerably simpler than the turbulent dynamics. The ‘steps’ that appear in decaying simulations (log10[flux] ≈−4) are associated with a time dependent (‘Floquet mode’) behaviourWaltz et al. (1998); these dynamics are too slow for these parameters to play a role in the transition to turbulence.

Figure 1: Heat flux versus time for simulations with ωs=0.12 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 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. Within the edge, however, we find an attractor, the edge state.
To analyze 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 ρ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 will fewer degrees of freedom. The properties of the edge state are qualitatively and quantitatively very similar for standard and narrow simulations.

Figure 2: Mean ϕ2 (averaged over y) at the midplane versus time and position in the
travelling wave frame x−vt for the edge state in narrow simulations with ωs=0.12. A periodicity over 3.2t0 is visible.

For narrow simulations, the edge state is found to be very close to a traveling wave. Detailed inspection (Fig. 2), however, reveals a small time oscillation, with period (3.2t0) 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 traveling 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 RMS variation from exact periodicity (when sampled once every 3.2t0) is found to be 0.5% over a period 64t0. 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 5: Figure 6: Figure 7: Figure 8: Mean of squared non-zonal potential ϕ2 (averaged over y) at the mid-plane versus time and x for ωs=0.15 (a) and 0.16 (b), and for an edge state with ωs=0.12. In (a), long-lived turbulence is seen in a slowly expanding region centered around the excitation front. In (b), turbulence is excited transiently over a period of 100t0, remaining localized near the traveling excitation front but then decays. The edge state (c) is much simpler and smoother.

The quasi-traveling wave solution in both narrow and standard simulations (Fig. 3) consists of a tilted finite ky traveling wave mode, fed by the gradient-drive, that produces a traveling zonal shear flow ahead (leftwards) of the pulse, that opposes the background shear flow (Fig. 4). The traveling perturbation strengthens the temperature gradient ahead of the pulse, and weakens the gradient behind, as expected from the energy transport equation, given the localized heat flux associated with the burst (the change in gradient in fig. 4 is of comparable size to the background gradient). Those two mechanisms would be compatible with a traveling wave in either directionMcMillan et al. (2009); Pringle et al. (2017), but when the nonlinearity in the simulation is turned off, the ky≠0 mode amplitude continues to propagate (not shown) in the same direction for 10a/cs due to the group velocity, which depends on the mean kx value and thus the sign of ωsMcMillan et al. (2009). Note that narrow simulations do not permit a vortex pairHorton and Hasegawa (1994) advection mechanism, where a spatially localised ‘blob’ self-advects across the domain, as they are dominated by one ky mode (unlike y-localised features seen elsewherevan Wyk et al. (2017)). Time snapshots (fig 3) of the mid-plane potential for narrow and standard simulations show similar tilting and localization 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.

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 (Fig. 8), 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 McMillan et al. (2009); Candy and Waltz (2003); van Wyk et al. (2016). Propagating phenomena have been frequently observed in tokamak turbulence simulations, especially when a background shear flow is imposed. Although these features are also present in the zero-shearing-rate limit McMillan et al. (2009), we see very clear propagating features for large shears, and these become more isolated as the shearing rate approaches the critical value. The propagating edge state has some features in common with the late-time nonlinear state (fig. 8), such as similar propagation velocities and spatial scales. Both are associated with a moving turbulence front supporting a moving zonal electric field that destabilizes the system in front of the turbulence front, so can both be seen as a traveling excitation wave. As in neutral fluid simulation, if the shear is increased beyond ωs=0.15 we observe in the tokamak case relatively long-lived turbulence that unpredictably decays to the quiescent state. It is clear from figures such as fig. 8 that for large shear (ωs≳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 and Champagne (1973), where these puffs travel in the direction of fluid flow, the bursts here travel is 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. However, if an ‘open’ boundary condition were applied, the system would become quiescent after the puff traveled to the boundary, regardless of initial condition. This sensitivity to boundary conditions is a consequence of the strong non-locality in these subcritical systems, which possess two (meta-)stable states.

Scaling of the transition threshold with background flow shear. — The scaling of transition threshold with shearing rate (Fig. 11) 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 exact value of the transition threshold depends on the initialization: we varied the parameters of a tilted wavepacket-like initialization to find the ‘most dangerous’ state with a minimal nonlinear instability threshold (a complete minimization would explore all possible initial states). 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. This should be contrasted with the typical situation in neutral fluid experiments, where the transition may involve several stages of linear growth chained together as a result of nonlinear effectsPringle et al. (2012).

A traveling-wave type edge state is found for all shear rates for the narrow simulations and for ωs>0.04 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.06 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.

{subfigure}

0.23
{subfigure}0.23

Figure 9: Figure 10: Figure 11: (a) Logarithm of maximum value of squared non-zonal potential maxx<ϕ(x,y)2>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.06 these results for both simulation types are similar.

Conclusions and discussions.—
The behavior of the edge of chaos is qualitatively similar to simple plasma-interchange (PI) modelPringle et al. (2017), strengthening the thesisMcMillan et al. (2009) that qualitative features in the dynamics are the consequence of fluid-like behavior, 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 state contains a quasi-traveling wave attractor (which is also seen in the PI model). The increasing amplitude of the edge state to levels comparable to the turbulent fluctuation level (fig 11) 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 traveling-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., 2017): 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 traveling wave destabilizes 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 ky modes. These propagation mechanisms appear to carry over to the avalanche/burst features in the fully turbulent state McMillan et al. (2009). 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 (fig. 11) for which a local perturbation 10% as large as the typical turbulence level is required to intitiate turbulence.

Acknowledgements.

Acknowledgments.—
B. McMillan is partially supported by EPSRC grant no. EP/N035178/1.
C. Pringle is partially supported by EPSRC grant No. EP/P021352/1.
B. Teaca is partially supported by EPSRC grant No. EP/P02064X/1.
Simulations were performed with the support of Eurofusion and MARCONI-Fusion. We acknowledge the authors of the GKW code for developing and distributing this software tool.