Subcritical fluctuations in rotating plasmas

Suppression of turbulence and subcritical fluctuations in differentially rotating gyrokinetic plasmas

Abstract

Differential rotation is known to suppress linear instabilities in fusion plasmas. However, numerical experiments show that even in the absence of growing eigenmodes, subcritical fluctuations that grow transiently can lead to sustained turbulence, limiting the ability of the velocity shear to suppress anomalous transport. Here transient growth of electrostatic fluctuations driven by the parallel velocity gradient (PVG) and the ion temperature gradient (ITG) in the presence of a perpendicular () velocity shear is considered. The maximally simplified (but most promising for transport reduction) case of zero magnetic shear is treated in the framework of a local shearing box approximation. In this case there are no linearly growing eigenmodes, so all excitations are transient. In the PVG-dominated regime, the maximum amplification factor is found to be with (safety factor/aspect ratio), the maximally amplified wavenumbers perpendicular and parallel to the magnetic field are related by , where is the ion Larmor radius, the ion thermal speed and the shear. In the ITG-dominated regime, is independent of wavenumber and , where is the ion-temperature scale length. Intermediate ITG-PVG regimes are also analysed and is calculated as a function of , and . Analytical results are corroborated and supplemented by linear gyrokinetic numerical tests. Regimes with for all wavenumbers are possible for sufficiently low values of ( in our model); ion-scale turbulence is expected to be fully suppressed in such regimes. For cases when it is not suppressed, an elementary heuristic theory of subcritical PVG turbulence leading to a scaling of the associated ion heat flux with , , and is proposed; it is argued that the transport is much less “stiff” than in the ITG regime.

pacs:
52.30.Gz, 52.35.Qz, 52.35.Ra, 52.55.Fa

1 Introduction

It has long been understood that anomalous transport in tokamaks is caused by turbulence at or just above the ion Larmor scale, which is powered by drift instabilities extracting free energy from background gradients. The ion-temperature-gradient (ITG) instability [40, 15, 31] and the resulting turbulence [17] have been identified as a particular culprit. With the advent first of gyrofluid and then gyrokinetic numerical simulations, ITG turbulence has been the main focus of numerical studies aiming to map out levels of turbulent transport expected for any given set of equilibrium parameters (mean profile gradients and the magnetic configuration) [30, 18, 29, 8]. These studies assumed that the differential (toroidal) rotation of the tokamak plasmas, caused by momentum injection from the heating beams, had a negligible effect on the turbulence and on the resulting transport. This, however, is known to be an inadequate approximation for many devices and configurations in which the rotational shear can be comparable to the typical turbulent rate of strain. Linear eigenmode analysis [3, 5, 20, 14] suggests that the ITG growth rates are reduced and can be completely quenched by shear. A reduction or even complete suppression of associated transport was therefore expected and indeed found numerically, at least for some parameter regimes [28, 39, 9] — a result that kindled high hopes for controlling turbulence with shear and achieving transport bifurcations to steeper temperature gradients [21, 48].

Since the tokamak rotation is (to lowest order in the gyrokinetic expansion) purely toroidal [27, 16, 10, 13, 2], strong perpendicular shear comes at the price of a (stronger by a factor of ) parallel velocity gradient (PVG) — which is itself a source of free energy and so can trigger a drift instability [11], which in turns gives rise to turbulence and anomalous transport [47, 19, 7]. This instability is in fact also suppressed by the shear, but only in the sense that no unstable eigenmodes survive at large enough velocity and small enough magnetic shear: transient linear amplification is still possible [45, 46, 35]. This transient amplification (illustrated in figure 1) can be sufficient to give rise to subcritically excited turbulence [7, 25, 26]. It is this turbulence that limits the effectiveness of differential rotation in suppressing anomalous transport. While transport bifurcations are still possible [25, 26, 37], finding them in the parameter space turns out to be quite a delicate task, mostly because the subcritical PVG turbulence is a largely unexplored phenomenon. In particular, it is not amenable to the usual “quasilinear” mixing-length-type arguments because those are based on the calculation of linear growth rates — and it is not obvious what should replace them for subcritical turbulence.

Figure 1: Time evolution of the squared amplitude (normalised to its initial value) of a pure-PVG-driven linear perturbation, obtained in a direct linear numerical simulation using the gyrokinetic code AstroGK [36] with , , , , . The effective growth rate for this case is shown in figure 4.

We take the rather obvious view that addressing this problem should start from a systematic understanding of the linear transient amplification of gyrokinetic fluctuations in the presence of velocity shear. The purest example is presented by the limit of zero (negligibly low) magnetic shear, which has the twin advantages of analytical simplicity and of being most amenable to transport bifurcations or at least reduced levels of turbulent transport — as suggested both by numerical experiments [25, 26] and laboratory measurements [32, 33]. Formally, this regime supports no growing eigenmodes at any finite value of velocity shear, so we can focus on transient amplification and its dependence on the parameters of the problem (the ion temperature gradient, the velocity shear, the safety factor and the device aspect ratio ) without the complications of dealing with the transition to linear eigenmode stability. This is the kinetic treatment of this problem, following in the footsteps of the fluid theory [35] (which was done for a finite magnetic shear; note that zero magnetic shear is formally a singular limit).

Since much of the transport modeling for fusion devices relies on gyrokinetic simulations, it is important to ascertain that the linear behaviour predicted analytically is reproducible in direct gyrokinetic simulations with standard numerical codes. In what follows, we do this using the AstroGK code [36], which, for all practical purposes, is the slab version of the widely used fusion code GS21 (running GS2 itself in slab mode produces similar results). The code solves (2) and (3) with nonlinearity switched off and a collision operator [1, 6] added to provide small-scale regularisation in phase space. These simulations demonstrate the degree to which our asymptotic results represent an adequate description of realistic parameter regimes and how the transition from the short-time to the long-time limit occurs, subject also to (minor) modifications of our results by collisions, finite resolution and numerical inaccuracies.

The rest of the paper is organised as follows. In § 2, we introduce gyrokinetics in a shearing box — the governing equations of the minimal model we have chosen to treat. In § 3, we consider a further simplified limit in which the ITG drive is negligible compared with the PVG drive and so we can concentrate on the essential properties of the latter, namely, work out how the transient amplification time and the amplification exponent depend on the velocity shear, and and which wavenumbers prove most prone to being amplified. In § 4, we generalise these considerations by including the ITG. Both analytical and (linear) numerical results are presented throughout. A qualitative summary of the linear regime and a comparison with the fluid limit treated in [35] are given in § 5. Finally, in § 6, a criterion for the onset of subcrtitical PVG turbulence is proposed, followed by a very crude heuristic theory of this turbulence and of the resulting heat transport — these considerations provide a version of the standard mixing-length “quasilinear” arguments suitable for transiently growing fluctuations and a set of scaling predictions in the spirit of [8]. In particular, we give a semi-quantitative form to the argument that heat transport must be much much less “stiff” in the presence of velocity shear [32, 33, 25, 26] than in in the standard ITG regime because the turbulence is no longer driven by the temperature gradient.

2 Gyrokinetics in a shearing box

Consider a plasma in a strong mean magnetic field. On the assumption of axisymmetry, this field can be expressed in the form , where (magnetic flux) and are two scalar functions determined by solving the mean MHD equilibrium equations and is the azimuthal angle with respect to which symmetry is assumed. It can be shown (see [2] and references therein) that if such a plasma rotates with some velocity ordered in the gyrokinetic expansion as similar in size to the speed of sound, this mean rotation must be purely azimuthal, the same for both species and with an angular velocity that depends on the flux label only: , where is the radial coordinate. Thus, each of the nested toroidal flux surfaces rotates at its own rate and there is a velocity shear set up by the variation of with . In this paper, we will be concerned only with the effect on the plasma stability of this differential character of the rotation. Formally, this effect is isolated by assuming the Mach number ( is the ion thermal speed) to be moderately low. Mathematically, this can be cast as an expansion in (subsidiary to the gyrokinetic expansion) where, while the mean velocity is ordered subsonic, , its scale length is ordered as . Under this scheme, effects such the Coriolis or centrifugal motion (as well as other, more obscure, ones) that scale with the magnitude of the mean velocity are negligible, while velocity gradients are retained.

Let us consider the vicinity of some flux surface and introduce a local orthogonal Cartesian frame that is moving with the flux surface and in which is the cross-flux-surface (“radial”) coordinate, is the coordinate in the direction of the magnetic field and completes the orthogonal frame (). Then the velocity field can be replaced by a pure linear shear flow:

(1)

where is the poloidal magnetic field (see A.1 for details; figure 2 illustrates the geometry). This model might be called the “flying slab approximation,” or, perhaps more conventionally, the shearing box — a common analytical simplification in the fluid dynamics of differentially rotating systems [23].

Figure 2: Magnetic field, velocity field (sheared) and the local Cartesian frame. The reference flux surface is shaded.

We will make three further simplifying approximations by assuming all fluctuations to be electrostatic, electrons to have purely Boltzmann response and the mean magnetic field to be (locally) straight and uniform, so it has neither curvature nor variation of magnitude nor shear (which means that the magnetic drifts can be dropped). With all these assumptions and working in the flying slab (shearing box), we can write the gyrokinetic system of equations in the following form (see A)

(2)
(3)

where is the ion gyrocentre distribution function, the equilibrium Maxwellian (of the ions), the nondimensionalised electrostatic potential ( is the ion charge), the ion Larmor radius, the ion number density, the particles’ peculiar velocity with respect to the mean flow, the position coordinate, the guiding-centre coordinate, and the angle brackets denote the gyroaverages holding constant the coordinate appearing in their subscript (the precise definitions are provided in A). We have defined the equilibrium gradients and two other standard parameters as follows

(4)

where and are the ion density and temperature, respectively, and is the azimuthal mean magnetic field.

Note that the velocity shear appears in two ways in (2): as perpendicular shear (multiplying ) and as parallel shear (multiplying in the right-hand side). The perpendicular shear rips apart the unstable fluctuations and will have a stabilising effect, whereas the parallel shear (the PVG) acts as a source of free energy in a manner analogous to ITG and drives a drift instability (§ 3.1). The relative size of these two effects is set by the value of .

2.1 Case of non-zero magnetic shear

Including a (locally) constant linear magnetic shear into the problem amounts to replacing in (2)

(5)

where is the scale length associated with the magnetic shear. This appears to introduce complications as we now have an “effective shear” that depends on the particle velocity . However, since the size of is constrained by the Maxwellian equilibrium distribution, this term can be neglected provided . Under this assumption, the theory developed below applies without modification. We note that the “shear Mach number” is precisely the parameter that is known to control linear stability and transient amplification of in the ITG-PVG-driven plasmas in the fluid (collisional) limit [35].

2.2 Shearing frame

The next step — standard in treatments of systems with linear shear — is to make a variable transformation that removes the shear terms ():

(6)

and similarly for . The Fourier transform can then be performed in the primed variables, so

(7)

where , , and (we denote ). As usual in the gyrokinetic theory, working in Fourier space allows us to compute the gyroaverages in terms of Bessel functions:

(8)

where , and we have shifted the origin of time: .

Finally, we rewrite the gyrokinetic system (2)–(3) in the new variables . Since we are interested only in the linear problem here, we will drop the nonlinearity. We also suppress all primes in the variables. The result is

(9)
(10)

where is the drift frequency and ; the argument of the Bessel function is .

2.3 Integral equation for the linearised problem

We integrate (9) with respect to time, assume the initial fluctuation amplitude small compared to values to which it will grow during the subsequent time evolution, rescale time , denote , use (10), in which the velocity integrals involving the Maxwellian are done in the usual way, and obtain finally the following integral equation for :

(11)

where is the normalised drift frequency, the normalised shear parameter, and

(12)

where , , and and are modified Bessel functions of the first kind.

Equation (11) is the master equation for the linear time evolution of the plasma fluctuations driven by the ITG (the term) and the PVG (the term).

3 Solution for the case of strong shear

We will first consider the maximally simplified case of pure PVG drive (strong shear). This is a good quantitative approximation to the general case if , which in terms of the basic dimensional parameters of the problem translates into

(13)

This is the regime into which the plasma is pushed as the flow shear is increased — under certain conditions, the transition can occur abruptly, via a transport bifurcation [25, 37, 26]. Besides being, therefore, physically the most interesting, this limit also has the advantage of particular analytical transparency (the more general case including ITG will be considered in § 4).

Thus, neglecting all terms that contain and , (11) becomes

(14)

3.1 Short-time limit: the PVG instability

Let us first consider the case in which the velocity shear is unimportant except for the PVG drive, i.e., we can approximate and so there is no time dependence in the Bessel functions in (14). We would also like to be able to assume so that the time integration in (14) can be extended to . Formally this limiting case is achieved by ordering and (in dimensional terms, this is equivalent to , and ). Physically, this regime is realised in the initial stage of evolution of the fluctuations or, equivalently, in the case of very weak shear but large .

Under this ordering, we can seek solutions to (14) in the form , where is the nondimensionalised complex frequency and its dimensional counterpart. The time integral in (14) can be expressed in terms of the plasma dispersion function, which satisfies [22]

(15)

With the aid of these formulae, (11) is readily converted into a transcendental equation for :

(16)

where , and .

Equation (16) is simply the dispersion relation for the ion acoustic wave modified by the PVG drive term. This point is probably best illustrated by considering the cold-ion/long-wavelength limit , , . Then , , and so, restoring dimensions in (16), we get

(17)

where in the last expression, we have restored dimensions and denoted , the sound speed. When is sufficiently large, the sound wave is destabilised and turns into the PVG instability. Note that it loses its real frequency in this transition.

   

Figure 3: Left panel: Growth rate (normalised by ) vs. . Here , and the two curves are for (red) and (brown). The growth rate becomes independent of for (see end of B.1) and all curves for large are very close to each other and similar in shape to (see right panel). The mode has no real frequency for ; for , it turns into a damped sound wave: the corresponding frequencies and the damping rate are shown only for , as thin blue (dashed) and red (solid) lines, respectively. The discrete points show growth rates calculated by direct linear numerical simulation using the gyrokinetic code AstroGK [36]. Right panel: Contour plot of vs.  and . Only positive values are plotted, black means . The red curve shows the stability boundary (64). The white line shows the boundary (65) between PVG modes (above it) and sound waves (below it).

Some further (elementary) analytical considerations of the PVG instability are relegated to B.1. Here it will suffice to notice that if the (dimensional) growth rate is scaled by and by , their mutual dependence is universal for all values of the velocity shear or , namely,

(18)

Once this rescaling is done, the dispersion relation (16) no longer contains any parameters (except , which we can safely take to be order unity). The growth rate, obtained via numerical solution of (16) with , is plotted in figure 3. The maximum growth rate is . This peak value is reached when and , i.e., at .

The conclusion is that, at least in the initial stage of their evolution, plasma fluctuations in a significant part of the wavenumber space ( and are constrained by (64)) are amplified by the PVG. This amplification does not, however, go on for a long time as the approximation we adopted to derive the dispersion relation (16) breaks down when (or if the dimensional units of time are restored) — this gives , so, realistically, after barely one exponentiation. The key question is what happens after that. We will see shortly that all modes will eventually decay and that the fastest initially growing modes are in fact not quite the ones that will grow the longest or get maximally amplified.

3.2 Long-time limit: transient growth

Let us now investigate the long-time limit, , (or, in dimensional form, , ). In this limit, the kernel involving the Bessel function in (14) simplifies considerably: we have , and so

(19)

Working to the lowest nontrivial order in , we can now rewrite (14) as follows

(20)

We will seek a solution to this equation in the form

(21)

where is the effective time-dependent growth rate (nondimensionalised) and the dimensional version of it (remember that time is scaled by ). Because of the exponential in the kernel, the memory of the time-history integral in the right-hand side of (20) is limited, so and we will be able to make progress by expanding

(22)

We will assume that this expansion can be truncated; the resulting solution will indeed turn out to satisfy , with all higher-order terms even smaller. Substituting (22) into (20), we obtain an implicit transcendental equation for :

(23)

This equation can be written in a compact form by invoking once again the plasma dispersion function (15): denoting , we get

(24)

— effectively, a time-dependent dispersion relation, reminiscent of the PVG dispersion relation (16).

Transient growth.

First of all, it is immediately clear from (23) that as time increases, (the real part of) the effective time-dependent growth rate must decrease and indeed go negative because the right-hand side has to keep up with the increasing left-hand side. Therefore, fluctuations will eventually decay. However, if is positive for some significant initial period of time, there can be a substantial transient amplification.

We can determine the time when the transient growth ends by setting in (24). We immediately find that as well and that

(25)

The dependence of on and — via and via the time normalisation factor of — tells us which modes grow longest. The interesting question, however, is rather which modes get maximally amplified during this transient growth.

Figure 4: Time evolution of the effective growth rate: the red (bold) line is obtained as a numerical solution of (26) and plotted vs.  (the time axis is logarithmic in base 10). The black (thin) lines are the asymptotics (27) and (32). The discrete points show the time evolution of the effective growth rate obtained in a direct linear numerical simulation using the gyrokinetic code AstroGK [36] with , , , , . The short-time-limit PVG growth rate for this case (obtained by solving (16)) is shown as a dotted horizontal line. The time is normalised using (25) for . The time evolution of the perturbation amplitude for this case is shown in figure 1.

Maximal amplification.

The total amplification factor is given by , where is the number of exponentiations experienced by the mode during its growth period. In order to determine this, we need to know the time evolution of up to . Using (25), it is convenient to rewrite (24) as follows

(26)

When , the solution is found by expanding the plasma dispersion function in :

(27)

This asymptotic is not valid when approaches . More generally, (26) has a solution , whose functional form is independent of any parameters of the problem. It is plotted in figure 4 together with the asymptotic (27) and with a direct numerical solution showing how the transition between the short-time (§ 3.1) and long-time limits occurs. The amplification exponent is easily found:

(28)

where we have used (25) for and computed the integral numerically. It is clear from (27) that the integral converges on its lower limit and is not dominated by it, so it does not matter that we cannot technically use (26) for short times.

      

Figure 5: Left panel: The amplification exponent vs.  and in the same AstroGK simulation as used to produce the discrete points in figure 4. The straight line shows the relationship between the wavenumbers given by the second formula in (29). Right panel: The maximum amplification exponent [ given by (28), maximised with respect to ] vs. . The dotted lines show the asymptotic [see (29)] and the asymptotic, the latter straightforwardly obtained from (28): , (but this is purely formal because the long-time conditions and will be broken in this regime, except for extremely long wavelengths). The discrete points show obtained via an AstroGK numerical parameter scan: , varying and holding all other parameters fixed as in figure 4 (note that the asymptotic results do not in fact depend on or , although the quality of the long-time asymptotic does).

According to (28), depends on both wavenumbers via only (a plot of this dependence will be given in figure 8). Assuming , we find that the amplification exponent is maximised for , giving2

(29)

These results are illustrated in the left panel of figure 5. The amplification time for the maximally amplified modes identified in (29) is, from (25),

(30)

where we have restored dimensions to make explicit the dependence of on .

Thus, we have learned that an entire family of modes, characterised by a particular (linear) relationship between and , given in (29), will eventually enjoy the same net amplification, even though, as follows from the results of § 3.1, they were not the fastest initially growing modes and some within this equally amplified family started off growing more slowly than others or even decaying. The more slowly growing modes are the longer-wavelength ones and, according to (30), they compensate for their sluggishness with longer growth times.

Note that (30) confirms that the long-time limit is analytically reasonable because for large , (30) formally satisfies and also , provided . The latter condition is marginally broken by the fastest initially growing modes: indeed, in § 3.1, we saw that they had , so for them, , not really a large number. For all longer-wavelength modes, is safely within the domain of validity of the long-time limit.

Note also that, as we show at the end of B.1, the time-dependent dispersion relation (26) and its consequences derived above can be obtained directly from the PVG-instability dispersion relation (16) simply by restoring its dependence, setting , taking the short-wavelength and long-time limit (, ), and assuming . This calculation underscores the fundamental simplicity of the physics of the transient amplification: perturbations initially destabilised by the PVG are eventually swept by the perpendicular velocity shear into a stable region of the wavenumber space.

Limits on short and long wavelengths.

We have seen that modes with parallel wavenumbers up to can be transiently amplified. From (29), we conclude that of these, the maximally amplified ones will have perpendicular wavenumbers up to , i.e., can be relatively large — unlike in the short-time limit treated in § 3.1, where the modes with grew the fastest (although large were also unstable).

It should be understood that, while there is no ultraviolet cutoff in our theory that would limit the wavenumbers of the growing modes (in either direction), such a cutoff does of course exist in any real system. In the parallel direction, those that were strongly damped in the short-time limit (see (64)) are unlikely to recover in the long-time limit. In the perpendicular direction, the cutoff in will come from the collisional damping, which, in gyrokinetics, contains a spatial diffusion (see, e.g., [1]), and from the electron Landau damping, which we have lost by using the Boltzmann electron response (see (3)) and which should wipe out large and .

On the infrared (long-wavelength) side of the spectrum, we have no cutoffs either. In a slab, these would be provided by the dimensions of the periodic box. In a real plasma, the cutoffs are set by the scales at which the system can no longer be considered homogeneous (in a tokamak, these are the equilibrium-gradient scale lengths and the minor radius for the perpendicular scales and the connection length for the parallel scales; we will need these considerations to fix transport scaling in § 6).

Significant amplification threshold.

If we maximise (28) without assuming (with the caveat that the long-time limit asymptotics are at best marginally valid then), we obtain a more general curve than (29), plotted in the right panel of figure 5. We may define a critical threshold for significant amplification: when . The role of this threshold will be discussed in § 6.1.

Long-time decay.

Finally, we obtain the long-time asymptotic decay law. Let us seek a solution of (26) such that and . Then

(31)

Since is only root-logarithmically large, the quality of this asymptotic is rather poor. If we insist on a more precise decay law, we can retain small corrections in (31) and get what turns out, upon a numerical test, to be a reasonably good approximation (see figure 4):

(32)

For the maximally amplified modes (see (29)), the dimensional damping rate is with given by (30). This tells us is that the decay is just slightly faster than exponential at the rate of order . The longest-wavelength modes decay the slowest, after having being amplified the longest.

4 Solution including ITG

Let us now generalise the results obtained in § 3 to include non-zero (i.e., non-negligible) density and temperature gradients. This means that we restore the terms involving and in the general integral equation (11).

4.1 Short-time limit: the ITG-PVG dispersion relation

The short-time limit, introduced at the beginning of § 3.1 for the case of pure PVG, is treated in an analogous fashion for the general ITG-PVG case. An analysis of the solutions of the resulting dispersion relation is useful in that its results assist physical intuition in ways relevant for some of the forthcoming discussion, but it is not strictly necessary for us to have them in order to work out how transient growth happens in the presence of the ITG drive. We have therefore relegated this analysis to B.2.

4.2 Long-time limit

We now continue in the same vein as in § 3.2 and consider the long-time limit (, ), in which we can simplify the kernels involving the Bessel functions in (11) by using (19) and also

(33)

This allows us to rewrite (11) in the form that generalises (20):

(34)

As in § 3.2, we seek a solutions to this equation in the form (21), taking and expanding the delayed potential under the integral according to (22). The result is the generalised form of (23):

(35)

Using again the plasma dispersion function (15) to express the time integrals and introducing the complex scaled frequency , we get

(36)

a time-dependent dispersion relation, which is the generalisation of (24).

Transient growth.

The general argument that the real part of (i.e., the effective time-dependent growth rate) must eventually decrease and so fluctuations will, in the end, decay, applies to (35) similarly to the way it did to (23) (see § 3.2.1), although this decay need not (and, as we will see, will not) be monotonic. The time when the transient growth ends is now determined as follows. Let and (real!). Then from (36) taken at , we find the real frequency by demanding that the imaginary part of the right-hand side vanish — this means that the coefficient in front of must be zero, because, for real , . This condition gives