A nonlinear approach to transition in subcritical plasmas with sheared flow

A nonlinear approach to transition in subcritical plasmas with sheared flow

Chris C. T. Pringle chris.pringle@coventry.ac.uk Applied Mathematics Research Centre, Coventry University, Coventry, CV1 5FB, United Kingdom    Ben F. McMillan b.f.mcmillan@warwick.ac.uk Department of Physics, Centre for Fusion, Space and Astrophysics, University of Warwick, CV4 7AL Coventry, United Kingdom    Bogdan Teaca bogdan.teaca@coventry.ac.uk Applied Mathematics Research Centre, Coventry University, Coventry, CV1 5FB, United Kingdom

In many plasma systems, introducing a small background shear flow is enough to stabilize the system linearly. The nonlinear dynamics are much less sensitive to sheared flows than the average linear growthrates, and very small amplitude perturbations can lead to sustained turbulence. We explore the general problem of characterizing how and when the transition from near-laminar states to sustained turbulence occurs; a model of the interchange instability being used as a concrete example. These questions are fundamentally nonlinear, and the answers must go beyond the linear transient amplification of small perturbations. Two methods that account for nonlinear interactions are therefore explored here. The first method explored is edge tracking, which identifies the boundary between the basins of attraction of the laminar and turbulent states. Here, the edge is found to be structured around an exact, localized, traveling wave solution; a solution that is qualitatively similar to avalanche-like bursts seen in the turbulent regime. The second method is an application of nonlinear, non-modal stability theory which allows us to identify the smallest disturbances which can trigger turbulence (the minimal seed for the problem) and hence to quantify how stable the laminar regime is. The results obtained from these fully nonlinear methods provides confidence in the derivation of a semi-analytic approximation for the minimal seed.

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

I Introduction

In plasma physics, as in neutral fluids, configurations exist which are linearly stable but allow long-lived turbulence to develop once a large enough displacement away from the laminar state is given. These are referred to as subcritical configurations Casson et al. (2009); McMillan et al. (2009); Roach et al. (2009); van Wyk et al. (2016); Johnson and Gammie (2005); Friedman and Carter (2015) and can play an important role in the analysis and engineering of transport in plasmas. As such, methods that characterize the near-laminar and turbulent states, their interface and the threshold required for the transition between states are valuable. The methods presented here to analyze the transition from a quiescent (laminar) state to a turbulent state are not specific to any particular scenario, but will be illustrated by examining a class of subcritical systems in which a pressure-gradient driven instability is suppressed by sheared flows. This scenario has instances in non-trivial situations such as in tokamak physics, where temperature-gradient driven drift instabilities are suppressed by poloidal zonal flows or system-scale toroidal flows. However, this scenario occurs more generally in non-stationary fluids such as accretion disks.

The dynamics of sheared flows in plasmas are somewhat surprising. On one hand, even a vanishingly small shear is enough to completely linearly stabilize the system. On the other hand, the properties of turbulence tend to vary smoothly as flow shear is imposed, as long as the system is initialized far from the laminar state. The connection between linear and nonlinear properties of these systems is therefore unclear. For example, although these systems are formally linearly stable, in the sense that linear perturbations decay at late time, perturbation energies can still grow by orders of magnitude in these regimesWaltz, Dewar, and Garbet (1998); Johnson and Gammie (2005); Squire and Bhattacharjee (2014) before being damped eventually. Conversely, a highly localized small perturbation may still be sufficient to push the plasma into the regime where nonlinear effects are important and permit the system to evolve towards sustained turbulence. The threshold perturbation amplitude depends on both the linear dynamics of the system and the nonlinear processes that balance the linear decay that would otherwise take place. Practically, a subcritical system with a large nonlinear instability threshold might be able to remain in a low transport state, but one with a small threshold might in practice always exhibit turbulence. This kind of consideration is taken into account in engineering pipes for transporting fluid, and may be required for predicting the behavior of subcritical plasmas.

To fully understand the dynamics of subcritical plasmas we need to move beyond linear analysis, i.e. linear transient amplification of small perturbations, and address these issue from the perspective of the full nonlinear problem. Previously, some subcritical tokamak plasma papers have looked at identifying a boundary in parameter space above which turbulence is permanent and below which it is transient (or even non-existent) in nature.Highcock et al. (2011); van Wyk et al. (2017) These works performed parameter scans of simulations started at very large initial amplitude and run for a determined period of time. We will approach the problem in the opposite direction. For fixed parameters, we consider the initial evolution of a smoothly varied range of initial conditions to separate those which lead to turbulence from those which relaminarize, an idea touched on previously in plasma researchvan Wyk et al. (2017); Friedman and Carter (2015). This enables us to identify the basins of attraction of the laminar and turbulent states and how these basins change as the parameters are changed.

While the laminar state itself is fully identified by the linearized system, identifying the extent of its basin of attraction is a nonlinear problem. Moreover, analyzing the type and properties of the solutions located on the boundary separating these states can offer insight into the transition to turbulence, providing an understanding of the processes leading to self-sustained turbulence. The states that occur on this boundary are often considerably simpler than the fully developed turbulent state, and in the model problem, a simple traveling wave is found on the edge, allowing a direct understanding of how linear and nonlinear terms balance. Intriguingly, these traveling waves echo features of the avalanche-like bursts that permeate the turbulent regime. These methods were designed to determine the threshold amplitude and shape of the smallest disturbance that triggers turbulence (i.e. the minimal seed), but the insight they offer into the dynamics of the transition to turbulence may well be of equal importance.

In this paper we will present two fundamentally nonlinear methods of analyzing the transition to turbulence. We will do so in the context of a simple fluid interchange model (section II), however they are general and applicable to any system with two (or more) linearly stable states. The first of these methods explores the boundary separating the basins of attraction of the two states (section III). This boundary behaves in a simple but nonlinear manner and is dominated by a traveling wave solution (section IV). The second method we present looks at which disturbances are most dangerous in terms of triggering transition (section V). As a means to capture the generic parameter dependencies and provide insight into the the nonlinear growth process, a semi-analytical model of the minimal perturbation amplitude is derived (section VI). Both approaches illustrate that the critical dynamics are localized in space. This is in sharp contrast with the linear mode of maximal transient growth which is a space-filling plane wave. In section VII we draw these results together to describe an entire transition scenario.

Ii The Plasma Interchange model

The plasma model considered here, which we denote PI (for plasma interchange model), is a fluid equation with as effective gravity term giving rise to an instability: it was introduced in Beyer, Benkadda, and Garbet (2000) and previously considered in McMillan et al. (2009) in order to explain the avalanches that arise in simulations of tokamak turbulence with a background shear flow. The equations are also similar in form to those encountered in local models of convection in gravitationally confined rotating plasmas (often subject to instabilities like the magneto-rotational instability)Johnson and Gammie (2005); Squire and Bhattacharjee (2014). It is closely related to the Hasegawa-Wakatani model, differing only in the linear terms that couple the density field to the temperature evolution. Physically, the instability arising in Hasagawa-Wakatani is of resistive-drift type, whereas the model here is an interchange instability. There are a wide class of ‘local’ plasma models such as gyrofluid or gyrokinetic models, with instabilities driven by field line curvature or gravity, and nonlinearities due to plasma advection, and much of the nonlinear behavior, such as the formation of zonal flows, is generic. We therefore claim that this model will represent some of the generic features of plasma turbulence in the presence of a uniform drive from curvature or an external force such as gravity, and a background shear flow. The model is formally derived by taking a local limit of the MHD equations and parameterizing the dynamics along the field line, so that the resulting degrees of freedom describe plasma motion in terms of quantities in the 2D plane perpendicular to the field.

The model evolves the fluid vorticity (as in generic fluid equations such as Navier-Stokes), which is self-advected, and the fluid density , which is advected by the velocity field. There is a gravitational force (representing the effective force due to magnetic field curvature) which induces a force on the fluid proportional to the density. As noted by Ref. Benkadda et al. (2001), the 2D model equations are equivalent to Rayleigh-Bernard convection equations, but unlike in the Rayleigh-Bernard problem, where the fluid is confined between two plates, the domain here is 2D and infinite. Because we have added a uniformly sheared background flow, this is similar to what is sometimes referred to as the Rayleigh-Bernard-Couette (RBC) model. In the RBC model the most unstable modes are streamwise-uniform, that is, with wavevectors perpendicular both to the background flow and the gradient, but this direction is perpendicular to the 2D plane chosen: the high effective stiffness of magnetic field lines (which is not explicitly modeled) ensures that the motion is spanwise-uniform and justifies the use of a 2D model. This model is therefore a generic representation of interplay between a gradient driven instability (Rayleigh-Taylor) and shear flows, but of a slightly different form to the classic fluid models. The model may be written


where the fields and are functions of spatial coordinates and as well as time , and is the convective derivative.

To make this model applicable to microinstabilities in plasma, where only a narrow range of wavenumbers are actually unstable, we follow references Benkadda et al., 2001; McMillan et al., 2009 and consider only the modes (which represents a zonal flow or density perturbation) and (representing a drift mode). This reduces the spatial dimensionality to 1D and is appropriate where the dynamics are dominated by a narrow range of unstable modes and coupling to the zonal flow. This is motivated by the dominant role of zonal perturbations in tokamak turbulence dynamics and by post-hoc qualitative similarity of the features of this model to gyrokinetic simulationsMcMillan et al. (2009). The model captures the linear drift wave instability process, zonal flow generation via Reynolds stresses and shear flow stabilization of turbulence by zonal flows. Normally, system-scale perturbations (at low ) would dominate in the Rayleigh-Bernard model, but these are absent by construction as only a single non-zero is retained (physically, the reason that longer perpendicular wavelength mode are stable for electrostatic drift instabilities is because Landau damping dominates over charge separation). The short wavelength compared to the system scale justifies the use of a local treatment with periodic boundary conditions.

The model comprises of four 1d equations corresponding to the zonal averaged density field (), the gradient of the zonal electric potential () and the complex wave potential ( and wave density (). These quantities evolve according to the equations


where is the background shear. These are solved in the domain subject to the boundary conditions , , and . Diffusion on the zonal density was been removed on the principle that a collisionless plasma was being simulated; it was nevertheless retained for because the long wavelength zonal flows otherwise build up excessively (compared to the gyrokinetic simulations), possibly because the tertiary instability is missing. In a collisionless setting, the diffusion term on the drift modes may be thought of as arising due to a dependent growth rate.

The behavior of this system is explored in reference McMillan et al., 2009: note that scaling transformations have been used to reduce the number of free parameters. In the absence of a background flow (), the linear dynamics may be captured straightforwardly in terms of a collection of non-interacting plane waves growing due to the interchange instability, but damped due to the diffusion term, so the growth rate is maximum when the wavenumber . However, for initially plane wave modes become tilted due to the sheared background flow (in the 2D plane), and the wavenumber increases linearly with time. Given any initial plane wave perturbation, the late-time wavenumber will eventually increase until the mode is damped. The system is thus linearly unstable for , but stable for all positive values of . This situation has a familiar spatial analogue, the convective instability, where a mode propagates through a finite unstable region, and is amplified for a certain time, but later propagates into a stable region, where it is damped.

Despite the linear stability, chaotic turbulent flow is observed at finite values of . For values of the shear below the turbulence fills the domain (figure 1, left), however above this value it becomes increasingly localized. This continues to the extent that for , burst like solutions are identified traveling across the domain at near constant rate with only slow variations of small amplitude (figure 1, right). If the shear is increased past this point, turbulence becomes increasingly transient and ceases to be sustainable.

Figure 1: Left: Contours of flux, , on a plot of space against time for the PI model for . The turbulence is chaotic and fills the domain, but not without structure as individual sub-patches advect across the domain. Right: Again contours of flux but this time with . Now the turbulence is localized and travels across the domain at a constant rate leaving small fluctuations in its wake.

Iii The edge of chaos

An important consequence of the linear stability of this system is that finite amplitude disturbances are required to trigger turbulence. In general, for any given disturbance there exists a critical amplitude below which the disturbance dissolves away and above which it leads to chaos. If the critical amplitude itself is chosen then the flow continues on in an inbetween state, never relaminarising and never becoming truly turbulent. This manifold is referred to as the edge of chaos (or laminar-turbulent boundary) and has been extensively studied in a range of subcritical shear flows.Itano and Toh (2001); Skufca, Yorke, and Eckhardt (2006) The edge of chaos itself has undergone more limited investigation in plasma models.Friedman and Carter (2015)

The edge of chaos is the boundary between the basins of attraction for the laminar and turbulent states. In a phase space visualization of the system, the edge forms a hypersurface. Any initial condition ‘outside’ or ‘above’ the edge will lead to turbulence, while those below the edge will laminarise.

We can locate the edge for the PI system by choosing a form of perturbation and iteratively refining its amplitude through bisection. If turbulence is observed the amplitude is reduced, while if the disturbance decays away the amplitude is increased. As the refinement increases we are able to track to initial conditions, one on either side of the edge. Initially their evolutions mirror each other, but eventually they begin to drift apart as one relaminarises and the other becomes turbulent. The more accurately the critical amplitude is determined the longer the separation of the two trajectories can be delayed. If the precise amplitude of the edge could be found the two trajectories would never separate. Instead, the limitations of numerical accuracy mean that divergence is inevitable. In order to be able to track the edge past this, the approach is restarted with the two states just after they begin to depart the edge by interpolating between them.

In figure 2, the results of the edge tracking calculation can be seen. After an initial transient, the simulation settles down. Due to a lack of diffusion, a background distribution of develops. Superimposed on top of this the dynamics become simple, an isolated traveling wave advecting at a constant rate.

Figure 2: Main: Time evolution of energy for the edge of chaos. The edge was found by bisection and the calculation restarted when numerical precision failed. Inset: A magnified version of one of the restart points. The states on either side of the edge can be clearly seen as can the interpolated state within the edge.

Iv The edge state

All initial conditions within the PI edge, for any given value of the background shear, evolve into the same state – henceforth referred to as the edge state. This simple traveling wave is an exact solution of the underlying equations (exact in the sense that they satisfy the equations, propagating at a constant rate without variation or fluctuation), and fully nonlinear in nature. That all initial conditions within the edge evolve into this edge state is reflective of the fact that this state is linearly stable within the edge. In an absolute sense it is unstable, however its only unstable direction points out of the edge so it does not affect the edge dynamics. The scenario of an attracting edge state has also been observed in a number of other systems, albeit usually only after restricting symmetries have been imposed.Schneider et al. (2008); Duguet, Willis, and Kerswell (2008)

The traveling wave can be calculated directly. It takes the form

and so is a solution to the equations

Equations (3)-(6) are substituted into the lefthand side to form a coupled eigenvalue problem. The equations are supplemented by applying a phase condition. We solve this by making an initial guess for the travelling wave. Solutions are then sought iteratively using a Newton-Raphson approach. As this is a relatively small dynamical system (we allowed degrees of freedom), the Jacobian can be explicitly calculated. As always, successful convergence of such an approach depends on making a suitable initial guess. For this problem we are able to do so simply by taking the end state of the edge tracking.

Figure 3: The exact traveling wave solution for . The solution is clearly localized to one section of the domain. All four components contribute to the traveling wave and combine to sustain it against diffusive decay.

In figure 3 we plot the traveling wave for . It is fully localized (in the sense that the solution does not change when the length of the periodic domain is increased). Having found this one solution we can then explore how it varies as the background shear changes. Beginning from our initial solution we set to be slightly changed and then use our previous solution as the initial guess in our Newton-Raphson routine at the new shear rate. When we do so we are then able to track the solution through phase space, as seen in figure 4 (main). Increasing leads us to find there is a maximum shear rate () for which this solution exists. Beyond this point the solution turns back. Having done so, rather than then continue back to zero shear, it then snakes upwards for as far as we continued it. During this progress, the form of the traveling wave changes drastically and its amplitude rapidly increases. It remains localized, but the spatial interface between the wave and the surrounding empty space becomes abrupt as shown in figure 4 (right).

Figure 4: Main: Tracking the traveling wave solutions as the background shear is varied. The maximum shear rate it exists at is , beyond which it doubles back. Inset: Plots of for the traveling wave solution at three different points. Only part of the domain is shown. As the traveling wave is traced upwards it increases in amplitude and the spatial interface becomes more abrupt.

The traveling wave solution appears to be a natural counterpart to the propagating bursts seen in fully developed turbulent systemsCandy and Waltz (2003); McMillan et al. (2009); Görler et al. (2011); van Wyk et al. (2016): the amplitude, velocity and spatial profile of the propagating bursts in the turbulent manifold of the PI modelMcMillan et al. (2009) are quite similar to the edge state near the critical shear . The states found in the turbulent phase, however, do not appear to be exact traveling waves, and may be periodic or relative periodic orbits: we made an attempt to find an exact traveling wave in the turbulent manifold but were not able to locate any.

V The minimal seed

The traveling wave calculated in the previous section dominates the behavior within the edge. In this sense it represents the typical amplitude of disturbance required to trigger a turbulent episode. The actual amplitude required will depend on the precise form of perturbation applied. Some forms are more efficient than others and may trigger turbulence at much lower amplitudes. With this in mind it seems pertinent to ask, what is the smallest disturbance capable of triggering turbulence? It may be of practical importance to know what the minimal disturbance required in order to construct a ‘maximum safe level’ of experimental fluctuation below which the system is guaranteed to remain laminar. Identifying this disturbance is equivalent to finding the point within the edge with the lowest amplitude. This state is called the minimal seed.Pringle, Willis, and Kerswell (2012)

To answer this, we begin by considering how disturbances inside the edge behave. We have already seen that for the PI model there is a simple attractor within the edge. The nature of this state means that for any initial condition chosen from within the edge, given enough time it will evolve into this simple traveling wave. The immediate implication of this is that no matter the initial amplitude of a disturbance within the edge, its eventual amplitude will be given by that of the traveling wave solution. If we denote a state, evolving in time, as , then for all disturbances within the edge




and is the traveling wave solution.

As a consequence of this, the minimal seed – the disturbance within the edge with the smallest amplitude – will also be the initial condition within the edge which maximizes the quantity


To identify this state we begin by considering the related question: of all possible disturbances of a given initial amplitude, which exhibits the most growth, , over asymptotically large times


If the initial amplitude is fixed to be that of the minimal seed, , then clearly will be maximized by the minimal seed itself. This can be seen as for the minimal seed (by definition the minimal seed is within the edge and so must evolve into the traveling wave), while all other disturbances are below the edge and so relaminarise implying for them. Even if only large rather than infinite times are considered, the minimum seed should still maximize this quantity.

The approach is hampered by the fact that we do not a priori know . Instead, if we maximize for a succession of increasing initial amplitudes we should be identify the point at which the edge is crossed by a sudden increase in the amount of growth produced. How sudden this increase is will depend on how large we take to be.Kerswell, Pringle, and Willis (2014)

The challenge is now to identify those perturbations which generate the most growth. Small disturbances which can exhibit large growth have historically been the subject of the study of transient growth. For the linear problem, identifying the disturbances which grow the most was originally done by considering the pseudo-spectra of the operator. More recently the problem is frequently considered by formulating a functional to be maximized (see reference Schmid, 2007 and references therein). The attraction of this approach is its flexibility and the fact that it can typically be solved using minor modifications to a standard time-stepping code. Importantly, it is straightforward to include the nonlinear terms in the calculation Pringle and Kerswell (2010); Cherubini et al. (2010); Monokrousos et al. (2011).

Let us define


The first term represents the amplitude of the state at time . All the other terms are constraints enforced by Lagrange multipliers. The first of these fixes the initial amplitude while the remaining terms impose the dynamical constraints that , , and all evolve in accordance with their respective nonlinear equations. Maxima of this functional are found by seeking zeros of its first variational derivatives. We will refer to as the direct variables, and as the adjoint variables. For the variational derivatives to be zero, the following conditions must hold. For the direct variables


The adjoint variables must satisfy


The two sets of variables are linked through the optimality and compatibility conditions


(a full derivation of these equations is included in appendix A).

These zeros are found via an iterative approach. An initial guess, , is made. This can be evolved forwards in time to find , guaranteeing equations (12)-(15) are satisfied. We can then use this to ‘initialize’ the adjoint variables , satisfying equation (21). The adjoint variables can then be evolved backwards to satisfy conditions (16)-(19) (note that the sign of the diffusive term indicates that this is the correct way to integrate these equations). At this point all the derivatives of are zero except (20), corresponding to , which tells us how to update our initial guess


where is chosen to maintain the initial state’s amplitude across iterations. This procedure is then repeated until is sufficiently small.

v.1 Linear transient growth

Before calculating the minimal seed for this problem, we first take the exploratory step of considering the maximally growing disturbance of infinitesimal amplitude.

The behavior of these is captured by the linearized equations


subject to the boundary conditions and . The equations separate if we make the change of variables leading to


Progress can be made analytically by making the ansatz


This ansatz arises more naturally when the 2D equations are presented in the Lagrangian frame (moving with the fluid), where an initially plane wave solution does not change wavenumber at late time, but the geometry secularly evolves with time (Johnson and Gammie (2005) denotes these ‘shearing wave’ modes shwaves).

The above ansatz leads to an equation for ,


with solution


Maxima and minima of this function can be found by seeking zeroes of the derivative of the exponent which occur when


In order to achieve maximum growth we require a minimum at yielding the relationship , leading to a maximum at . Putting this back into our equation for reveals the maximum growth being


A similar argument leads us to conclude must always monotonically decay, as is expected from inspection of equation (25).

We confirm the linear result using the numerical method outlined previously. The procedure smoothly converges to the optimally growing disturbance (figure 5). As expected this takes the form of a single Fourier mode in and with a phase shift of between them, and nothing in or . The optimal disturbance is independent of , however the amount of growth and how long it takes does vary with the shear, as seen in figure 6. These scalings match with those derived above.

As an aside, this calculation also illustrates the linear instability of the system at as the amount of growth becomes infinite as long as .

Figure 5: The linear energy growth optimal for . From left to right: convergence of the growth rate as the algorithm iterates; decay of the amplitude of the gradient vector as the algorithm iterates; the real part of for the linear optimal.
Figure 6: The amount of growth (left) and time taken for the growth to be achieved (right) by the linear optimal as a function of . The lines are the predicted scalings of and and the crosses are the numerical results.

The linear transient dynamics of this model are similar to those in tokamak geometry at zero magnetic shearHighcock et al. (2011). They are simpler, however, than the general tokamak caseWaltz, Dewar, and Garbet (1998); Bokshi et al. (2016), where the transient dynamics take the form of repeating periods of growth followed by damping (Floquet modes), as a geometrical coupling (magnetic shear) provides a means to couple modes of different radial wavenumber. In cases where the system is linearly stable, the maximum amplitude is reached after the first growth period, so the linear optimal would be sensitive only to behavior over this period; this may also be true for the nonlinear problem.

One approach for characterizing linear transient growth for tokamak instabilities is described in Roach et al. (2009). An effective growth rate is defined by where is the time taken for the mode to grow by some predetermined factor . This can be applied to any initial condition, even if it eventually decays. The method here can also determine such effective growth rates, and choosing to optimize the initial state uniquely defines it.

v.2 Nonlinear optimals

We now turn our attention to finding the minimal seed by calculating optimally growing disturbances with finite amplitude. In theory this is a simple extension of the linear problem, but in practice it suffers from a number of technical problems that need to be understood.

The first is that although the adjoint equations are in fact linear in the adjoint variables, they do require that the forward variables are stored for all times as they also appear in the equations. For large systems (or long time runs) this can become memory intensive. This can be worked around by introducing checkpointing where the forward state is only saved a smaller number of discrete times (say , , ). During the backwards run the forward variables can be reconstructed in a piecemeal fashion with each section of the direct variables recovered when it is needed by integrating forward again.

Choosing the target time is also a challenge. For the linear problem we optimized for to maximize the growth. Although that can also be done here,Rabin, Caulfield, and Kerswell (2012) it is easier – and more consistent with our argument – to simply fix the target time as something large. Although it is tempting to set to be as large as can be feasibly managed, this leads to complications and in particular convergence can be very slow (in terms of the number of iterations required) and perhaps not even possible. This is because as most states will evolve towards zero after an initial transient, can be susceptible to numerical noise and the gradient vector becomes a poor approximation of the true one. For our purposes we settled on taking which is large enough for the dynamics to develop but not for damaging levels of viscous decay to occur.

Iterative approaches are not guaranteed to converge. For small initial amplitudes, convergence should be straight forward and the result is expected to be the same as the linear one (or at most a small modification to it). Once larger amplitudes are encountered, however, it seems less clear what to expect. In general, if the initial states being considered evolve into a chaotic regime then we don’t expect convergence at all. Small changes to the initial state lead to large changes in the energy growth and the problem effectively becomes non-smooth. If we consider a system for which the edge separates chaos from smooth deterministic behavior, then we expect subcritical amplitudes (below the edge) runs to converge and supercritical amplitudes (above the edge) to fail to converge. If the edge itself is chaotic then critical amplitude runs will also not converge.

This gives a simple criterion for identifying the critical amplitude by running the algorithm at a sequence of amplitudes and identifying where convergence ceases to be possible. This abrupt cut off becomes smoothed when a finite optimization time is taken. As the critical amplitude is approached, chaos takes longer and longer to develop until this time frame is larger than the chosen optimisation time. The approach is further hindered when either the non-trivial state or the edge is not chaotic. We have already seen the latter is the case for the PI model, and for some larger values of the former is also true. If the non-trivial state is deterministic then convergence is possible above the edge as smoothness is not lost. This is not a problem as the amount of growth should be distinct from that seen below or within the edge. If the edge is deterministic (i.e. there is a simple traveling wave within the edge as here) then for finite times disturbances close to, but on either side of, the edge will both evolve to something approximating the edge state. Although given enough time these would diverge away from the edge state, if they are arbitrarily close to the edge then the amount of time required is arbitrarily large. The amount of growth on either side of the edge after more moderate times will be comparable and the only way to determine whether the edge has been crossed is to use the optimal as the initial condition for a longer time simulation.

v.3 Critical amplitudes

In figure 7 we show the results of the algorithm for for two different initial amplitudes – and . Although the amplitudes are very similar, fundamentally different results are obtained. The lower amplitude smoothly converges. The larger amplitude starts to follow the same path before diverging. By checking how each of the states the algorithm iterates through evolve in time, we see that this divergence corresponds to the point at which we begin to see turbulence in the system. The last few states that are iterated through at all lead to turbulent episodes and so are turbulent seeds above the edge. All of the states found for are below the edge, including the final converged optimal. We conclude that the critical amplitude, , is sandwiched in between these two values. In this way we bound the minimum amplitude required to trigger turbulence.

Figure 7: The nonlinear energy growth optimal at . From left to right: convergence of the growth rate as the algorithm iterates; decay of the amplitude of the gradient vector as the algorithm iterates; the amplitude evolution of the nonlinear optimal and the first turbulent seed found by the algorithm. In all cases the purple line corresponds to and the green line corresponds to .

Both the nonlinear optimal and the first turbulent seed (the first initial condition found by the algorithm that initiates turbulence) are localized and show similar structure (figure 8). The localization is to be expected. This allows the disturbance to have local amplitudes large enough to nonlinearly self-interact (a necessity for turbulence) while simultaneously keeping the total amplitude small. Despite the minimum seed being localized in space, it still initiates space-filling turbulence (figure 9). After an initial transient (),the flow approaches the traveling wave solution () before growing into the fully turbulent state (). This turbulence spreads in both directions to fill the entire domain. The amount of time this takes is a function of how close to the critical amplitude the initial amplitude is. The initial growth towards traveling wave takes an amount of time approximated by the time-scale of the linear transient growth. As it approaches the traveling wave, how close it comes to that solution (and hence how long the next stages take) depend on the initial amplitude as well as the leading eigenvalues of the traveling wave solution.

Figure 8: The nonlinear energy growth optimal (, purple) and the turbulent seed (, green) at . The two states are both strongly localized and are similar in form, the principal difference being in the fields.
Figure 9: Evolution of the first turbulent seed for . The seed is localized and its initial evolution reflects this as it approaches the localized traveling wave solution within the edge. From here it then evolves into higher energy turbulence which spreads to fill the domain. The contours are based upon the flux, , but at each time it is normalized by the instantaneous maximum value of the flux at that time: . This normalization is chosen as the turbulent seed grows by several orders of magnitude during it evolution. Contours are , , and of the instantaneous maximum flux.

The picture changes a little as we move to higher values of the background shear. At , the observed turbulence is much smoother and more slowly varying. This is also true of the behavior within the edge. Because of this, initial conditions near the edge spend much longer near to the edge and our iterative approach converges above the edge even for these moderate choices of (figure 10). Due to this care must be taken to identify the point at which the edge is crossed by examining individual runs rather than simply finding the point at which convergence is no longer possible.

Figure 10: The nonlinear energy growth optimal at . From left to right: convergence of the growth rate as the algorithm iterates; decay of the amplitude of the gradient vector as the algorithm iterates; the amplitude evolution of the nonlinear optimal and the first turbulent seed found by the algorithm. In all cases the purple line corresponds to and the green line corresponds to .

Because we are able to converge so close to the edge, we are able to much more accurately pin down the minimal seed. The fact the energy growth optimals on either side of the edge are so similar in form (figure 11) strongly indicates that the minimal seed, which is sandwiched between them, must also be very similar.

Figure 11: The nonlinear energy growth optimal (, purple) and the turbulent seed (, green) at . The two states are barely distinguishable and are both strongly localized.

By repeating this approach for a range of values of , we are able to see how the minimum amplitude of the edge scales with the background shear. In figure 12 we plot the amplitudes of the minimum seed and the traveling solution as functions of . We are unable to continue the minimum seed calculation to values of quite as high as the traveling wave goes because turbulence starts to become difficult to sustain here. The method adopted can only be expected to work when turbulence is sustained (otherwise after long times all disturbances must return to the laminar) and there is a clear differentiation between the amplitude of the observed turbulence and that of the edge. If these amplitudes become too similar then the optimization algorithm has no reason to find seeds that lead to turbulence.

Figure 12: The amplitudes of the minimum seed (solid) and of the traveling wave solutions (long dashes) as functions of the background shear. Also included is the predicted critical amplitude from the semi-analytical analysis (short dashes) found in section VI. Unsurprisingly the minimal seed triggers turbulence at amplitudes significantly below the edge state. The semi-analytic approach agrees well with the minimal seed, especially for small values of where the approach is better justified.

One might attempt to trigger turbulence with either the linear optimal or the most unstable linear mode, but this would fail as both these disturbances are plane waves which must eventually decay, even if they transiently attain amplitudes much larger than typical turbulence levels (although small amounts of noise in the initial perturbation may be enough to allow a transition). This is an indication of the potential for an approach based simply on transient linear growth to be misleading in, for example, identifying the most dangerous perturbations.

Note that the nonlinear stability of plane wave solutions is a generic feature of periodic shearing-box systems with convective nonlinearitiesJohnson and Gammie (2005); Squire and Bhattacharjee (2014). Even for systems for which a plane-wave linear optimal can trigger turbulence, the critical amplitude will scale with the system size, unlike that of the nonlinear optimal, which usually becomes independent of system size. Only the nonlinear analysis captures the fact that turbulence triggering is a local process.

Vi Semi-analytic approximation

The minimal seed calculation performed in the previous section provides a means to calculate the critical amplitude but does not directly explain which mechanisms are responsible for allowing a small initial perturbation to evolve to a nonlinearly saturated state. In this section we present a semi-analytic approach to approximating the minimal seed, and hence provide a simple closed form estimate for the critical amplitude, thus developing some insight into the processes occurring during the initial transient.

The approach revolves around the assumption that for turbulence to be sustained, the background shear needs to be ameliorated (this is justified partly by examining the time evolution of in earlier results). This occurs when is of a similar size to . This is exactly the point where the linear and nonlinear terms in the time evolution of and become comparable. We assume, for the sake of obtaining a rough estimate, that up to this point we can use the linear time evolution for and as in section V.1.

We consider the evolution of an initial disturbance, . The wave density is then given by (in our previous notation this gives ), while the remaining two components, and are taken to be zero. After applying the Fourier transform


we know this new field evolves as


As the disturbance evolves, its effective wave number changes linearly with time. This can be scaled out by introducing the new variable in which case


We are interested in the evolution of , governed by


This will be dominated by the behavior when is at its largest. As previously seen, the amplification factor in equation (33) has a maximum at , . Close to this we can approximate the behavior of by


plus higher order terms which we have dropped, along with the ’s, for clarity.

We make the ansatz


corresponding to the plane wave linear optimal but modulated by a Gaussian and so localized in space. The choice of is made here, both because it simplifies the algebra, and also because this gives the most localization of the initial condition in real space while still exciting mostly modes with near-maximal linear amplification. With this choice the behavior close to its maximum amplitude is


which in real space gives


Inserting this into equation (34) we find


We estimate the maximum value that achieves by neglecting the diffusive term on the RHS and first integrating with respect to . Setting will give us the most shear, which is in turn evaluated by differentiating with respect to to get


For a broader wavepacket with the same norm, the derivatives would have lead to a smaller prefactor; an explicit optimization of would balance this against better localization in wavenumber space for maximum transient growth. balances the background shear, , when


In figure (12) we include this predicted amplitude for transition. The agreement between this prediction and the amplitude of the minimum seed is good, especially at low which fits with our assumptions. As a means of estimating the critical amplitude, this approach captures several key aspects of the dynamics. Firstly, the scaling of the true critical amplitude seems to be dominated by the transient growth factor, indicating that linear transient growth is the main mechanism able to amplify the initial seed perturbations. Secondly, the spatial localization of the initial condition (the full-width at half-maximum of ) and transition state is well predicted (to within a factor of 2). Thirdly, the peak initial wavenumber of the minimal seed is similar to that of the semi-analytical initialization.

It has been shown that there is a correlation between the region where turbulence may be sustained and and the level of linear transient amplification for some tokamak casesvan Wyk et al. (2017), and this kind of semianalytical treatment might be able characterise the role of nonlinear processes in the transition from turbulence to the laminar regime (other approaches are also promisingConnaughton, Nazarenko, and Quinn (2015)). Note, however, that the properties of the transition to turbulence and those of the saturated turbulent state are not directly related: for example, the amplitude where the nonlinear transition occurs increases with , whereas the average turbulence amplitude decreases with .

The major benefit of this semi-analytical approach is that it allows a quantitative understanding of the mechanisms allowing small perturbations to evolve towards a nonlinearly self-sustaining state, and provides scaling estimates. The full minimal seed calculation is a necessary preliminary step, from which the active mechanisms can be determined. The assumption, for example, that the initial evolution of occurs due to the shape of the amplified wave-packet was based on inspecting the time-evolution of the minimal seed: this is rather different to a common assumption elsewhere that the modulational instability is key to zonal-flow growth during nonlinear saturationDiamond et al. (2005).

Vii Conclusions and discussions

We have presented two methods useful in the analysis of the transition to turbulence in linearly stable systems. Such subcritical configurations are exhibited by many plasma systems of interest, e.g. the suppression of temperature-gradient driven drift instabilities by poloidal zonal flows in tokamaks. For these subcritical configurations, linear approaches struggle to give meaningful results, while the nonlinear approaches presented here offer fresh insight.

The edge of chaos is the boundary between the basins of attraction of the laminar state and the turbulent one. By tracking trajectories confined to this boundary, we found that it is dominated by a traveling wave solution which is an attractor within the edge. This traveling wave is localized in space, even for shear rates at which turbulence is space filling. This is indicative of fact that one can initiate global turbulence with local disturbances. Similar behavior is seen in classical shear flowsDuguet, Schlatter, and Henningson (2009); Avila et al. (2013). By tracking the traveling wave through phase space we are able to quickly track the typical amplitude of the edge. Further, the maximum shear rate for which turbulence is sustained seems to correspond to the the maximum shear rate at which the traveling wave exists. It should be noted that the traveling wave ceasing to exist does not have to exactly correspond to the point at which turbulence ceases to exist,Avila et al. (2011) but it is expected to be related. At this point, something fundamental has to change in the nature of the edge for this system and hence for the system as a whole.

If the edge state represents the typical turbulence inducing disturbance, the minimal seed is the most dangerous disturbance – the smallest disturbance capable of triggering turbulence. As with the edge state, this seed is localized. This should be expected as it allows the disturbance to be of low total amplitude while still being locally large enough to exhibit nonlinear behavior.

The play off between local and global amplitude is explored in the semi-analytic approach we put forward which qualitatively captures the key dynamics. The semi-analytic approximation is shown to offer a good agreement with the nonlinearly computed minimal seed amplitude. It may offer a computationally efficient approach in systems where performing a full nonlinear optimization is not feasible (an important consideration in kinetic models). However, the principal outcome is a quantitative understanding of the mechanisms that allows small perturbations to evolve towards a nonlinearly self-sustaining state.

The two approaches (edge tracking and minimal seed calculation) link to each other. In figure 13 we plot the evolution in phase space of the two states sandwiching the minimal seed at , previously shown in figure 11. As expected both of them initially display closely matching evolutions as they track their way along the edge. Eventually they start to diverge, one relaminarising while the other exhibits turbulent behavior. The trajectories depart one another as they reach the traveling wave solution, marked with a cross.

Figure 13: Phase space plot of flux () against amplitude for the transition to turbulence at . The orange line is the nonlinear energy growth optimal for while the purple line is the optimal for , shown in figure 11. Both initially evolve in a very similar manner, approaching the traveling wave for this shear, marked ‘’. As they approach the edge state the two trajectories diverge to the laminar state and turbulent state respectively.

The future task is to apply these methods to the full plasma problem. Edge tracking is a straight forward tool that can be readily implemented within already existing codes. Although it seems unlikely to lead to the discovery of attracting solutions embedded within the edge, unstable solutions may be identified. Close approaches to these exact solutions within the edge can be found by seeking moments of near-recurrence.Duguet, Willis, and Kerswell (2008) Newton-Krylov methods based on well-established routines such as GMRESSaad and Schultz (1986) reduce the computational intensity of solution tracking. The probability of successful convergence can be enhanced by using so-called globally convergent methods.Dennis and Schnabel (1996) Calculating minimal seeds for large problems is also computationally intensive. A possible way to reduce the overhead would be to instead consider smaller domains in which turbulence may even be space filling. In the case of shear flows the minimal seeds for such reduced problems appear to closely mirror those of the full, non-periodic problem.Pringle, Willis, and Kerswell (2015) For the PI model considered, the difference in amplitudes between the edge state and the minimal seed is not excessive. This is somewhat atypicalPringle and Kerswell (2010); Rabin, Caulfield, and Kerswell (2012); Duguet et al. (2013) and it seems likely that for the full plasma problem this energy gap will increase drastically.

Despite these computational hurdles, the methods presented here seem ripe for applying to tokamak problems.McMillan, Pringle, and Teaca () The limitations of linear approaches are becoming increasingly clear and so fully nonlinear methods allow for a fresh take. Ultimately, stabilizing subcritical plasma flow requires a means of quantifying the nonlinear stability of a system. Both edge tracking and minimal seeds offer a meaningful way of doing so. Endeavors to incorporate these into classical control have already shown promise.Kawahara (2005); Rabin, Caulfield, and Kerswell (2014)

Viii Acknowledgments

We thank Zaid Iqbal for useful discussions. C. C. T. Pringle is partially supported by EPSRC grant no. EP/P021352/1. B. McMillan is partially supported by EPSRC grant no. EP/N035178/1. B. Teaca is partially supported by EPSRC grant No. EP/P02064X/1.

Appendix A Derivation of the adjoint equations

In section V we defined the functional


which we sought to maximize over all possible choices of initial condition. In order to find these optimals, we used the variational derivatives found here.


The majority of these terms are in an appropriate from, but the terms in the form

require manipulation. The time-derivative terms are all similar and in general can be dealt with as


The diffusive term yields


The last step is achieved by applying the boundary conditions




The terms that incorporate the background shear give us




having again applied boundary conditions. We now consider


The final term is