Transport Bifurcation Induced by Sheared Toroidal Flow in Tokamak Plasmas

Transport Bifurcation Induced by Sheared Toroidal Flow in Tokamak Plasmas


First-principles 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 zero-magnetic-shear 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 reduced-transport 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 quasi-linear 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.


I Introduction

A major obstacle to the development of a working magnetic confinement fusion device is the transport of heat by pressure-gradient 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 high-flow 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 1 and temperature are constant on each flux surface. In the current study, which considers the Cyclone Base Case regime (19), the magnetic flux surfaces are concentric and circular in cross-section, so that , where is the minor radius of the torus and the poloidal magnetic field. Thus, may be used as the flux label, and the magnetic field may be written as .

We allow an equilibrium plasma flow , of the same order as the ion thermal velocity


where is the temperature of the ions and is their mass 2. It can be shown that such a flow is toroidal, as any poloidal component will be quickly damped (16); (17). The flow can then be expressed as , where is the angular velocity of the flow (which must be constant on a flux surface).

The state of each species of charged particles can be described by its distribution function, , whose evolution is given by the Vlasov-Landau equation:


where is the velocity coordinate and the collision operator. In order to describe the turbulent fluctuations which give rise to the transport, the Vlasov-Landau 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


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:


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 3 that to lowest order, is a local Maxwellian in the frame of the equilibrium flow , and that, to first order in ,


where is the charge of species and is the perturbed electrostatic potential. The non-Boltzmann 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 non-negligibly (4)4.

Under these combined assumptions, it can be shown that evolves according to the following gyrokinetic equation:


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:


where is a gyroaverage at constant . The electrons are treated through a modified Boltzmann response (26):


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,


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:


where is the equilibrium density scale length, the equilibrium temperature scale length. The perpendicular velocity shear and the magnetic safety factor are:


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


and the second perpendicular coordinate as


where is the poloidal angle5, then transforming to a frame rotating with angular frequency and keeping terms to the correct order in the gyrokinetic expansion, we find for any perturbed quantity :


The gyrokinetic equation (7) may now be written


In this investigation, we vary the perpendicular flow shear , and . Other parameters are kept fixed:


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 ion-ion 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:


so that


All fluxes will be reported in dimensionless units by normalising them to the gyro-Bohm values


The flow shear will be everywhere normalised to ; hence at this point we rescale accordingly:


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 ITG-driven turbulence is . When is increased above this threshold, the heat flux increases rapidly up to more than a hundred times the gyro-Bohm 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 6.

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.

Figure 1: ITG instability and turbulence at zero flow shear: (a) linear growth rate; (b) turbulent heat flux vs the normalised temperature gradient . The case of , further explored in Fig. 3, is marked by *.
Figure 2: Turbulent heat (a) and toroidal angular momentum (b) fluxes (normalised to gyro-Bohm values) as functions of flow shear for different values of the ion temperature gradient ; (c) turbulent Prandtl number vs flow shear (for cases where the heat flux is non-zero).

iii.2 Momentum Flux

The momentum flux is zero when (as might be expected in an up-down 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


where the turbulent viscosity and the turbulent heat diffusivity are



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 7, peaks at and then drops to approximately 0.6 when . For those intermediate values of where the turbulent fluxes are reduced or tend to 0, rises sharply, reaching just above 0.7 for low values of . The location of this sharp rise varies with temperature gradient similarly to the location where the turbulence is suppressed.

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.

Figure 3: Transiently growing linear modes at . (a) The heat flux as a function of time, normalised to its value at , where , so chosen to skip the short initial transient associated with a particular choice of initial condition. All modes grow transiently for . (b) Duration of transient growth from to the peak value of vs flow shear.
Figure 4: Two measures of the strength of the transiently growing linear modes at : (a) the effective growth rate of the mode and (b) the number of exponentiations of the heat flux during the growing phase, both vs. the flow shear. The apparent discontinuity at in (a) is a result of taking the average growth rate of the heat flux over rather than the initial or peak growth rate; for all the average will include a period where the growth rate tends to zero.

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 finite-Larmor-radius 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 non-zero 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 self-generated 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 instability-driven 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 ITG-driven 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


Fig. 5(b) shows the spectrum of strongly PVG-driven 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 it8.

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)


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):


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 9 — in sharp contrast with the very rapid increase of the nonlinear fluxes at high flow.

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.

Figure 5: The spectrum of turbulent fluctuations for normal and subcritical turbulence: (arbitrary units) vs and for (a) (ITG turbulence with no flow shear) and (b) (subcritical, strong PVG-driven turbulence)10. Also shown in (b) is the line of maximum transient amplification (28), as calculated in (14).

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) 11. However, Fig 2(b) shows the momentum flux at constant ion temperature gradient, which would not, in fact, stay constant during this process. Indeed, as soon as flow shear began to suppress the turbulent transport of momentum, it would also suppress the turbulent transport of heat, causing an increase in the temperature gradient as well as in the flow gradient. This increase in the temperature gradient would restore the turbulence to its former levels. Nevertheless, we will show in this section that when neoclassical transport is also taken into account, a bifurcation is possible.

v.2 Inverting the Problem

What can actually be controlled in a steady-state 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 set-up 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:


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 high-flow-shear 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:


where are the input gradients for the nonlinear simulation labelled , and the weights and are calculated so as to satisfy for all :


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:


where the neoclassical fluxes are calculated as


and the neoclassical thermal diffusivity and viscosity are


The formulae for the neoclassical diffusivity and viscosity in the large-aspect-ratio, concentric-circular-flux-surface, banana () regime, were taken from (17). In (17), the ion-ion collision frequency is defined as:


(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


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

Figure 6: Momentum flux divided by the heat flux (each normalised by the corresponding gyro-Bohm estimate) vs the flow gradient at a constant value of the heat flux plotted using (a) interpolation from the data explained in Section V.3 and (b) the parameterisation from Section VI. Also shown are the neoclassical contribution to the momentum flux (dotted line) and the momentum flux at constant heat flux without the neoclassical contribution (dashed line).

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 PVG-driven subcritical turbulence at high flow shears; instead it predicted a bifurcation resulting in a non-turbulent 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 high-performance 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

Figure 7: Turbulent heat flux vs the temperature gradient for different values of the flow shear, showing (a) low-shear values , (b) a close up of the low region in (a), and (c) high-shear values

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 low-stiffness low- region to the high-stiffness intermediate- region that there appear to be two distinct thresholds. The first threshold is the transition from no turbulence to non-zero 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 low-stiffness region only exists for : between the first threshold joins the second and the low-stiffness region disappears.

At high flow shear, , there are two principal differences to the case with low flow shear. Firstly, the low-stiffness 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 re-enforces 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 finite-magnetic-shear simulations of (10). The principal differences between the zero-magnetic-shear case considered here and that case are that we have a higher value of the critical gradients for all , and that we find a low-stiffness region at low values of — a feature that seems to be absent at finite magnetic shear.

Figure 8: The dependence of the heat transport stiffness on the flow shear in the low- and intermediate- regions, measured using simulations close to both thresholds (points, see Fig. 7) and using the parameterised model of Section VI.1.2 (lines).
Figure 9: The first and second critical thresholds, and , measured using simulations close to both thresholds (points, see Fig. 7), and using the parameterised model of Section VI.1.2 (lines).

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:


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,


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:


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.

Figure 10: The modelled heat flux (lines) shown along with the simulated data points for (a) low flow shear and (b) high flow shear. Legends as in Figures 7(a) and 7(c).

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:


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

Figure 11: Heat flux vs the temperature gradient at a constant ratio of the momentum flux to the heat flux (in units of ) plotted using (a) interpolation from the data (Sec V.3) and (b) the parameterised model (Sec VI). Also shown is the neoclassical contribution to the heat flux. Interpolation in (a) is impossible near this neoclassical line, where the contours are closely spaced and is multivalued. (b) also shows the maximum possible temperature gradient at low for a given (labelled “max ”). Fig. (c) plots the same curves as (b), showing both against and , and the curves projected on the - plane, illustrating the increase in the flow shear along each curve of constant . Points A and B in (a) correspond to points A and B in Fig. 6(a).

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 :


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 gyro-Bohm 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:


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


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 self-consistently 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.

Figure 12: The region in which bifurcations can occur, calculated using the analysis of (22) and the parameterisation of the fluxes described in Section VI. The cross indicates the location of the bifurcation described in Section V.5. The dashed line represents the optimum value of for a given value of , Eq. (51).
Figure 13: The regions (shaded) in which transitions can happen, as a function of the beam power and the beam particle energy, for three different devices. The dashed lines indicate, for each device, the value of for a given which would lead to the optimum temperature gradient, as described in Section VII.5. The points indicate typical values of and for each device. The projected and for MAST Upgrade were taken from Ref. (53).
() () () () ()
DIII-D 2.7e+03 9.0e+19 1.7e+00 1.7e+00 6.1e-01
JET 2.6e+03 6.3e+19 2.9e+00 2.4e+00 9.3e-01
MAST 0.6e+03 3.9e+19 8.0e-01 5.3e-01 5.6e-01
Table 1: Typical plasma properties from the ITER Profile Database (52). The symbol denotes the minor radius of the device. The temperature was calculated as the mean of the core (TI0) and edge (TI95) temperatures.

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 flow12. We have extended the range of flow shears up to 20 times the linear ITG instability growth rate calculated at zero flow shear. We have considered the zero-magnetic-shear regime. We have used a low value of the safety factor, , which reduces the strength of the destabilising parallel velocity gradient drive relative to the stabilising perpendicular velocity shear. In so doing we have established the following key results:

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

  2. 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.

  3. 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.

  4. 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.

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

  6. 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.

  7. 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.

  8. 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.

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 HPC-FF 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”.


  1. For the density, this is only true if the toroidal flow is subsonic, which is what we assume below, see Eq. 6.
  2. 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.
  3. With the additional assumption that the like-particle collision frequency satisfies .
  4. This can be systematised by assuming that the magnitude of the flow is small but the gradient of the flow is large in a Mach-number ordering subsidiary to (4), namely (6) where .
  5. 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.
  6. This is partly owing to the fact that at finite magnetic shear growing linear eigenmodes exist for non-zero 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.
  7. 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 ITG-driven turbulence with zero flow shear, reducing the overall diffusive transport.
  8. 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 5-10%. 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) .
  9. This behaviour approximately agrees with the analytical result of (14), which predicts as .
  10. It should be noted that as the sign of