Transport Bifurcation Induced by Sheared Toroidal Flow in Tokamak Plasmas
Abstract
Firstprinciples numerical simulations are used to describe a transport bifurcation in a differentially rotating tokamak plasma. Such a bifurcation is more probable in a region of zero magnetic shear than one of finite magnetic shear because in the former case the component of the sheared toroidal flow that is perpendicular to the magnetic field has the strongest suppressing effect on the turbulence. In the zeromagneticshear regime, there are no growing linear eigenmodes at any finite value of flow shear. However, subcritical turbulence can be sustained, owing to the existence of modes, driven by the ion temperature gradient and the parallel velocity gradient, which grow transiently. Nonetheless, in a parameter space containing a wide range of temperature gradients and velocity shears, there is a sizeable window where all turbulence is suppressed. Combined with the relatively low transport of momentum by collisional (neoclassical) mechanisms, this produces the conditions for a bifurcation from low to high temperature and velocity gradients. A parametric model is constructed which accurately describes the combined effect of the temperature gradient and the flow gradient over a wide range of their values. Using this parametric model, it is shown that in this reducedtransport state, heat is transported almost neoclassically, while momentum transport is dominated by subcritical PVG turbulence. It is further shown that for any given input of torque, there is an optimum input of heat which maximises the temperature gradient. The parametric model describes both the behaviour of the subcritical turbulence (which cannot be modelled by the quasilinear methods used in current transport codes) and the complicated effect of the flow shear on the transport stiffness. It may prove useful for transport modelling of tokamaks with sheared flows.
pacs:
I Introduction
A major obstacle to the development of a working magnetic confinement fusion device is the transport of heat by pressuregradient driven turbulence (1). Experimental results indicate that a shear in the equilibrium plasma flow can reduce, or even eliminate, this turbulence (2); (1); (3). Several numerical studies (gyrofluid (4), quasilinear (5), and gyrokinetic (6); (7); (8); (9); (10); (11)) have reproduced this effect and established that it results from the shear in the flow velocity component perpendicular to the equilibrium magnetic field. It has also been noted that a shear in the flow velocity component parallel to the field, which leads to a linear instability (12); (13); (14), can have the opposite effect and increase the turbulence amplitudes (7); (9); (10); (11). Some numerical studies have concluded that, in the regimes considered, when the effect of parallel flow shear is included the turbulence is never fully suppressed (6); (7). This contrasts with experimental results where almost neoclassical, i.e., collisional (15), levels of heat transport have been observed (2); (3).
The ratio between the destabilising parallel and stabilising perpendicular components of the flow shear is set by the angle of the flow with respect to the magnetic field. In general, friction between trapped and passing particles and magnetic pumping keep the equilibrium flow nearly toroidal, with a large component of destabilising parallel flow shear (16); (17); (18). Nonetheless, two recent papers (10); (11) which investigated values of the flow shear several times larger than the linear growth rate of the ion temperature gradient (ITG) instability (one of the principal drivers of turbulence in fusion devices) found that for relatively low temperature gradients there is a range of flow shear values where a purely toroidal flow can completely quench the turbulence.
The question remains how large shears in the equilibrium toroidal flow can be achieved, given finite available sources of angular momentum. In Refs. (10); (11), the turbulent flux of angular momentum was calculated and found to be large at moderate flow shears; if it is too large, increasing the input of angular momentum will not generate strong enough shear in the flow. Ref. (10) considered the standard Cyclone Base Case (19), with the normalised inverse magnetic gradient scale length and a range of values of the temperature gradient and the velocity shear. In this paper, which expands the results reported in (11), we consider a similar regime to that described in (10), with the difference that here the magnetic shear . This choice is motivated by the experimental studies which have found that internal transport barriers (local regions of very steep flow and temperature gradients) tend to form in regions of zero or low magnetic shear (3); (20), and by previous numerical results which show a lower ion thermal diffusivity at zero magnetic shear (21). The results reported below show that the range of flow velocity gradients and ion temperature gradients where the turbulence is suppressed is much larger than in the case of finite magnetic shear and, therefore, that a transition in both the flow and temperature gradients from low to high values is more readily achieved (such a transition is in principle also possible with , but in a much smaller region of parameter space (22)). A positive feedback between the increase in the flow gradient and the suppression of the turbulence provides the mechanism for a jump from low to high flow gradients, with a simultaneous jump in the temperature gradient. We show that this transition results in a state where the transport of heat is nearly neoclassical, whilst the transport of momentum remains largely turbulent.
In a fusion device, the quantities that can be controlled are the input of heat and toroidal angular momentum, which in steady state are proportional to their outgoing fluxes. By varying these quantities, it is possible to vary the equilibrium gradients. In contrast, in the local gyrokinetic simulations used in this study the input parameters are the gradients, and the fluxes are calculated from the output. In order to demonstrate the existence of a bifurcation in the gradients at fixed values of the input fluxes, we use local gyrokinetic simulations to map out the dependence of the turbulent fluxes on the temperature and flow gradients, and then we add the neoclassical fluxes and invert this numerically determined dependence to find the gradients as functions of the total fluxes. We first do this by straightforward interpolation in the parameter space, then propose a simple parameterisation of the fluxes that fits the data and can be used to predict in which parameter regimes bifurcations can occur.
The rest of this paper is organised as follows. In Section II, we present the gyrokinetic system of equations used in our simulations and describe the numerical setup. In Section III, we report a numerical scan in two parameters: the flow shear and the ion temperature gradient, calculating the turbulent heat and momentum fluxes over wide intervals of these parameters. In Section IV, we describe the subcritical turbulence, driven by the ion temperature gradient (ITG) and the parallel velocity gradient (PVG), that is responsible for the turbulent fluxes. In Section V, we interpolate our results and determine how the temperature and flow gradients depend on the heat and momentum fluxes: a transport bifurcation is obtained as a result. In Section VI, we construct a parametric model of the transport. This allows us, in Section VII, to study the effect of the transport bifurcation on the temperature gradient and investigate the range of heat and momentum flux values for which we expect transitions to reduced transport to occur. In Section VIII we summarise and discuss the results.
Ii Equations and Numerics
ii.1 Gyrokinetics with Velocity Shear
We use the gyrokinetic approximation (23) to calculate the way in which the equilibrium gradients of the temperature and flow affect the transport of heat and momentum by turbulent fluctuations. For a detailed description of highflow gyrokinetics, the reader is referred to (18); (24); (25) and references therein. Here a brief summary is presented.
We consider a toroidal plasma with an axisymmetric equilibrium magnetic field .
The field lines lie on closed flux surfaces,
which can be labelled by the poloidal magnetic flux
contained within each surface.
The magnetic field can be written as
,
where is the toroidal angle,
is the radial coordinate measured from the central axis of the torus (the major radius),
and is the toroidal magnetic field.
Owing to the fast motion of the particles along the magnetic field lines,
equilibrium quantities such as the density
We allow an equilibrium plasma flow , of the same order as the ion thermal velocity
(1) 
where is the temperature of the ions and is their mass
The state of each species of charged particles can be described by its distribution function, , whose evolution is given by the VlasovLandau equation:
(2) 
where is the velocity coordinate and the collision operator. In order to describe the turbulent fluctuations which give rise to the transport, the VlasovLandau equation can be expanded by splitting into an equilibrium part and a perturbed part . The latter is further split into averaged and fluctuating parts; so
(3) 
where the angle brackets denote a spatial and temporal average over fluctuations. It is assumed that the perturbations are small and that the gyrokinetic ordering (23); (18) holds, so that:
(4) 
where
and are the typical parallel and perpendicular wavenumbers of the fluctuations,
is the growth rate of the fluctuations,
and are the gyrofrequency and Larmor radius of the ions, respectively,
and is scale length of the variation of .
It is then possible to show
(5) 
where is the charge of species and is the perturbed electrostatic potential. The nonBoltzmann part of the distribution function, , is the gyrocentre distribution, independent of the gyrophase and a function of time , the gyrocentre position , the magnetic moment and the energy of the particle (18). The gyrocentre position , where is the position coordinate, is the gyrofrequency of species , the unit vector is the direction of the equilibrium magnetic field, and is its magnitude.
The turbulence can be affected both by the gradient of the flow and by its magnitude,
the latter through the Coriolis and centrifugal drifts.
Here we neglect the effects of the magnitude of the flow
(the consequences of this simplification are discussed briefly in Section VIII)
but allow the gradient of the flow to be of order the growth rate of the fluctuations,
,
as is necessary for the flow gradient to affect the turbulence nonnegligibly (4)
Under these combined assumptions, it can be shown that evolves according to the following gyrokinetic equation:
(7) 
where and are the equilibrium density and temperature of species , respectively, is a gyroaverage at constant , , is the magnetic drift velocity, is the parallel (peculiar) velocity, is the sign of , is the mass of species , and is the speed of light. We have further simplified the problem by assuming purely electrostatic perturbations (), so the system is closed by the quasineutrality condition:
(8) 
where is a gyroaverage at constant . The electrons are treated through a modified Boltzmann response (26):
(9) 
where the overbar denotes averaging over the flux surface. Thus, Eq. (7) is only solved for ions: .
Knowing , we can calculate the turbulent heat and angular momentum fluxes associated with a given minor radius and given values of the gradients,
(10) 
(11) 
ii.2 Numerical Setup
The nonlinear gyrokinetic equation (7) is solved using the gyrokinetic code GS2 (27); (28). Its solution depends on (among other quantities) the equilibrium gradients: , and (the latter of which is the gradient the of parallel component of the background flow). To simplify the treatment of these gradients, we consider a small region of the plasma known as a flux tube (29), which encompasses a set of magnetic field lines, extending to include several turbulence decorrelation lengths in each direction. This allows equilibrium quantities to be expanded around a reference flux surface near the centre of the domain, labelled by . Thus we may write:
(12)  
(13)  
(14) 
where is the equilibrium density scale length, the equilibrium temperature scale length. The perpendicular velocity shear and the magnetic safety factor are:
(15) 
respectively, where is the major radius of the magnetic axis. For a given equilibrium magnetic field, the solution of a flux tube simulation and, therefore, the values of and , is a function of the three parameters , and . Defining the local radial coordinate as
(16) 
and the second perpendicular coordinate as
(17) 
where is the poloidal angle
(18) 
The gyrokinetic equation (7) may now be written
(19) 
In this investigation, we vary the perpendicular flow shear , and . Other parameters are kept fixed:
(20) 
The density gradient, the magnetic safety factor and the aspect ratio are chosen to conform to the Cyclone Base Case, following (19) and (11), while, as we explained above, the magnetic shear is set to 0.
Note that the value of effectively controls the strength of the PVG drive (the term proportional to on the right hand side of (19)) for a given velocity shear . The relatively low value of in the Cyclone Base Case reduces the destabilising effect of the PVG — compared for example to the Waltz Standard Case where (7).
Collisions are included by means of a model collision operator, which includes the scattering of particles in both pitch angle and energy (30); (31). The numerical ionion collision frequency .
The effects of the perpendicular flow shear are included in the code as follows (32). Expanding , where is the coordinate along the field line, the perpendicular flow shear terms in (19) can be eliminated by adding time variation to the radial wavenumber:
(21) 
so that
(22) 
All fluxes will be reported in dimensionless units by normalising them to the gyroBohm values
(23) 
The flow shear will be everywhere normalised to ; hence at this point we rescale accordingly:
(24) 
The resolution of the majority of simulations was — the number of gridpoints in the and spectral directions, in the spatial direction, the average number of pitch angles (i.e. the number of values of for a given , which varies with the parallel coordinate), the number of gridpoints in energy space and the two values of . This grid provided sufficient scale separation in both spatial directions perpendicular to the magnetic field for calculating the turbulent transport, and sufficient resolution in velocity space to calculate the velocity integrals in Maxwell’s equations with the required accuracy. A discussion of the parallel resolution is given in Section IV.
Iii Turbulent Fluxes
iii.1 Heat Flux
Using simulations in the manner described in Section II, the heat and momentum fluxes were calculated for values in the range , vs flow shear values in the range .
Fig. 1(b) shows the heat flux vs the ITG in the absence of flow shear. When , the critical temperature gradient above which there is ITGdriven turbulence is . When is increased above this threshold, the heat flux increases rapidly up to more than a hundred times the gyroBohm value.
Fig. 2 shows the turbulent heat and momentum fluxes vs the flow shear for different ITG values.
As the flow shear is increased from 0, the heat flux initially responds weakly,
either increasing or decreasing slightly.
As the flow shear approaches 1, the heat flux dips sharply for all values of the ITG.
For moderate temperature gradients (),
the turbulence is suppressed altogether and the heat flux drops to zero.
The suppression of turbulence also happens at finite magnetic shears (10),
but for a significantly narrower range of
For larger , the heat flux starts to rise again. This increase is due to the PVG — we have verified that this effect disappears if the term proportional to on the right hand side of eq (19) is artificially set to 0. This revival of the turbulent transport at large shears was also observed in (10), but was absent from the quasilinear study conducted in (5), which showed the turbulence being completely suppressed above a sufficiently large toroidal shear.
iii.2 Momentum Flux
The momentum flux is zero when (as might be expected in an updown symmetric case, (33); (34)). As increases, so does the momentum flux — a result of the turbulent viscosity. However, when the flow shear reaches values at which it starts to suppress turbulence significantly, the trend is reversed and a remarkable situation arises in which increasing flow shear reduces the transport of momentum. This behaviour persists until reaches even larger values, and the turbulence is reignited by the PVG drive. The positive correlation between the velocity shear and the momentum flux is then reestablished. It is the existence of the window of suppressed momentum transport at moderate and that will enable the transport bifurcation analysed in Section V.
iii.3 Turbulent Prandtl Number
There is a clear correlation between the heat and the momentum flux, which is best quantified in terms of the turbulent Prandtl number
(25) 
where the turbulent viscosity and the turbulent heat diffusivity are
(26)  
(27) 
respectively.
The Prandtl number obtained from our simulations is plotted in Fig. 2(c).
There is a similar basic trend for all values of :
rises from approximately 0.5 when
Although the Prandtl number does vary with both and , this dependence is relatively weak compared to the individual dependence of and on these parameters. Thus, approximating everywhere keeps it within approximately 25% of the true value for the entire range of . This will prove convenient in constructing a model of turbulent transport presented in Section VI.
Iv Subcritical Turbulence
Before studying the implications of the results presented above, we first discuss the nature of the fluctuations that give rise to the behaviour of the turbulent fluxes discussed in Section III.
With , the ITG instability exhibits robust linear growth, with a threshold of (Fig. 1(a)). However, for any finite value of the flow shear there are no growing linear eigenmodes in the system: all modes grow only transiently before decaying (Fig 3(a)). While formally this means that is a singular limit, there is is in fact no physical discontinuity in behavior: as , the duration of the transient growth, , tends to infinity (roughly as ; see Fig. 3(b)). There is therefore a continuous transition at from a transient mode that grows for an infinitely large time, to a growing eigenmode. The lack of growing eigenmodes when stems from the wavenumber drift (21): as time tends to infinity at finite flow shear and zero magnetic shear, the radial wavenumbers increase inexorably through the shearing of the modes, until they are extinguished by finiteLarmorradius effects.
This situation differs somewhat from the case with finite magnetic shear (10); (37), where there are linearly growing modes for a finite range of . At finite magnetic shear, growing eigenmodes are possible because while leads to an effective dependence of the radial wavenumbers on time, makes them dependent also on the position along the field line. Therefore, for modes moving with a speed proportional to , the magnetic shear cancels the effect of the velocity shear on (13); (4). When the magnetic shear becomes very small, or the flow shear very large, the required mode velocity becomes much greater than the thermal speed of the ions and the cancellation is no longer possible because the mode cannot travel so fast (13). As a consequence, growing linear eigenmodes only exist for nonzero magnetic shear and only up to a finite value of the flow shear (10).
In a standard picture of ITG turbulence with (19), the turbulence is driven by unstable linear modes. With the exception of a narrow interval of temperature gradients where selfgenerated zonal flows suppress the turbulence (the Dimits shift), the presence and the amplitude of the turbulence are largely determined by the presence and the growth rate of those unstable modes. In the present case, by contrast, we have strong levels of turbulence sustained nonlinearly in a parameter regime where there are no linearly unstable modes, a phenomenon first discovered in tokamak turbulence in (10) and (11), but well known in various hydrodynamical contexts (38). This phenomenon is known as subcritical turbulence.
Subcritical turbulence differs from standard instabilitydriven turbulence in several important ways. Firstly, because there is no linear instability, turbulence will not grow from initial perturbations of arbitrarily small amplitude (10). Fluctuations must be initialised with sufficient amplitude (generally of the order of their amplitude in the saturated state) in order for turbulence to be sustained; thus, the absence of sustained turbulence in a numerical experiment may merely indicate that the initial amplitude is insufficient, not that the plasma is quiescent.
The second difference concerns the scales at which the turbulent energy resides. In ITGdriven supercritical turbulence (), these scales are those where the linear drive injects energy — this tends to correspond to and relatively narrowly concentrated around values of a fraction of unity (at least for low values of and moderate temperature gradients; see (39) and Fig. 5(a)). In subcritical turbulence, the preferred wavelengths appear to be those at which the amplification of the transient modes is strongest. In the case of turbulence where the PVG drive is dominant, Ref. (14) has shown that the amplification is maximised along a line in Fourier space where
(28) 
Fig. 5(b) shows the spectrum
of strongly PVGdriven turbulence at high flow shear ().
We see that although the result (28)
is based on linear theory and
is derived for slab geometry,
it appears to describe the peak of the spectrum reasonably well.
The spectrum extends to much higher
parallel wavenumbers than in the standard ITG case;
consequently, higher parallel resolution is necessary to resolve it
Finally, faced with subcritical turbulence, we are left without an intuitively obvious way of estimating its saturation level and, consequently, the turbulent fluxes. It is the presence of linear eigenmodes with a defined growth rate (positive or negative) which has enabled the quasilinear modelling of the heat flux in many situations without resorting to full nonlinear simulations ((40) was an early exposition; see (41) for an overview and (42); (43) for recent work). When the growth of all modes is transient, such models will not work in their current form. The question arises as to which characteristics of the linear transient growth are relevant in the resulting nonlinear state. The problem is as yet unsolved, but here we consider two natural measures of subcritical driving.
Let us first define the effective transient growth rate (9)
(29) 
This effective growth rate is plotted in Fig. 4(a). It initially decreases with increasing flow shear, but then starts to increase for as the PVG drive becomes significant. The effective growth rate seems to correlate with the behaviour of the heat flux (Fig. 2(a)): it first decreases and then increases with increasing flow shear. However, the minimum in the effective growth rate (for as shown in Fig. 4(a)) occurs at whereas the minimum in the heat flux for the same occurs at .
Besides the rate of growth, what is likely to be decisive in determining whether the finite transient growth of perturbations proves sufficient to maintain a steady level of turbulence is the time that this growth lasts. A natural diagnostic that accounts for both effects is the total transient amplification of the modes before they start to decay (14); (38); (44):
(30) 
It drops precipitously as is increased from 0,
reaches a minimum at
(which does not coincide with , where is minimal)
and then appears to increase very gently with
Thus, while the linear behaviour shows that the nonlinear fluxes are somewhat determined by the vigour of the underlying linear drive, a quantitative model of subcritical PVG turbulence that would allow the heat fluxes to be predicted from the characteristics of transient linear growth remains an unsolved problem.
V Transport Bifurcation
v.1 Possibility of Bifurcation
In Section III,
we demonstrated the existence of a wide interval in
in which a sheared toroidal equilibrium flow completely suppresses turbulent transport.
If such a suppression could be achieved experimentally in a tokamak,
confinement of energy would be dramatically improved.
Unfortunately, while it is possible to specify the flow shear in numerical simulations,
in experiment it can only be varied indirectly by applying torque,
and the effect of that torque is strongly dependent
on how quickly the angular momentum escapes from the plasma.
If only a limited amount of torque can be injected and the momentum flux is too large,
it could be impossible to achieve flow shears that are large enough to suppress the turbulence.
Fig 2(b), however, suggests an intriguing possibility.
For all temperature gradients, the momentum flux at first increases,
reaches a maximum and then decreases
before increasing again at higher flow shears.
There is, therefore, a parameter window in which
increasing the flow shear decreases the transport of momentum.
If this were to happen, the momentum would build up,
increasing the flow shear, which would further decrease the transport of momentum.
Such a positive feedback could lead to a bifurcation and so it might seem that high values of flow shear could be reached without excessive input of momentum,
a possibility which was discussed in (45)
v.2 Inverting the Problem
What can actually be controlled in a steadystate experimental situation is the flow of heat and momentum through a particular surface. This means that we must switch from using and as independent parameters to using the total fluxes of heat and momentum, and respectively.
Let us consider an experimental setup where the heat and the momentum are injected by beams of neutral particles. We assume that the energy and momentum from those neutral beams are deposited uniformly across the torus. We further assume that the beams are tangential to the magnetic axis and we ignore corrections of order the aspect ratio of the device. Then:
(31)  
(32)  
(33) 
where is the number of neutral beam particles injected per unit time, is the beam particle velocity, is the beam power, is the beam particle energy and is the volume fraction of the plasma enclosed by the flux surface. Thus, the total heat flux is determined by the beam power (31), whereas the ratio of the total momentum angular flux to the total heat flux, , is determined by the beam particle energy (33). We want to know whether, by varying and , it is possible to reach, or to trigger a transition to, a highflowshear regime where the turbulent transport is suppressed. Thus, it is necessary to convert the dependence of and on and , determined from local simulations, to a dependence of and on and .
v.3 Interpolation
The gyrokinetic code GS2 gives the fluxes as a function of and . Given the expense of the nonlinear simulations necessary to do this, it is computationally challenging to use GS2 as a root finder to invert the problem and find the gradients as a function of the fluxes (46); (47). Instead we interpolate within the set of data points described in Section III (as well as additional simulations in parameter regions of particular interest) to obtain the fluxes as continuous functions of the gradients; these functions can then be inverted numerically to give and as functions of the fluxes.
Interpolation in a multidimensional parameter space is a nontrivial operation. A standard technique is to use radial basis functions (48), which weigh each data point by its distance in parameter space from the point of interest (after the parameter space has been normalised to ensure that variation occurs on the same scale in each coordinate). There are many choices of the function, or kernel, which is used to calculate the relative importance of each data point. Here we choose a linear kernel (equivalent to linear interpolation in the case of only two points) (48). Using this, the values of the fluxes at a point can be calculated as follows:
(34)  
(35) 
where are the input gradients for the nonlinear simulation labelled , and the weights and are calculated so as to satisfy for all :
(36) 
where and are the values of the fluxes obtained from simulation .
v.4 Neoclassical Transport
If the turbulent transport is successfully suppressed, neoclassical (collisional) transport becomes important. The total fluxes are the sum of the neoclassical and turbulent contributions:
(37) 
where the neoclassical fluxes are calculated as
(38)  
(39) 
and the neoclassical thermal diffusivity and viscosity are
(40)  
(41) 
The formulae for the neoclassical diffusivity and viscosity in the largeaspectratio, concentriccircularfluxsurface, banana () regime, were taken from (17). In (17), the ionion collision frequency is defined as:
(42) 
(where is the Coulomb logarithm). Thus, the neoclassical transport is a function of the temperature and density, and scales with these quantities in a different way to the turbulent transport. To determine the ratio between the neoclassical and turbulent transport it is necessary to choose specific values for the temperature and density. Here, as in (11), we take . Later, when we consider the case of specific tokamaks, we will use typical values of and to estimate .
Using the results (40) and (41) the neoclassical Prandtl number can be calculated from the parameters given in (20) and is found to be
(43) 
Thus the neoclassical Prandtl number is much smaller than the turbulent Prandtl number (see Fig. 2(c)). In other words, the neoclassical transport of momentum is much smaller than the neoclassical transport of heat, in contrast to turbulent transport, for which the fluxes of momentum and heat are comparable. While the formulae (41) and (40) are approximate, we emphasise that the qualitative result of the following section is not affected by small changes in the values of and , provided that the property continues to hold. In particular, this means that this qualitative result is not affected by small changes in the value of .
v.5 Bifurcation
In Fig. 6(a), the momentum flux is once more plotted against , but this time keeping the heat flux constant at . The local maximum in the momentum flux still exists. Starting at this local maximum (point A), if the torque on the plasma was to be increased at constant , the flux of momentum could only increase if the flow shear were to jump to a much higher value (point B) where the PVG instability would drive turbulent momentum flux. A bifurcation is manifest.
The inclusion of the neoclassical fluxes is critical in obtaining this result. To illustrate this, Fig. 6 also shows the momentum flux at constant heat flux without the neoclassical contribution to the fluxes, i.e., with . If, in such a synthetic “purely turbulent” transport case, the flow shear is increased at constant , the temperature gradient increases to maintain the turbulent heat flux, and, as shown in Fig. 6, increases monotonically with . The key property of neoclassical transport that helps change this behaviour and produce a bifurcation is that the neoclassical Prandtl number is much smaller than the turbulent Prandtl number. This means that as turbulent transport is suppressed, the system goes through a regime where the neoclassical contribution to the heat flux is significant, while the neoclassical contribution to the momentum flux is not. So as is increased at constant , the turbulence is suppressed, and the temperature gradient starts to rise, but as this happens more of the heat flux is transported neoclassically, and so the feedback loop breaks down: it is no longer necessary to increase the levels of turbulence to maintain the same . The same is not true of the momentum: the neoclassical viscosity cannot make up for the lack of turbulence, and so the transport of momentum peaks and then falls.
At high flow shears the turbulent momentum flux increases rapidly once again; this turbulent flux is driven by the PVG, as discussed in Section III. As a result, we observe that the bifurcation results in a turbulent state at B, with a significant turbulent momentum flux. This is a substantial difference from the situation envisioned in (5), where a transport bifurcation was predicted using a reduced quasilinear model. Their model did not (and, being a quasilinear model, could not) predict the existence of the PVGdriven subcritical turbulence at high flow shears; instead it predicted a bifurcation resulting in a nonturbulent state where all transport was neoclassical. Thus a full nonlinear analysis is necessary to describe accurately the reduced transport state produced by the bifurcation that we find in a turbulent plasma.
Restoring dimensions, we find that the bifurcation results in a jump in the flow shear from to . It is instructive to consider whether such values of shear appear in current devices. Taking the measured values of the perpendicular flow shear, and , in highperformance discharges in the JET and MAST tokamaks (20); (49), we estimate that flow shears of up to approximately have been recorded in both devices, and thus that values of shear of the order examined in this paper have been observed.
Vi Parameterised Model
Through simple interpolation, without making any assumptions about the way the fluxes depend on the gradients, we have shown the existence of a bifurcation in the gradients at constant fluxes. However, even interpolating one line of constant requires a very large number of data points in the region of the line. To produce the required data set for Fig. 6(a), we had to perform about 350 nonlinear gyrokinetic simulations, each run to saturation. In order to extend our understanding of the bifurcation, and of the range of flux values for which similar bifurcations can occur, we will first consider the behaviour of the turbulent fluxes in further detail and then use this analysis to construct a simple parameterised model of the turbulent fluxes and as functions of and .
vi.1 Modelling
Dependence of on
We wish to construct a parameterised model of as a function of and . In Section III we described the dependence of on ; we now consider the dependence of on . Both experimentally (50) and numerically (19), it is usually found that increases very sharply with — a property known as stiff transport. A recent experimental study (50) has indicated that increasing the flow shear at low magnetic shear might have the effect of reducing the stiffness. Fig. 7 shows that the effects of flow shear on are in fact quite complex. Let us consider the cases of low flow shear, , and high flow shear, , separately.
For the case of , shown in Fig. 7(a), the threshold in above which turbulence can be sustained nonlinearly increases rapidly with , as the perpendicular shear suppresses the ITG instability. Above this threshold there are broadly three regions: low, intermediate and high .
At high , far above the threshold, for any value of flow shear the heat flux eventually asymptotes to the same dependence as it has at . In other words, as the ITG drive becomes very strong, the effect of flow shear becomes negligible.
At intermediate , the heat flux rises rapidly from low values to join the universal, high asymptotic. As increases, because the threshold rises, the heat flux rises more rapidly with above the threshold; thus at intermediate values of , the stiffness increases with ; this is shown in Fig. 8.
At low , just above the threshold, the heat flux rises very slowly, that is, the stiffness is very low (see Fig. 8). This low region is shown in detail in Fig. 7(b). Such is the sharpness of the transition from the lowstiffness low region to the highstiffness intermediate region that there appear to be two distinct thresholds. The first threshold is the transition from no turbulence to nonzero turbulent transport. Above the first threshold turbulence is present but rises slowly with (the low region). Above the second threshold rises rapidly (the intermediate region). These thresholds are plotted in Fig. 9. The lowstiffness region only exists for : between the first threshold joins the second and the lowstiffness region disappears.
At high flow shear, , there are two principal differences to the case with low flow shear. Firstly, the lowstiffness region is not present, and there is only one threshold. Secondly, the critical temperature gradient for turbulent heat transport starts to decrease as the PVG drive reenforces the ITG drive. When , the PVG drive is strong enough to drive turbulence at very low : the threshold drops to zero. As the threshold drops, the transport stiffness at intermediate values of decreases, in a mirror image of the case at low flow shear. At high , the heat flux still asymptotes to the universal curve.
Thus, increasing can both increase and decrease the nonlinear thresholds, and both increase and decrease the stiffness. The rise and fall both in the thresholds and the stiffness can also be seen in the finitemagneticshear simulations of (10). The principal differences between the zeromagneticshear case considered here and that case are that we have a higher value of the critical gradients for all , and that we find a lowstiffness region at low values of — a feature that seems to be absent at finite magnetic shear.
Parameterisation of
In (22) a simple model was used to characterise the behaviour of the turbulent heat flux, and describes the qualitative behaviour of the heat flux reasonably well. Here, we describe a more complex model that reproduces most of the features described above and is quantitatively close to the interpolated fluxes. We consider only the low and intermediate regions, as the bifurcation occurs near the boundary between these regions. In order to describe in these regions, we have to parameterise the two thresholds and the transport stiffness .
We parameterise the first and second critical thresholds as linear and quadratic functions of respectively:
(44)  
(45) 
where is the nonlinear threshold for turbulence at and the parameters , and are chosen to fit the data. The modelled thresholds are plotted next to the measured thresholds in Fig. 9. As with the observed thresholds, the first threshold joins the second between .
Next we parameterise the transport stiffness . The measured values of in the first and second regions are shown in Fig. 8. Between the first and second thresholds (the low region), is modelled as a constant, . In the intermediate region, remembering the observation that broadly rises and falls with the nonlinear thresholds, we allow the stiffness to depend on . Thus,
(46)  
(47) 
The parameters are , and . The model of is shown in Fig. 8.
Thus is parameterised as a piecewise linear function of . It is zero below the first threshold, has gradient above the first threshold and gradient above the second:
(48)  
The model is compared with the data points in Fig. 10. To reproduce the bifurcation of Section V.5 in a quantitatively correct manner, it is in fact sufficient to represent accurately the two thresholds and the low stiffness region. However, by also parameterising the variation of , we have provided a model which uses six parameters to describe the (complicated) behaviour of the heat flux over a wider parameter regime with reasonable accuracy.
vi.2 Modelling
While the turbulent Prandtl number does vary with and , as discussed in Section III.3, this variation is relatively weak. Thus we model by assuming a constant turbulent Prandtl number:
(49) 
The same choice was made in (22).
vi.3 The Modelled Bifurcation
In Fig. 6(b) the model is used to replot the angular momentum flux at constant . It has the same shape as the interpolated curve in Fig. 6(a), and is quantitatively very close to it at low flow shear (). The agreement at higher flow shears is less good, which reflects the fact that the second critical gradient is not really a quadratic. Nonetheless, we consider this degree of precision adequate. The model is used below to explore further properties of the transition without the need for extremely large numbers of additional simulations.
Vii The Reduced Transport State
In Section V, a transition was described leading to a reduced transport state. Here we use the parametric model that was designed in Section VI to describe the properties of this state.
vii.1 Heat Flux at Constant
Fig. 6 showed the effect of varying whilst keeping constant. In contrast, Fig. 11 demonstrates the effect of varying whilst keeping constant. At high and low temperature gradient (marked (I) in Fig. 11(b)), increasing the heat flux has the effect of increasing the temperature gradient as might be expected: this is the ordinary turbulent regime. However, as is lowered below there arises a counterintuitive situation where decreasing has the effect of raising the temperature gradient. This is of course the combined effect of the neoclassical Prandtl number and the flow shear as described in Section V, and also discussed in detail in (22). In this region (marked (II) in Fig. 11(b)), becomes comparable to , but as the heat flux is lowered, the system cannot drop to a neoclassical state because of the smallness of the neoclassical transport of momentum. Thus the constancy of the ratio (i.e., the large finite flux of momentum) keeps the system in a turbulent state; the temperature gradient and the flow shear rise, and the curve representing vs becomes flatter and curls up in order to maintain its distance from the neoclassical line. At large temperature gradients and flow shears (region (III) in Fig. 11(b)), the momentum flux increases rapidly due to the PVG drive and the high flow gradient (see Fig 2(b)) and so the curve rolls over and resumes its original downward trend. Eventually the heat flux asymptotes to the neoclassical value.
The last part of this trend, while physically obvious, can only be seen using the modelled heat flux (Fig 11(b)), as it is not feasible to interpolate in the region near the neoclassical line where the contours of constant are very closely spaced and becomes multivalued.
vii.2 Temperature Gradient Jump
Interpolation cannot yield the temperature gradient after the transition directly: as explained above, the contours of constant become too closely spaced as they approach the neoclassical line in Fig. 11(a). However, the temperature gradient of the final state, labelled B on both Fig. 6 and Fig. 11(a,b), can be calculated indirectly by using the value of from Fig. 6(a) and rearranging (25) and (37) to give as a function of , , and :
(50) 
Taking , yields a temperature gradient at point B of . The temperature gradient at point A (the other side of the jump) is 7.4. Alternatively, using our model, the values of the temperature gradient can be read from Fig. 11(b) which results in an identical jump of to . If we compare with the case of zero momentum input and zero flow shear (point C on Fig. 11(a)), we find that flow shear has enabled a total increase in the temperature gradient of a factor . If this were reproduced across the whole device, it could increase the core temperature by a factor of around 10, but this would require a more detailed 1D model (see (51)), rather than the 0D model presented here. Note that while the jump in temperature gradient caused by the bifurcation (the transition from A to B) is a striking feature, a comparable contribution to the increase of the temperature gradient arises from the incremental suppression of the turbulent transport by the velocity shear (the difference between A and C).
vii.3 Neoclassical Heat Flux, Turbulent Momentum Flux
It was noted earlier that the reduced transport state produced by the bifurcation is still turbulent, with a momentum flux far greater than its neoclassical value (Fig. 6). In Fig. 11(b), it can also be seen that the reduced transport state does not lie exactly on the neoclassical line; it does however lie very close to it, which implies that the turbulent heat flux is much smaller than neoclassical. Thus, the bifurcation takes the system into a state where the heat is mostly transported neoclassically, but the angular momentum is mostly transported by turbulence. The small ratio of to reflects the fact that the turbulence levels in the reduced transport state are small. The dominance of over given the low levels of turbulence reflects the fact that the flow gradient is large and the neoclassical momentum transport is very low compared to the neoclassical heat transport ().
vii.4 Bifurcation by lowering
Starting from point A, if were to be reduced at constant , the system would again be forced to jump to point B. In effect, what would happen is that the heat flux would become principally neoclassical; the momentum flux would drop because , and a bifurcation would ensue in the manner described in Section V.5. Thus, a decrease in the input heat could lead to a higher temperature gradient. We note, however, that is normalised to the gyroBohm value: . Thus, decreasing could correspond to decreasing the deposition of heat, but it could also correspond to an increase in temperature or density (51).
vii.5 Optimum
Fig. 11(b) shows that for every value of , there is an optimum that gives a maximum temperature gradient (the dotted line at in Fig. 11(b)). The maximum temperature gradient increases with , but only slowly. The optimum value is the value of such that:
(51) 
and is plotted as a function of in Fig. 12.
vii.6 Transition Region
Finally, we will show that there is in fact only a limited range of both and where bifurcations can happen in the way described above. To illustrate, we observe that no transition can occur along the contour in Fig. 11(b): if is kept constant at this value, decreasing from greater than its optimum value of increases the temperature gradient smoothly up to its maximum value. The existence of a bounded region where such transitions can happen is studied in detail by the authors of (22), who calculate, using a simple transport model, the region in which transitions occur, and derive a criterion necessary for their existence. We apply the analysis of (22) using the model for the turbulent fluxes given in (48) and (49), to determine the range of values of and for which bifurcations of the form we have described can occur. In essence, bifurcations can only occur when there exist multiple values of and that give rise to the same values of and . From Fig. 11(b), it is clear that this is only possible for values of where there is a minimum in the curve of constant , and for values of which lie between this minimum and either a maximum in the curve or the point where the curve intercepts the neoclassical line. Thus the region in which transitions are possible is bounded by the curve
(52) 
and by the neoclassical line. This region is plotted in Fig. 12.
In order to give a clearer meaning to this diagram, we use equations (31  33) and typical properties of plasma devices taken from the ITER Profile Database (52) (and listed in Table 1) to replot the region in which transitions can happen in terms of the input beam power and beam particle energy, and . To generate these plots, we also calculate the collision frequency selfconsistently using (42), to take account of the varying strength of the neoclassical transport in each device. The transition regions for each device are displayed in Fig. 13, along with a point indicating the typical values of and in each device. Fig. 13 shows that in which transitions can happen lie within an order of magnitude of the typical values. It should be stressed that with the simplified magnetic geometry (Section II) and model of neoclassical transport (Section V.4) used in this study, and with the many assumptions about the way the energy and momentum are deposited in the plasma (Section V.2), closer agreement could not be expected. In particular, the assumptions of Section V.2 are likely to overestimate the applied torque and hence overestimate the beam energy needed for a transition.
()  ()  ()  ()  ()  

DIIID  2.7e+03  9.0e+19  1.7e+00  1.7e+00  6.1e01 
JET  2.6e+03  6.3e+19  2.9e+00  2.4e+00  9.3e01 
MAST  0.6e+03  3.9e+19  8.0e01  5.3e01  5.6e01 
Viii Conclusions
We have determined the way the turbulent transport of heat
and toroidal angular momentum
is affected by the equilibrium gradients of temperature and velocity
over a wide range of those two parameters,
in the case of a plasma with a sonic background flow

The plasma is linearly stable to microinstabilities for all nonzero values of flow shear but subcritical turbulence can be sustained nonlinearly across much of this stable region.

For velocity shear (we remind our readers that is normalised to ), increasing reduces the levels of turbulence driven by the ion temperature gradient. For , the destabilising effect of the parallel velocity gradient becomes dominant and increasing the flow shear increases the levels of turbulence and turbulent transport.

For low flow shears, , and low heat flux, , flow shear reduces the thermal transport stiffness. For low flow shears and higher , flow shear increases the thermal transport stiffness. At high flow shears, , the thermal stiffness is reduced as the PVG drive becomes dominant.

For the turbulent Prandtl number stays within the range 0.5–0.8 across a wide range of temperature gradients. Thus, the turbulence transports heat and momentum with comparable vigour.

For , there is a large region around where the turbulence is completely quenched.

Owing to the existence of this region, and the fact that the neoclassical Prandtl number is much lower than the turbulent Prandtl number, it is possible to trigger a bifurcation to a regime of high velocity and temperature gradients by either (i) increasing the momentum input at constant heat input or (ii) lowering at a constant ratio of heat to momentum input.

In the high gradient (reduced transport) state that is reached via this bifurcation, the heat is principally transported neoclassically, whereas the momentum is principally transported by turbulence.

For any given input of toroidal angular momentum, there is an optimum input of of heat which maximises the temperature gradient. Increasing the input of heat above this optimum can reduce the temperature gradient by increasing the levels of turbulence.
Acknowledgements.
We are grateful for helpful discussions with I. Abel, P. de Vries, W. Dorland and G. W. Hammett. This work was supported by EPSRC (EGH, FIP), STFC (AAS), EURATOM/CCFE Association (MAB,CMR,SCC) and the Leverhulme Network for Magnetised Plasma Turbulence. The views and opinions expressed herein do not necessarily reflect those of the European Commission. Computing time was provided by HPCFF and by EPSRC grant EP/H002081/1. Much of this work was carried out at the Isaac Newton Institute for Mathematical Sciences, Cambridge, UK, during the programme “Gyrokinetics for Laboratory and Astrophysical Plasmas”.Footnotes
 For the density, this is only true if the toroidal flow is subsonic, which is what we assume below, see Eq. 6.
 The definition of the thermal velocity given in (1) is chosen in accordance with the analytical works (18); (17) cited in this paper; maintaining correspondence with (10); (11), where the ion thermal velocity was defined as , results in various factors of in the normalisations of the output from simulations.
 With the additional assumption that the likeparticle collision frequency satisfies .
 This can be systematised by assuming that the magnitude of the flow is small but the gradient of the flow is large in a Machnumber ordering subsidiary to (4), namely (6) where .
 It should be noted that because the sign of in this paper is opposite to the sign of the poloidal field in GS2, the coordinate , defined in (17), also has the opposite sign to the GS2 coordinate of the same name.
 This is partly owing to the fact that at finite magnetic shear growing linear eigenmodes exist for nonzero flow shear, whereas at zero magnetic shear there are no such eigenmodes except when the flow shear is also zero; this is discussed further in Section IV.
 The low value of for small has been observed before and is sometimes referred to as a shear pinch (35); (36). It occurs because, at low , the perpendicular flow shear can give rise to a contribution to the viscous stress that has the opposite sign to that of ITGdriven turbulence with zero flow shear, reducing the overall diffusive transport.
 Studies of the effect of increasing the parallel resolution showed that at values of flow shear , the parallel resolution used for the parameter scan described in Section III (14 gridpoints) was insufficient and led to errors in the critical temperature gradient (above which subcritical turbulence could be sustained) of order 510%. We do not consider this error significant enough to merit a repeat of the parameter scan with substantially higher parallel resolution (although in principle such an exercise would be useful) .
 This behaviour approximately agrees with the analytical result of (14), which predicts as .
 It should be noted that as the sign of