A Measuring skewed distributions

Symmetry breaking in MAST plasma turbulence due to toroidal flow shear


The flow shear associated with the differential toroidal rotation of tokamak plasmas breaks an underlying symmetry of the turbulent fluctuations imposed by the up-down symmetry of the magnetic equilibrium. Using experimental Beam-Emission-Spectroscopy (BES) measurements and gyrokinetic simulations, this symmetry breaking in ion-scale turbulence in MAST is shown to manifest itself as a tilt of the spatial correlation function and a finite skew in the distribution of the fluctuating density field. The tilt is a statistical expression of the “shearing” of the turbulent structures by the mean flow. The skewness of the distribution is related to the emergence of long-lived density structures in sheared, near-marginal plasma turbulence. The extent to which these effects are pronounced is argued (with the aid of the simulations) to depend on the distance from the nonlinear stability threshold. Away from the threshold, the symmetry is effectively restored.

Keywords: Flow shear, tokamak turbulence, Beam Emission Spectroscopy, gyrokinetic simulations

1 Introduction

Mean flow shear associated with toroidal differential rotation is believed to play an important role in the suppression of turbulence and of the resulting transport in tokamak plasmas [1, 2, 3], most notably in the enhanced confinement regimes of the H-mode [4, 5] and Internal Transport Barriers (ITBs) [6, 7, 8]. A previous study of MAST plasmas [9] found that there was a significant correlation between higher flow shear and higher temperature gradients, suggesting reduced transport; the same trend has been reported in numerical simulations [10]. However, pinning down experimentally a precise manifestation of the effect of flow shear on the local properties of the turbulence, which, presumably, causes the transport, proved to be difficult [11] and remains so. In this paper, we show, both experimentally and numerically, that the flow shear does significantly alter the properties of plasma turbulence. Besides being an essential step towards learning why and how it does so and, therefore, how turbulent transport could be manipulated in practice, it is also of fundamental physical interest to quantify and understand how the statistical properties of a turbulent plasma respond to flow shear.

Let us start by observing that, in the absence of flow shear, the local dynamics of perturbations to an underlying Maxwellian equilibrium is constrained by a type of reflection symmetry, which follows from the up-down symmetry of the magnetic configuration [12, 13]. In the presence of flow shear, this symmetry will be broken and we shall see that the signature effects of the shear can be understood in terms of this symmetry breaking. Let us explain this more quantitatively.

Tokamak plasmas’ local departures from equilibrium are described by the gyrokinetic theory (for a recent review, see [14]). The gyrokinetic equation determines the evolution in time, , of the perturbed distribution function of the gyrating particles, , where is the (“radial”) coordinate perpendicular to a flux surface, is the coordinate along a field line, is the “binormal” (or “poloidal”) coordinate, and () are the velocity components perpendicular and parallel to the magnetic field; the third velocity-space coordinate, the gyrophase angle, is averaged over in gyrokinetics. In an up-down symmetric equilibrium, with no toroidal rotation, the gyrokinetic equation is invariant with respect to the transformations [12, 13]:


where is the fluctuating part of the electrostatic potential, is the perturbed parallel magnetic potential, and is the perturbed parallel magnetic field.

1.1 Symmetry of field distribution

If a stochastic (turbulent) local state of a tokamak plasma is subject to the symmetry (1), an immediate consequence is that positive and negative amplitudes of , and, therefore, of the density perturbation (which is the most experimentally accessible fluctuating field), must occur with the same probability, i.e., the distribution of must be symmetric (even with respect to positive and negative ). A finite flow shear breaks the symmetry of the gyrokinetic equation under the radial reflection , allowing this distribution to become skewed towards either over- or underdensities.

In what follows, we will show both experimentally (Section 2.4) and numerically (Section 3.2) that this is indeed what happens, although how pronounced the skew is depends very strongly on the distance from the nonlinear stability threshold, with symmetry essentially restored far from it (Section 3.3). Our experimental proxy for the density field is the fluctuating intensity field measured using the Beam Emission Spectroscopy (BES) diagnostic [15] on the Mega-Amp Spherical Tokamak (MAST) [16, 17], while numerically the density field is computed using gyrokinetic simulations of a MAST plasma [18] that employ the local flux-tube code GS2.

1.2 Symmetry of correlation function

Another consequence of the symmetry (1), this time relating to the spatial properties of the turbulent fluctuations, is that the two-point spatial correlation function of (for example) the fluctating density field in the plane must be symmetric (even) with respect to reflection in the radial coordinate . This symmetry too is broken when flow shear is present: the correlation function develops a finite tilt.

Indeed, consider the convective time derivative of the distribution function, which incorporates toroidal rotation (in the gyrokinetic equation, this derivative is equal to the evolution operator incorporating various drifts, nonlinearities and drive terms [14]):


where the toroidal rotation is locally approximated by , we assume that we are in a frame rotating with the flux surface (), and is the flow shear (see the Appendix of [18] for a concise description of this term in the local flux-tube gyrokinetic equation). Under the transformations (1), the two terms in (2) have opposite signs, hence the breaking of the reflection symmetry.

The effect of the shear term in (2) on the properties of turbulence is best elicited by making the coordinate transformation into the “shearing frame” (e.g., [19, 20, 21])


under which the convective derivative (2) becomes simply . It is then legitimate1 to represent the perturbed distribution (and all other fluctuating fields) as a Fourier sum


whence, by comparing the exponents, the “Eulerian” wave number can be expressed as a function of the “Lagrangian” wave number and time as


Thus, in the presence of flow shear, the radial wave number grows secularly in time. If a typical perturbation has a life time , then, using (5), we can estimate the typical value of the radial wave number of turbulent fluctuations to be


This tendency for perturbations with a given to have a certain non-zero (whose relative sign to is set by the sign of the shear ) manifests itself as a tilt in the two-point correlation function, with an angle


The presence of a non-zero tilt is a spatial manifestation of the symmetry breaking due to flow shear.

Tilting of the spatial correlation function by flow shear has been observed before, mainly in the edge of tokamak plasmas, with reports of the elongation [22] and the breaking apart [23] of turbulent perturbations. Similar observations have also been made in linear devices [24]. In the core region of a plasma, on the conventional-aspect-ratio tokamak DIII-D, a small tilt () of the spatial correlation function has been reported and attributed to flow shear [25]. Below, we will show, again both experimentally (Section 2.3) and numerically (Section 3.1), that flow shear on MAST can tilt the correlation function to large angles , making the symmetry breaking very clearly manifest and identifying the tilt as a clear signature of the presence of flow shear.

If , the life time will presumably be independent of and so


will be proportional to the flow shear. Generally, at flow shears that are dynamically significant, will have a nontrivial dependence on . Various theories can be, and have been [26, 27, 20], constructed that provide predictions for the dependence of the life time on the flow shear. Testing such models experimentally and numerically can in principle be done by measuring the tilt and then, by inverting (6), estimating the life time of the turbulence. We will give an example of such an analysis in Sections 2.3.4, 3.1.1, and 3.3.1, although it will fall short of validating any specific theory due to a certain deficit of currently available experimental data. Numerically, we will discover that, like skewness, the life time of the perturbations and, therefore, the tilt angle are strong functions of the distance to the threshold, with the life time becoming shorter and the tilt gentler as this distance is increased (Section 3.3).

The structure of the rest of this article is as follows. In Section 2, we consider three experimental cases that illustrate the effect of flow shear. In Section 3, the same analysis as has been performed on the experimental data in Section 2 is repeated on the turbulent density field from gyrokinetic simulations, the results of which are used to gain a more detailed (and more physical) understanding, through studies varying both the flow shear and the ion-temperature gradient systematically, how these two parameters affect the skewness of the distribution of and the tilt of its correlation function. Finally, in Section 4, we conclude with a summary of our findings, a discussion of their implications and of possible directions for future work. Supplementary technical information can be found in Appendix A, where we discuss spurious sources of skewness in the experimentally measured distributions of the fluctuating intensity field, as well as the effect of the MAST BES instrument functions on this distribution. A similar, very extensive analysis regarding the tilt and other parameters of the spatial correlation function, can be found in [28].

2 Symmetry breaking in experimentally measured plasma turbulence in MAST

2.1 The MAST BES system

The BES system on MAST [15] consists of a radial-poloidal array of avalanche-photodiode (APD) detector channels [29], which image the South-South (SS) heating beam;2 the separation between the detectors is approximately in either direction. The detectors register Doppler-shifted emission, digitised at a frequency of , from collisionally excited neutral-beam atoms. The fluctuating part of the intensity of this emission is proportional to the local fluctuating number density of the plasma (at fixed mean density) [30], thus enabling measurement of a turbulent field.

The spatial localisation of the BES measurements within the beam line allows us to construct the spatial correlation function of the fluctuating density field. When this spatial correlation function is calculated (see Section 2.3.2), we assume that the turbulence is homogeneous. Since there is considerable variation of some of the equilibrium quantities across the full radial extent of the BES array (see Figure 1), we do not use the full array, but, for the entirety of this analysis, only consider a subarray of radial-poloidal channels, which is centred at the major radii given in Table 1.

The BES diagnostic is sensitive both to the local turbulent density field and to global MHD activity [11]. Since we are only interested in the former, we always select for analysis time intervals during which the magnetic signal, measured using a Mirnov coil at the outboard mid-plane of MAST [31], is below a certain threshold. The time windows that we analysed for the three representative cases described in Section 2.2 are: for the BLM case, ; for the DLM case, ; for the IFS case, . These time windows are long () compared to the correlation time of the turbulence in all of these cases (see Table 2). As there are also detector channels in a subarray, this means that a statistically large number of values of the turbulent density field is sampled.

2.2 Experimental example cases

Because the detector separation of the BES system on MAST is quite close to the typical radial correlation length of the turbulence, it is quite hard to find long intervals of BES data that would be sufficiently resolved spatially in order for reliable two-dimensional correlation functions to be obtainable (how this is done and how it is determined whether resolution limits are crossed is explained in [28]). There is, therefore, relatively little data available and so it is not currently feasible to have an extensive, fully resolved parameter scan of MAST turbulence in a broad range of values of flow shear.3 Thus, with what we have, a practical strategy for studying the effect of flow shear on turbulence is to identify a few representative cases and compare them. Namely, we would like to examine cases where there manifestly is or is not flow shear present and ask whether they are different in a clearly measurable and qualitatively understandable way.

Figure 1: Cubic-spline-fitted equilibrium profiles for the BLM, DLM and IFS cases (described in Section 2.2): (a) electron density, measured using Thomson scattering (TS) [32, 33], (b) ion temperature, measured using Charge-Exchange Recombination Spectroscopy (CXRS) [34], (c) electron temperature (from TS), (d) toroidal velocity (from CXRS), and (e) -profile, reconstructed by EFIT constrained by the motional-Stark-effect (MSE) [35] measurements. The vertical solid lines indicate, for each of the three cases, the centre of the BES viewing location for the subarray being used (the width of the subarray is ). The vertical dashed lines indicate the position of the magnetic axis in each case. The black curve in the inset to (d) is the cubic-spline fit used for the calculation of (9).
Quantity Symbol BLM DLM IFS EGK
Shot number Shot
Time into shot
Major radius at centre of BES view
Normalised flow shear (9)
Electron-temperature gradient
Ion-temperature gradient
Electron-density gradient
Electron-ion temperature ratio
Minor radius/minor radius of LCFS
Safety factor
Magnetic shear
Ion-ion collisionality
Electron-ion collisionality
Derivative of elongation
Derivative of triangularity
Shafranov shift
Mach number
Ion gyroradius
Time normalisation
Pitch angle of magnetic field
Ion temperature
Electron temperature
Electron density
Toroidal velocity (from CXRS)
Toroidal velocity (from BES)
Table 1: Local equilibrium parameters for each of the experimental cases analysed (Section 2.2). The BLM and DLM cases are the main comparison cases discussed in Sections 2.3 and 2.4; the IFS case is first introduced in Section 2.3.5. The EGK case gives the parameters for the gyrokinetic simulations analysed in Section 3 and is also discussed there. Profiles of , , , , and are plotted in Figure 1 for the BLM, DLM and IFS cases. The numerical values for these three cases are average quantities over the width of the BES subarray, whilst the values for the EGK case are taken at the centre of the BES array (those were the values used in simulations).

The answer is that they are and to demonstrate this, we compare measurements made using the BES system on MAST from two times in the shot #28155. This shot is particularly useful for distilling the effect of flow shear because at a locked mode occurs, braking the rotation of the plasma [36]. The two times that we consider are before the locked mode (BLM) at and during the locked mode (DLM) at . Measurements of the toroidal velocity using the charge-exchange recombination spectroscopy (CXRS) diagnostic are plotted in Figure 1(d) for these two times. Before the locked mode, an ITB is present at the major radius , resulting in a steep gradient in the toroidal velocity, however, the BES measurement location in that experiment is outside of the ITB at , where the flow shear is less strong. During the locked mode, the ITB no longer exists and the toroidal rotation profile flattens. The BES view is centred at for this latter analysis. In the centre of the BES viewing location, a magnetic island forms, within which there is no gradient in the toroidal rotation: see the inset of Figure 1(d). Therefore, by studying this shot, we have access to two clearly different regimes: one with (BLM) and one without (DLM) the flow shear.

The BLM case is quite typical of the sheared cases that can be found in the BES measurements data base—and most of this data base consists of sheared cases. Indeed, the BES diagnostic requires the SS neutral beam in order to measure the turbulence in MAST and the tangentially directed beam provides a torque to the plasma, causing it to rotate. Therefore, there is almost always some level of flow shear acting on the turbulent fluctuations measured by the BES. The DLM case is, thus, rather special, but that is the point: turbulence in MAST is rarely without a flow shear and so any such situation provides a very useful opportunity for a clean comparison. A third experimental example, which will be introduced in Section 2.3.5, is an intermediate-flow-shear (IFS) case, taken from shot #27278 at , with the BES measurement taken at . It represents a less clear-cut situation, which we include to highlight the range of behaviour that one encounters in MAST. While we do not claim to have identified all of the interesting cases, we did analyse a large number of intervals and the examples we have chosen do appear to us to be a worthwhile selection.

The profiles of the electron number density , ion () and electron () temperatures and the safety factor , in addition to the toroidal velocity already discussed above, are plotted in Figure 1 for all three cases. In Table 1, a full set of equilibrium parameters is provided: the values quoted there are averages over the width of the BES subarray that we use for measurements (see Section 2.1). All three cases are plasmas with double-null-divertor (DND) geometries, and, therefore, the equilibria must be up-down symmetric, making the symmetry considerations of Section 1 relevant.

Together, the three cases described above provide examples across a typical (for MAST) range of values of the normalised flow shear


where is the safety factor, is the distance from the magnetic axis to the BES viewing location, is the ion thermal speed, the ion temperature, the Deuteron mass, the minor radius of the last-closed flux surface (LCFS), the toroidal rotation frequency, the toroidal velocity and the major radius at the BES viewing location (these quantities can be found in Table 1). The definition (9) is the same as that used as an input parameter for the gyrokinetic simulations using GS2, however, it is different from the physical flow shear by a factor of , where is the total magnitude of the magnetic field, and its poloidal component. In the large-aspect-ratio limit, this factor is close to unity, but, for a spherical tokamak such as MAST, there is a significant difference. For consistency with the standard definition of used in simulations [18] (Section 3), we will always use (9) as our flow-shear parameter.

Note that, from inspecting the raw velocity profile in Figure 3(d), we know that the flow shear in the DLM case is approximately zero across half of the BES subarray, however, the value for the DLM flow shear given in Table 1 is nearly two thirds of the flow shear for the IFS case. This apparent discrepancy occurs because a cubic-spline fit (see inset of Figure 1(d)) is performed in order to calculate the radial derivative of the toroidal velocity, and the spline fit underestimates the flatness (or, equivalently, the width) of the region due to the magnetic island. The resulting value of the flow shear, calculated as the average over the BES subarray, depends sensitively on the point at which the spline fit determines the minimum to occur, as the sign of the flow shear changes at this point, resulting in significant uncertainty in the calculated value of the flow shear in this case.

2.3 Two-dimensional spatial structure of the turbulence

Instantaneous turbulent density field

Figure 2: Snapshots of the raw (but 2D-spline-interpolated) BES fluctuating-intensity signal from shot #28155 at: (a) , (b) , (c) , (d) . Times (a) and (b) occur during the BLM time period (significant flow shear) and times (c) and (d) occur during the DLM time period (no flow shear). The correlation functions of these two cases are given in Figure 3, which include these snapshots in the temporal average.

Before discussing statistical properties of turbulence, let us inspect the instantaneous spatial structure of the fluctuating density field during the time windows being studied. In Figure 2, example snapshots are plotted of the BES-measured intensity amplitude at certain chosen times. Each channel has been normalised by its own mean, so these are essentially snapshots of the relative density perturbation .

While we know (and it will also become very obvious in Section 2.3.2) that the BLM case [Figures 2(a) and (b)] has flow shear and the DLM case [Figures 2(c) and (d)] does not, it is not easy (indeed, well-nigh impossible) to deduce that by inspection of such snapshots or by watching the full time-evolution sequences for these two cases. Indeed, comparing Figures 2(a) and (c), we see that similar untilted (and so seemingly unsheared) structures occur in both the BLM and the DLM cases; conversely, comparing Figures 2(b) and (d), we also see that similar tilted structures also occur in both cases.4 Therefore, the morphology of particular realisations of the turbulence does not reveal in any especially glaring way that turbulence is different between the two cases. Statistically, however, it is very different, as we are about to discover.

Spatial correlation function

We now characterise the average spatial structure of turbulence using the spatial correlation function of the BES-measured intensity signal. Denoting by the time-dependent intensity measured by the BES detector channel , we split this signal into mean (averaged over time) and fluctuating parts,


The two-point correlation matrix of the fluctuating intensity field is then


where , and are the radial and poloidal (vertical) locations of the detector channel , respectively, and the angle brackets again denote time averages (the auto-correlation functions for each detector, calculated for , are corrected for photon-noise effects as described in [11] and at the beginning of Section 2.4). Assuming the turbulence is (approximately) spatially homogeneous across the BES subarray used for our analysis (see Section 2.1), its two-point correlation function depends only on the relative distances between channels and . We reconstruct this correlation function from the correlation matrix (11) by considering the relevant ranges of values of and and using the following fitting function


The parameter is used to account for global modes [11], as well as for fluctuations in the neutral beam [37]. As we have selected times during which there is no MHD activity (see Section 2.1), the value of should be small. Indeed, in all measured cases it is at most a few percent, therefore also indicating that neutral-beam fluctuations do not affect the measurements.

The fit (12) allows us to extract the spatial correlation parameters of the turbulence: the radial, , and poloidal, , correlation lengths and the radial, , and poloidal, , wave numbers. Operationally, in order to constrain the fit (12), we fix the product , which is determined from the two-time, single-point correlation function of the intensity field (see Section 2.4 and Appendix A of [28]). To obtain correct values of the spatial correlation parameters of the true density field, one must take account of the finite spatial resolution of the BES system, which is quantified in terms of the point-spread functions (PSFs) of its detector channels. How to do this systematically is worked out in [28]; here , , , and are all corrected for PSF effects using this technique.

Physically, we would like to work with the correlation function in terms of point separations in the plane perpendicular to the mean magnetic field, where is the coordinate perpendicular to the flux surface (so, in the outboard midplane of the tokamak, the same as the radial coordinate) and is the “binormal” coordinate, perpendicular both to the field and to the direction. Assuming that the correlation length of the fluctuating density field along the magnetic field is long compared to , we then convert the correlation function (12) to


where , , , , and . These correlation lengths and wave numbers are given in Table 2. The spatial correlation function (13) is plotted in Figures 3(a) and (b) for the BLM and DLM cases, respectively.

Figure 3: Spatial two-point correlation function (13) of the fluctuating density field for (a) the case before the locked mode (BLM, with flow shear) and (b) the case during the locked mode (DLM, no flow shear), both described in Section 2.2. Table 2 shows the parameters of these correlation functions calculated by fitting (12) to the spatial correlation function (11), correcting for PSF effects, and transforming into the coordinates perpendicular to the magnetic field. These correlation functions should be compared to the correlation functions of the numerically simulated turbulence in Figure 8.
Tilt: /deg
Table 2: Correlation parameters for BLM, DLM, and IFS cases, all described in Section 2.2. The correlation functions (13) with these parameters are shown in Figures 3 and 4. All lengths are normalised by the ion gyroradius ( is the ion thermal speed, the ion cyclotron frequency), all times to ; the values of both of these can be found in Table 1. The rms amplitude of the fluctuating density field is calculated from the rms amplitude of the fluctuating intensity field, as defined by (16), using the method described in Section 7.1 of [28]. The skewness is calculated using (17) for the distribution of . The life time is defined in (14); the correlation time was calculated using the CCTD method [38, 39].

Symmetry breaking

The difference between the BLM case (with flow shear) and the DLM case (without flow shear) is manifest: the BLM correlation function is highly tilted whereas the DLM one is very nearly even in . The tilt can be defined in terms of the tilt angle (7), where we now use the and parameters of the fitting function (13). The values of are given in Table 2. Thus, in line with the theoretical discussion in Section 1.2, the statistical symmetry of the plasma turbulence under radial reflection is broken in the presence of flow shear.

Correlation times and flow shear

As promised in Section 1, we may calculate the effective “life time” of a sheared structure, defined, in view of (6), by


The values of this quantity for our experimental cases are given in Table 2. We also give there the correlation time of the fluctuating density field determined from its time-correlation function using the cross-correlation time delay (CCTD) method [38, 39]. We see that the values of and , while not wildly disparate, cannot really be declared to agree. We are going to be relaxed about this problem: first, because both the CXRS value of and the CCTD value of are subject to large errors (regarding the reliability of the CCTD method in the presence of non-negligible PSFs, see Section 9 of [28]), and secondly, because (6) is obviously a relation allowing an arbitrary constant of order unity. However, if the measurement of is viewed as reliable, the life time defined by (14) can be meaningfully used as an effective measure of the (Lagrangian5) correlation time for comparative studies of turbulence at different equilibrium-parameter values, as it will be in Sections 3.1.1 and 3.3.1.

Correlation functions of MAST turbulence generically exhibit a tilt when flow shear is present, with a range of values of —we have checked this in a number of experimental cases, not shown here. This is a reliable feature, to the extent that it could be thought of as a way to find out whether a flow shear is present.6 Given that the values of might not be computed very precisely from CXRS measurements of the rotation profiles, one might wonder whether in fact a more reliable (and dynamically relevant) measure of the local flow shear would be obtained from the tilt angle of the spatial correlation function and the correlation time computed from the time-correlation function, viz., (ignoring here the difference between Lagrangian and Eulerian correlation times; note that could also be the local shear associated with zonal flows, as will be discussed in Section 3.1.2).

Case of intermediate flow shear

Figure 4: Same as Figure 3, but for the IFS case, whose normalised flow shear (9) is smaller than in the BLM case—the tilt of the correlation function is also smaller (see Table 2).

The measurement uncertainties discussed above and a relatively limited range of values of for which well-resolved BES data is available7 preclude us from being able to provide here a straightforward vs.  parameter study. Generally speaking, larger will give larger , but the spread is large (in Section 3.3, we will argue that this is because is a function not just of but also of how far from the stability threshold the system is). The case of relatively low, but measurable, values of shear is perhaps an interesting benchmark as in this case we expect the tilt to be small. We have found such a case, case IFS, described in Section 2.2 and documented in Tables 1 and 2. This has normalised flow shear (9) that is smaller than that in the BLM case (with its other equilibrium parameters not differing in any particularly remarkable way from the BLM ones).

The correlation function for the IFS case is shown in Figure 4. Its tilt is in between that of the BLM and DLM cases and its is shorter than for the stronger-sheared BLM case (see Table 2). Clearly, two points a parameter scan do not make, but this is the best one can do with current data. The salient message is that the transition from stronger flow shear (BLM) to no flow shear (DLM) is gradual in parameter space, with the tilt angle passing through a range of values—and a more comprehensive parameter study ought to be on experimentalists’ agenda. The most interesting outcome from such a study would be the dependence of on , telling us how the dynamical properties of turbulence change with shear (cf. Section 3.3.1).

2.4 Distribution of fluctuation field

In Section 1.1, we argued that the breaking of the radial-reflection symmetry by the flow shear can lead to the breaking of the symmetry (evenness) of the distribution of the fluctuating density and, therefore, BES-measured intensity field. Here we show that this is indeed the case.

We consider the fluctuating part of the intensity field (10) and, for each detector channel , normalise it by its own standard deviation


where we have subtracted the mean square amplitude of the fluctuating part of the photon noise signal (determined from the signal at each channel during the calibration of the instrument [11]). The angle brackets signify averages over time. The normalisation to is required in order to account for a degree of variation of this quantity across the BES subarray used for our analysis. The total standard deviation (total rms fluctuation amplitude) over all channels in the subarray is


This quantity is proportional to the total rms fluctuating density field , whose value can be reconstructed from it by correcting for PSF effects (see Section 7.1 of [28]) and is given in Table 2. The distribution of the normalised fluctuating intensities can be affected by several different known experimental effects, including MHD activity, radiation spikes, and PSF effects. A discussion of these effects and how they have been accounted for in our present analysis is given in Appendix A.

Symmetry breaking

Figure 5: Distribution of the fluctuating intensity field for (a) the case before the locked mode (BLM, with flow shear) and (b) the case during the locked mode (DLM, no flow shear), both described in Section 2.2. For each detector channel, is normalised by its own root-mean-square (standard deviation) value (15). The black-dashed line in each case gives the unit normal distribution. The distribution function and its skewness for the DLM case (b) were calculated from 2 (rather than 5, as in all other cases) sets of 4 poloidal BES channels, located at  m and  m, where the toroidal rotation profile is completely flat, as seen in Figure 1(d).

In Figure 5, we show the probability distribution of the normalised fluctuating intensities for the BLM and DLM cases. In the case without flow shear (DLM), the distribution is even with respect to positive and negative (and very nearly normal). In contrast, the case with shear (BLM) exhibits a relatively small, but measurable preponderance of positive , i.e., overdensities are statistically somewhat more common than underdensities. This can be quantified as a positive value of the skewness of the distribution8


For the BLM case, ; for the DLM case, , but this is similar to the skewness in the background emission signal, as well as to the size of the error in the determination of the skewness that occurs due to PSF effects (see Appendix A.2), so cannot be considered significant.

Figure 6: Same as Figure 5, but for the IFS case (lower flow shear than for BLM; see Table 2).

In Figure 6, we show the fluctuating-intensity distribution for the IFS case, with flow shear lower than BLM, already discussed in the context of its lower tilt in Section 2.3.5. We see that the distribution is again skewed, with its skewness lower than for the BLM case (but still measurable), just as its tilt was (to identify a quantitative dependence between the skewness and flow shear, we will need to resort to simulations, in Section 3.2, as we do not have a sufficient range of data for a proper parameter scan). What is more interesting is that the sign of the skew is negative, i.e., here, underdensities are privileged compared to overdensities. This emphasises the fact that, while the presence of a flow shear allows such asymmetric distributions, theory (at least in its current state) does not predict the sign of the asymmetry—we do not know how the sign of the skew depends on local equilibrium parameters.

Finally, we must stress a very strong caveat. The fact that the symmetry (1) is formally broken by flow shear does not necessarily imply that a skewed distribution should result, it merely allows it (in contrast to the case of zero flow shear, when a skewed distribution would clash with theory [12, 13]). Here we have demonstrated, as a proof of principle, that skewed distributions can indeed be found experimentally when flow shear is present. However, skewness is not as robust and inevitable effect of flow shear as tilted correlation functions appear to be and we have found many experimental cases in which flow shear is present according to CXRS measurements but the skewness of the fluctuating-intensity distribution is so small that, in view of measurement uncertainties, it cannot be declared present beyond reasonable doubt. In Section 3.3, we will see, by analysing numerical simulations, that how skewed the distribution is depends quite strongly on how far from the nonlinear stability threshold the system finds itself, with only near-marginal cases showing significant skewness. Thus, while a tilt (possibly small) in the correlation function might be viewed as a signature of the presence of flow shear (as argued in Section 2.3.4), a noticeable skewness of the fluctuating-density distribution should perhaps be viewed as a signature of a sheared turbulent state being close to threshold. The physical reason for this is probably that, close to the threshold, the turbulent state in MAST is dominated (as shown in [18] and discussed in Section 3.2.1) by a finite number of long-lived, finite-amplitude coherent structures, whereas farther from the threshold, a more (evenly) distributed sea of fluctuations emerges.

Skewness vs. tilt: two-component turbulence?

Figure 7: Same correlation function as in Figure 3(a), for the BLM case, but this time conditioned on (a) lower intensities (core of the distribution, ) and (b) higher intensities (tail of the distribution, ), in the way described in Section 2.4.2. The parameters of the fitting function (13) in each case are given in Table 3. The tail is manifestly less tilted than the core and also has a larger radial correlation length. These correlation functions should be compared to the ones for numerically simulated turbulence in Figure 13.

The above discussion and the greater robustness of the tilted correlation functions than of the skewed distributions as signatures of flow shear suggest that these features describe the response to the flow shear of different types of turbulent fluctuations. Indeed, intuitively, as perturbations are sheared by the flow, their radial wave number increases according to (5) and so then does . This will push them into a damped region of wave-number space (see, e.g., [20]), quickening their demise. It is these perturbations, which are on their way out, that would statistically contribute to the tilt in the correlation function. In contrast, one might imagine that rare but particularly strong perturbations whose presence skews the distribution of the fluctuating field, might be less vulnerable to (or simply not yet affected by) such destruction by shear and so would turn out to be less tilted. This view is perhaps reinforced by the observation, which we will be in a position to make when we consider the numerically simulated turbulence in Section 3.2, that the rms fluctuation level is significantly lower in the turbulence with flow shear than without it and so the tail of the distribution represents not some additional population appearing alongside the “core turbulence” but rather the “survivors”—sufficiently large structures that can create sufficiently large local velocity gradients to overcome the suppressing effect of the shear; these large structures would diffuse and disappear rapidly in the strong unsheared turbulence.

These speculations can be experimentally tested, if in a rather tentative manner, by constructing a correlation function conditional on the fluctuation amplitude. Indeed, the turbulent density field is continuously (and quickly) advected by the flow past our, relatively small, BES subarray, which can thus be thought of as sampling different types of turbulent structures at different times, rather than capturing a fully statistically representative spatial domain containing all types of structures at all times. Thus, we collect all snapshots of the turbulent field (such as those discussed in Section 2.3.1) in which the fluctuating intensities at all detector channels are below a certain threshold, specifically, , and calculate the spatial correlation function of this collection, representing the “core” of the distribution, in the way described in Section 2.3.2. We repeat the same procedure for a complementary set of snapshots, in which at least one detector channel sees , to extract the correlation parameters of the “tail” of the distribution. This analysis is necessarily crude, because the reduced number of snapshots (especially for the tail) reduces the quality of the statistics and the size of the BES window may have an effect that we cannot study systematically. Proceeding with it nevertheless, we find that the tilt angle for the core of the distribution in the BLM case is noticeably larger than in the tail (and also that the radial correlation length in the tail is much larger than in the core; see Table 3). For the DLM case, there is no significant difference between the tail and the core, the tilt angle is essentially zero for both. This result appears to support the line of reasoning proposed above, although obviously a much more extensive study of a large number of well-resolved cases with good statistics is called for.

The need for such a study to make our argument fully convincing (or otherwise) is accentuated by the disconcerting failure of the IFS case to agree with the core vs. tail distinction seen in the BLM case: in fact, we find that the tilt angle in the tail is larger than in the core (see Table 3). We do not know why the IFS case has this behaviour: whether, for example, it might be related to the lower (dynamically unimportant?) flow shear, or to the preponderance of under-, rather than overdensities in the distribution of the fluctuating field, or to the presence of zonal flows.9 In Section 3.2.2, we will see that the trend in numerical simulations appears to be broadly consistent with the BLM case.

Case BLM core BLM tail IFS core IFS tail
No. of snapshots 1434 123 3065 206
Tilt: /deg
Table 3: Correlation parameters of the fitting function (13) for the BLM and IFS core and tail conditional correlation functions (for the BLM case, they are shown in Figure 7). The errors given here are for the fitting procedure only. The parameters for the overall correlation functions can be found in Table 2. Note that the number of snapshots does not represent the number of statistically independent instances: neighbouring snapshots are separated by , so it takes approximately 20 snapshots for a given structure to pass through the BES subarray (due to the apparent poloidal motion caused by toroidal rotation [39]).

Cognizant of all the uncertainties and of the clearly insufficient amount of experimental evidence on offer, we conclude this section with a cautious conjecture that turbulence under the action of flow shear consists of a tilted, sheared, statistically more ubiquitous, but smaller-amplitude component and less tilted, rarer, but stronger perturbations, which skew the distribution function of the fluctuating field. We shall test this idea further in the next section.

3 Symmetry breaking in numerically simulated gyrokinetic turbulence

Having established experimentally, at least as a proof of principle, that the presence of the flow shear breaks both the spatial symmetry of the turbulence and the symmetry of the distribution of the fluctuating density perturbations, we now turn to gyrokinetic simulations, both to check that modelling can reproduce what theory predicts and the experiment has revealed and to seek clearer guidance on the parameter dependence of these effects. These numerical results give us a degree of confidence in the validity of our interpretation of the experimental data, which, especially in what concerns the skewed distributions (Section 2.4), was perhaps not on its own sufficient to build a rock-solid case.

In this section, we use a set of local, flux-tube, electrostatic, two-species, ion-scale numerical simulations using the gyrokinetic code GS210 that were originally performed for a different, independent study of the transition to turbulence in MAST [18]. The experimental case whose equilibrium parameters served as input for these simulations has previously been used for validation and verification exercises involving both NEMORB [40] and GS2 [41]. This experimental case is labelled EGK and its parameters are given in Table 1. While the BES data from this case has been studied extensively using 1D reconstruction of various correlation parameters [40], it unfortunately turns out not to be suitable for a full-scale 2D analysis as developed in [28], because calculation of and, therefore, the tilt, suffers from insufficient resolution. However, it is not wildly different from our BLM case11 and, in any event, our aim here is to determine qualitative trends rather than to claim precise agreement between simulations and experiment (although a more limited 1D analysis of the BES data for the EGK case does show respectable agreement with GS2 results [41]).

The obvious advantage accorded by (local flux-tube) simulations is that one can change individual equilibrium parameters at will, while keeping other parameters constant, and thus look for trends. We will take this route by analysing a sequence of simulations in which the flow shear is varied compared to the nominal EGK case—and confirm that the gyrokinetic simulations are able to recover, qualitatively, both the tilted correlation functions (Section 3.1) and the skewed distributions (Section 3.2) observed in the experiment.

Run name GKa5 GKa4 GKa3 GKa2 GKa1 Marginal
Table 4: Various statistical characteristics of turbulence in a series of simulations corresponding to the equilibrium parameters of the EGK case (Table 1, except in these runs) and a sequence of values of flow shear . The Marginal case is the EGK case itself. The rms fluctuation amplitude is calculated using (16), but replacing with ; the summation is over all grid points . Similarly, the skewness is calculated for the fluctuating density field using (17). The spatial correlation functions for the GKa5, GKa3, GKa1, and Marginal cases are plotted in Figure 8. The fluctuating-density distributions for the GKa5, GKa1, and Marginal runs (the latter two skewed) are shown in Figure 12. The conditional correlation functions (CORE and TAIL) are discussed in Section 3.2.2 and plotted in Figures 13 (run GKa1) and 14 (Marginal run).

Figure 8: Directly calculated spatial correlation functions of the fluctuating density field for 4 of the runs documented in Table 4, with, from (a) to (d), value of the flow shear decreasing from the experimental (EGK) value to . The spatial domain of the simulation is regularly gridded, therefore multiple values of the correlation function (11) (with instead of , being the grid point) occur for each pair; these values are then averaged to produce the spatial correlation function that is plotted (it is also averaged over time, typically several hundred , i.e., tens of correlation times). These correlation functions are to be compared with experimental correlation functions with and without flow shear in Figure 3. The correlation parameters for the fitting function (13), approximating the true correlation functions plotted here, are given in Table 4. The red contours in these plots correspond to the fitting function (13) with these parameters, showing the quality of the fit (and thus supporting its use for experimental data; see [28] for further discussion).
(a) vs.  (b) vs. 
Figure 9: (a) Life time of perturbations, defined by (14), , and normalised by , vs. flow shear for a number of values of the ion-temperature gradient (distinguished by different shapes of the data points); the data points are coloured according to the value of the gyro-Bohm-normalised heat flux (see discussion in Section 3.3). (b) Same, but here is plotted vs.  with data points coloured according to the value of (this plot includes data from the full set of simulations carried out in [18] and so covers a larger number of values of ITG than (a)); the vertical dashed line indicates the experimental value of for the EGK case (Table 1).

However, it must be clear that this does not, in fact, constitute a parameter scan that one might expect to achieve experimentally in a real plasma. Indeed, the suppression of the turbulence caused by increasing the flow shear will cause the plasma equilibrium to respond: as the turbulent transport is reduced, assuming constant ion heat flux, the ion-temperature gradient (ITG) must increase [10, 9]. Generally speaking, experimentally one expects to find turbulent plasma states near the stability threshold (which, in the instances of ion-scale MAST turbulence that we are considering is a nonlinear stability threshold, the turbulence being subcritical [18]). As a crude simplification, we can think of this threshold (or “zero-turbulence manifold” [10]) as a boundary in a two-dimensional parameter space consisting of and ITG, although obviously other equilibrium parameters are in fact also involved [19, 10, 9, 42]. From the point of view of interpreting experimental results, the key questions are (i) near which part of the zero-turbulence manifold the experimental case under consideration is located and (ii) how close to the threshold it is. The second question is very crucial indeed, because we know from simulations [18] (and will further confirm in what follows) that the properties of our turbulence depend on the distance from threshold very sensitively. To investigate these matters, in Section 3.3, we will consider a set of numerical simulations in a range of values of both flow shear and ITG and argue that both the tilt and the skewness are strong functions of the distance from threshold, i.e., that the symmetry-breaking effects associated with the flow shear fade away and the symmetry is restored if one probes deep into the strongly turbulent regime far from the stability boundary.

3.1 Tilt vs. flow shear

We start by considering a parameter scan in flow shear, with 6 different values of (see Table 4) and all other equilibrium parameters the same as those of the EGK experimental case, listed in Table 1 (except, in these runs, , which is within the experimental uncertainty of the measurement [18]). The largest value of the flow shear in this scan, , corresponds to the EGK case itself—and at even larger values, there is no turbulence (all perturbations decay), so we will refer to the simulated EGK case as Marginal.

For the purposes of this analysis, we adopt the definition of the fluctuating density field analogous to the definition (10) of the fluctuating intensity in the experimental measurements reported in Section 2: at each grid point , it is the density perturbation at this point (extracted directly from the simulations—it is the velocity integral of the perturbed particle distribution function, which is being solved for by GS2), with a running time average, calculated over a moving interval, subtracted, so the time average of the fluctuating field is zero by definition. Note that this corresponds to using a  kHz high-pass frequency filter (which experimentally was needed anyway; see Appendix A.1) and that this procedure will remove any long-lived zonal density component (see Section 3.1.2).

The spatial correlation functions of the fluctuating density fields calculated in this fashion are shown in Figure 8 and the correlation parameters (lengths and wave numbers) defined by the fitting function (13)12 are documented in Table 4. It is manifest that, as the flow shear increases, the radial-reflection symmetry is broken and the correlation function develops a tilt, which becomes quite large as the Marginal case is approached. Qualitatively, the difference between the zero-flow-shear, untilted correlation function for run GKa5 and the sheared, tilted ones for runs GKa3 or GKa1 is very similar to the difference between the DLM and BLM cases shown in Figure 3.

Tilt and life time

Quantitatively, we find that the tilt increases approximately linearly with , except near the stability threshold. Recalling the discussion at the end of Section 1.2 and in Section 2.3.4, we conclude that, except near the stability threshold, the life time (14) of the turbulent perturbations is independent of the flow shear—presumably because dynamically the flow shear is not very important away from the threshold. This point is illustrated by Figure 9(a). The data for other values of the ion temperature gradient also shown in this figure supports this conclusion and will be further discussed in Section 3.3.1.

Tilts at zero flow shear as signature of zonal flows

(a) without high-pass filter (b) with high-pass filter
Figure 10: Perturbed density field from run GKa5 () integrated over shown vs.  and time (a) without subtracting moving time average and (b) with moving time average subtracted (as for all fluctuating density fields used in Section 3). The presence of a zonal density component is manifest. The red boxes show regions for which conditional correlation functions shown in Figure 11 and discussed in Section 3.1.2 were calculated.

Figure 11: Same correlation function as Figure 8(d), for run GKa5 (), but calculated for two spatial subdomains: (a) with positive zonal density and (b) with negative zonal density, as discussed in Section 3.1.2.

As we saw in Figure 8(d), when the mean flow shear is zero, the correlation function of the fluctuating density field is untilted, in agreement with the symmetry (1). We stress, however, that this is only true if the correlation function is calculated over a sufficiently large spatial domain. Because drift-wave turbulence is prone to developing zonal flows [43, 44], locally a flow shear can be present—and indeed is, in the simulations that we are analysing. As these simulations have kinetic electrons, they can, and do, develop a zonal density perturbation on ion scales, so the presence of a zonal field is detectable directly from the measured density field, in the form of a -independent, long-lived component. Because zonal fields have a long life time, compared to the life time of the turbulence, they can only be seen if we remove the high-pass frequency filter effectively imposed on our fluctuating density field by subtracting the running time average (as stated at the beginning of Section 3.1). This is illustrated in Figure 10, where we show the -integrated perturbed density field for the zero-flow-shear run GKa5 with and without this filter.

Experimentally, the high-pass frequency filter may be difficult to abandon (see Appendix A.1), however, the zonal shear can be detected indirectly via correlation functions of the fluctuating density field:13 if those are calculated over a spatial domain that is smaller (radially) than the wave length of the zonal flow, they will be tilted by the same mechanism as in the case of a mean flow shear. This is illustrated in Figure 11, which shows correlations functions for the zero-flow-shear run GKa5 restricted to regions of positive and negative zonal density, as marked in Figure 10,—as is obvious from the tilts of the correlation functions, these regions turn out also to be regions of positive and negative zonal shear (in this context, we should perhaps recall the presence of tilted instantaneous structures in the zero-flow-shear DLM case, noted in Section 2.3.1). The overall correlation function in Figure 8(d) is an average over several such regions—the opposite tilts average out and the overall radial-reflection symmetry is preserved.

Since turbulence with mean flow shear will also develop zonal flows, a similar analysis applied to it will show a certain spread in the tilt angle of the correlation function depending on whether the zonal shear locally subtracts from or adds to the mean flow shear. This can matter if the mean flow shear is sufficiently weak.

3.2 Skewness vs. flow shear

Figure 12: Distribution of the fluctuating density field , normalised by its own standard deviation, for (a) GKa5 (no flow shear), (b) GKa1 (with flow shear, close but not at the nonlinear stability threshold), (c) Marginal (at the threshold) runs (see Table 4). The black-dashed line in each case gives the unit normal distribution. These distributions should be compared with experimental ones with and without flow shear in Figure 5.

The numerical scan in flow shear also reveals that the symmetry of the distribution of the fluctuating density field is broken by the flow shear: the skewness of the distribution increases with the flow shear, as documented in Table 4; the distributions for a run with no flow shear (for which in accordance with the theoretical expectations outlined in Section 1, there is no skew) and runs with two values of the flow shear, one at the nonlinear stability threshold and one just off it, are shown in in Figure 12. The skew in the distribution for the Marginal case in Figure 12(c) is extremely pronounced, much more so than anything that we have observed experimentally, but the distribution for the case that is only slightly offset from the stability threshold, shown in Figure 12(b), is very similar qualitatively to the experimental BLM case in Figure 5(a).14 Generally speaking, in many runs with various values of and , we only find distributions with a pronounced tail on the overdensity side in cases that are very close to the threshold (and in all cases that were simulated, we do find asymmetry in favour of over-, rather than underdensities), while the majority of the skewed distributions are similar to that in Figure 12(b). This gives us a degree of confidence that the skewness that we have found in the experimentally measured distributions in Section 2.4 is indeed due to flow shear. It is likely that the fact that this skewness is not very large (and indeed not always measurable) is due to the rather sensitive dependence of it on the distance to the threshold: note the precipitous decline in skewness from the Marginal run to run GKa1, even though the difference in the values of for these two cases ( and , respectively) is small and experimentally would be within measurement errors on the velocity profile (see further discussion in Section 3.3).

Physics of skewed distributions

Besides encouraging the physical interpretation of the experimental evidence in terms of symmetry breaking associated with flow shear, numerical simulations give us a crucial insight as to how, dynamically, this symmetry breaking occurs. In [18], for which these simulations were originally carried out, it was shown that the transition to turbulence in these flow-sheared MAST configurations occurs via emergence and accumulation (as equilibrium parameters move deeper into the unstable regime) of long-lived, intense coherent structures, which occupy a small fraction of the volume. This is because turbulence in these equilibria is subcritical—while the system is formally linearly stable, perturbations do grow transiently and, for finite perturbations, this can lead to non-zero nonlinearly saturated turbulent fluctuation levels and heat fluxes. Since finite amplitudes are required for such a state to persist, the only way for the system to have small overall heat fluxes as the threshold is approached is by restricting finite amplitudes to a small fraction of the volume, hence the spatially sparse, but intense structures. Clearly, in terms of the probability distribution of the amplitudes, the presence of such structures and the fact that they give rise to a non-zero overall heat flux, implies skewed distributions. Away from the threshold, the structures become more numerous, overlap, interact, break up and dissolve into a more volume-filling turbulent sea of “eddies,” characteristic of standard, unsheared drift-wave turbulence—resulting in approximately symmetric distributions.

Another important observation afforded us by the numerical simulations but unavailable experimentally (because equilibrium parameters in experiments cannot in general be changed independently) is that, other things being equal, increasing flow shear does have a suppressing effect on the rms amplitude of the turbulent density field: this is documented by the values in Table 4 (it is not evident in Figure 12, because these are distributions of normalised to its own standard deviation). Thus, as we anticipated in Section 2.4.2, the tail of the skewed distributions found for the near-threshold, sheared turbulence is, in a sense, a remnant of a wider unsheared/far-from-threshold distribution, containing large-amplitude perturbations that might have existed in a more unstable local equilibrium configuration (although, as argued above, these perturbations have a rather distinct character once only a dwindling number of them are left with a sole responsibility of maintaining turbulent transport; see also Section 3.2.2).

Skewness vs. tilt: two-component turbulence in simulations

Figure 13: Same correlation function as Figure 8(b), for run GKa1, but this time conditioned on (a) lower densities (core of the distribution, standard deviations) and (b) higher densities (tail of the distribution, standard deviations) in the way described in Section 3.2.2. The correlation parameters for each case are given in Table 4. These correlation functions can be compared to Figure 7 for the experimental BLM case.

Figure 14: Same as Figure 13, but for the Marginal run (see Table 4 and Figure 8(a)). The bottom row has correlation functions calculated in the standard way from the high-pass-filtered fluctuating density field (as explained at the beginning of Section 3.1); the top row has correlation functions of an unfiltered field, i.e., including the zonal density perturbation. The latter manifestly dominates the core population but not the tail.

In Section 2.4.2, we conjectured that turbulence in the presence of flow shear is a two-component mix of strong but rare perturbations giving rise to the skew in the distribution and a wider sea of weaker ones, which constitute the core of the distribution and are tilted, more or less passively, by the flow shear. In other words, the tilt and the skew represent symmetry breaking by flow shear occurring in different types of perturbations.

To test this picture, we introduced conditional correlation functions, designed to measure the spatial correlation parameters associated with the core and the tail of the distribution. Here we perform the same analysis on the numerically simulated fluctuating density fields, with the same condition that frames containing at least one instance of larger than 3 standard deviations are declared to belong to the tail (for the purposes of this selection, the full artificially large  cm GS2 computational domain is split into  cm subdomains—this is comparable to the  cm BES subarray that we used in Section 2.4.2 and necessary in order for conditioning on the maximum value of within a subdomain to pick out statistically different patches of turbulence). The results are documented in Table 4 and conditional correlation functions for the core and tail of two of the runs with flow shear are shown in Figures 13 (run GKa1) and 14 (Marginal run). We see that fluctuations in the tail of the distribution are always less tilted than those in the core. This is similar to the behaviour exhibited by the experimental BLM case in Section 2.4.2.

Another qualitative similarity between simulated and real sheared turbulence is that the perturbations in the tail of the distribution, in both cases, have longer radial correlation lengths , not just lower tilt (i.e., smaller ). Table 3 shows this for both BLM and IFS cases (even despite the failure of the latter to follow the trend as far as its tilt angle is concerned) and Table 4 confirms the same trend for our numerical runs with varying flow shear; it is also quite obvious from the plots of the conditional correlation functions in Figures 13 and 14.

Comparing these two figures, we note also a sizeable increase of the poloidal correlation length in the core of the distribution as the threshold is approached (see Table 4). Interestingly, in the Marginal run, the core fluctuations are so weak, that the zonal density perturbation (in our definition of the fluctuating density field, it is removed by time averaging/high-pass frequency filtering; see Section 3.1.2) becomes comparable to them—so much so that, as we see in Figure 14, it dominates the conditional core correlation function (although not the total one or the tail) if the high-pass frequency filter is removed (suggesting what is perhaps a promising experimental direction of enquiry).

3.3 Restoration of symmetry far from threshold

In our discussion of both experimental and numerical results so far, we have repeatedly anticipated the argument that what truly matters for the degree to which the symmetry of the turbulence is broken is how far from the (nonlinear) stability threshold the system is. Thus, changing the flow shear in Sections 3.1 and 3.2 led to significant changes in behaviour because this took the system closer to or farther from the marginal case.

Run name Marginal GKb1 GKb2 GKb3 GKb4 GKb5 GKb6 GKb7
Table 5: Same as Table 4, but for and a sequence of values of the ion-temperature gradient .

This notion is easily confirmed if we fix the flow shear at what in the previous sections was the marginal value, , and change instead the ion-temperature gradient (ITG) . A sequence of runs within such a parameter scan is documented in Table 5. The trend that emerges is very clear: as the ITG is increased from the marginal value (), both the tilt angle of the correlation function and the skewness of the distribution of the fluctuating densities diminish (the decline in skewness is particularly strong). Thus, the symmetry is broken less strongly.15 The situation in this sense is rather similar to what happened in the scan (Table 4) with decreasing flow shear. One might argue that, as turbulence becomes stronger away from the threshold (with the distance from the threshold measured either in or in ), the flow shear ceases to be an important factor and the symmetry formally broken by it is effectively restored.

Tilt, flow shear and life time

It is worth discussing a little further how the effective restoration of symmetry happens for the tilt of the correlation functions. We argued at the end of Section 1.2 and again in Section 3.1.1 that the flow shear did not need to be dynamically important in order to induce a tilt in the correlation function. In such a “passive” tilting scenario, the tilt is simply proportional to , with the constant of proportionality being the life time of the perturbations, as per (8) or (14). This life time, defined by (14) and given in Table 5, becomes shorter as the ITG is increased. This is also manifest in Figures 9(a) and (b)—the latter figure is a scatter plot of for a large set of values of and against the gyro-Bohm-normalised ion heat flux , which is a good measure of distance to threshold (it is zero at the threshold and increases steeply away from it [18]). The reason correlation times shorten is that larger ITG implies stronger drive and so higher turbulent amplitudes: “eddies” turn over faster.16 Shorter then implies smaller tilt angle () farther from marginality.

In physical terms, the flow shear ceases to be important if it is smaller than the local, fluctuating velocity shear acting on a turbulent perturbation and caused by its fellow perturbations. This nonlinear fluctuating shear must, in saturation, be comparable to some appropriate linear rate at which the perturbations are replenished by the ITG drive. In this sense, the competition between the life time of the turbulence and the flow shear as the threshold is approached is reminiscent of the so-called Waltz rule [48]: that in order to affect (suppress) the turbulence, the flow shear must be similar to, or greater than, the linear growth rate (we recall, however, that the turbulence that we are dealing with here is subcritical [18] and so formally there is no linear growth rate; what effective quantity from the linear theory should be used instead is a matter of current research interest [20, 10, 41]).

Tilt and skewness vs. distance from threshold

(a) Tilt vs.  and (b) Skewness vs.  and
(c) Tilt vs.  (d) Skewness vs. 
Figure 15: (a) Tilt angle of the correlation function and (b) skewness of the distribution of the fluctuating density field vs. flow shear and ITG in a range of values of these parameters around the nonlinear stability threshold. The tilt and the skewness are replotted in (c) and (d), respectively, vs. the gyro-Bohm-normalised ion heat flux, with data points coloured according to the value of . These plots contain data for and ITG , covering the entire database of the runs carried out in [18]. In (d), inverted triangles mark the cases where skewness is negative.

Figures 15(a) and (b) show the tilt angle of the correlation function and skewness of the distribution of the fluctuating density field for a range of values of the flow shear and ITG around the nonlinear stability threshold and departing from it—this extends the and scans documented in Tables 4 and 5. The conclusion remains the same: away from the threshold, symmetry is effectively restored—particularly quickly for the distribution functions (as we argued in Section 3.2.1, the restoration of symmetry of the distribution of fluctuating densities far from the threshold occurs because intense coherent structures that dominate the transition to turbulence and are responsible for the skewed distributions become too numerous for individual survival, interact and break each other apart [18]).

As we already did in Section 3.3.1 for the life time of the perturbations, we quantify the dependence of the measures of asymmetry (tilt and skewness) on the distance to threshold in terms of their dependence on the gyro-Bohm-normalised heat flux , which vanishes at the threshold but increases strongly and monotonically away from it [18] and thus constitutes a good order parameter. Scatter plots of the tilt and the skewness vs.  are shown in Figures 15(c) and (d), respectively, for the entire data base of numerical simulations originally carried out in [18]. We see that, modulo some scatter, both the tilt and the skewness are approximately functions of the distance to the threshold (to be precise, Figure 9(b) suggests that it is the life time , rather than the tilt, that is a nearly single-valued function of this distance, at least when measured sufficiently far from the threshold). The case of is the exception to this rule, as expected theoretically, being subject to the exact symmetry (1): namely, at zero flow shear, both the tilt and the skew are essentially zero regardless of how marginal, or otherwise, the system is.

4 Discussion

To summarise, we started with a premise that, as the presence of the flow shear appears to be quite strongly correlated with larger temperature gradients both in MAST [9] and in numerical simulations [10, 18], the effect of the shear on the local structure of plasma turbulence must be detectable. Theoretically, we argued (in Section 1) that flow shear would break the reflection symmetry (1) [12, 13] and that statistically this symmetry breaking would be actualised in the form of skewed distributions (Section 1.1) and tilted correlation functions (Section 1.2) of the fluctuating density field. We did then find both of these signatures of symmetry breaking experimentally (Sections 2.4.1 and 2.3.3, respectively), at least as a proof of principle—in specific examples of turbulent density fields in MAST with and without flow shear.

For a number of reasons to do with spatial resolution of the MAST BES system and with the range of equilibrium parameters available, we cannot, at this point, have a broad parameter-scan study and so we sought reassurance, validation and further physical insight from a set of numerical simulations of MAST turbulence performed by [18] and covering a range of flow shears and ion-temperature gradients in a MAST-relevant equilibrium configuration. Again we found both skewed fluctuating-density distributions (Section 3.2) and tilted correlation functions (Section 3.1), qualitatively similar to the experimental ones.

With the aid of the simulations, we were able to establish that, generally speaking, the symmetry breaking was most pronounced in the vicinity of the (nonlinear) stability threshold, whereas away from it, symmetry was gradually restored (Section 3.3.2)—essentially because flow shear became dynamically less and less relevant far from the threshold (Section 3.3.1). The dependence of the symmetry breaking on the distance from threshold—especially of the skewness of the fluctuating-density distributions—turned out to be quite strong. Experimentally, while we expect to find turbulent states in the general vicinity of the stability boundary (as indeed confirmed in [18]), it is not reasonable to expect them to be exactly on it. The degree to which symmetry is found to be broken can then be viewed as an indicator of how close to the threshold any given measured instance of turbulence is.

We have also argued that tilted correlation functions provide a versatile diagnostic that could be used to probe the effective life time (correlation time) of turbulent structures (Sections 2.3.4, 3.1.1, and 3.3.1) and the local zonal shear (Section 3.1.2). Furthermore, when conditioned on the amplitude of the fluctuating density, these correlation functions appeared to reveal, both experimentally (Section 2.4.2) and numerically (Section 3.2.2), that sheared, skewed, near-threshold turbulence could be thought of as consisting of two components: longer-lived, spatially larger, more intense, but statistically rarer structures and a general sea of smaller, weaker, but statistically volume-filling fluctuations. Indeed, that spatially sparse coherent structures dominate the transition to turbulence in the physical situation considered here was shown by [18] (see the summary of their argument in Section 3.2.1)—what we have found here is that the statistical signature of these structures is the breaking of symmetry (skewness) of the distribution function of the fluctuating density, or, conversely, that the physical way in which the system takes advantage of the breaking of symmetry is the emergence of these structures. The tilting of the correlation functions is a better understood and, arguably, more straightforward phenomenon (Section 1.2) affecting both strong and weak perturbations, the latter more than the former (see Sections 2.4.2 and 3.2.2).

It is quite clear that, experimentally, a much more extensive study is called for. The trends revealed by the numerical simulations serve as a guideline to the range of questions that can be asked and of parameter dependences that could be verified or falsified. It would be interesting to learn in what circumstances fluctuating-density distributions are skewed towards underdensities, rather than overdensities (exemplified by one of our experimental cases; see Section 2.4.1). The nature and role of zonal flows generally and in subcritical turbulence with flow shear in particular (such as the MAST turbulence that we have investigated), remain a fascinating subject, which perhaps could be tackled experimentally using some of the tools developed here (see Section 3.1.2) or other methods [45]. More generally, a full experimental (as well as numerical and theoretical) understanding of the transition to turbulence and of the resulting turbulent states in subcritically unstable, sheared tokamak turbulence is a worthwhile goal, important for practical purposes of optimising tokamaks and perhaps learning how to manipulate stability boundaries, but also appealing to the more curiosity-driven angels of our nature.


This work was carried out within the framework of the EUROfusion Consortium. It received funding from the Euratom research and training programme 2014-18 under grant agreement No 633053 and from the RCUK Energy Programme [grant number EP/I501045]. The views and opinions expressed herein do not necessarily reflect those of the European Commission. Financial support was also provided to MFJF by Merton College, Oxford. YcG was supported by the National R&D Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT and Future Planning (Grant No. 2014M1A7A1A01029835) and by the KUSTAR-KAIST Institute, KAIST, Korea. The work of AAS was supported in part by grants from UK STFC and EPSRC. We thank M. Barnes, S. Cowley, W. Dorland, E. Highcock and C. Roach for many inspiring discussions of plasma turbulence and its modelling.

Appendix A Measuring skewed distributions

Experimental measurement of skewed (or otherwise) distributions of the fluctuating intensity fields is complicated by a number of extraneous effects, which we detail and describe how to correct for in Appendix A.1. In Appendix A.2, we extend the work of [28] to determine the effect of PSFs on the skewness of the fluctuating-field distributions.

Figure 16: Probability distribution of the intensity perturbation (unnormalised) for (a) the DLM case (see Section 2.2) for BES channels, as in in Figure 5(b); (b) the same as (a) but spike-filtered; (c) the same as (a) but spike- and then bandpass-filtered in the range ; (d) the detected signal of the BES diagnostic when the shutter is closed, measured during shot #27367 at , using the channels inner BES subarray; (e) the same as (d) but spike-filtered; (f) same as (d) but spike- and then bandpass-filtered in the range . Note the different -axis scales in (a-c) and (d-f).

a.1 Accounting and correcting for spurious skewness

The distribution of fluctuating intensities of the raw BES signal can be skewed for two main reasons not associated with the symmetry breaking of the turbulent field: the presence of high-energy radiation and the slow variation of the intensity signal. We discuss these effects by considering Figure 16(a), where the distribution of (unnormalised) fluctuating intensities from the raw BES signal for the DLM case (see Section 2.4) is plotted.

First, let us consider the smaller (in this case) of the two effects: the low-probability tail on the right-hand side of the distribution in Figure 16(a). This tail is caused by high-energy radiation (neutron, gamma ray, or hard -ray) impinging on the avalanche-photodiode detector (APD) of the BES camera and producing spikes in the time series that are isolated to a single detector channel.

We can observe this effect more clearly by analysing BES data from another MAST shot (#27367), where the shutter over the optical system of the BES diagnostic was closed, and so no beam-emission photons were able to reach the APD. The distribution of unnormalised fluctuating intensities for this “shutter-closed” case is plotted in Figure 16(d). It can be seen to be made up of two components: the first is a Gaussian distribution of photon (electronic) noise, the second is a tail of positive due to the high-energy radiation, causing large positive skewness .

The sum of two randomly distributed independent variables (here, beam-emission photons and high-energy radiation) has a probability distribution that is the convolution of the probability distributions of the two variables (see, e.g., [49]). Therefore, formally, knowing (having measured) the probability distribution of the radiation spikes, it is possible to deconvolve this probability distribution from an experimentally measured probability distribution. However, this is not possible in practice because the actual number of large-amplitude radiation spikes is small (typically at most one or two per channel in a time window). As a result, we resort to a cruder method, but an effective one, of removing the radiation spikes from the time series by identifying large (above a certain threshold) differences between the intensity at one point in time and the next, and then replacing this large-intensity value with the value of the neighbouring point17 [15]. The result of this “spike filtering” can be seen in Figure 16(e), where the large tail of the distribution in Figure 16(d) has been successfully removed by this method.

However, our spike filter is not able to remove all of the radiation spikes, especially ones with smaller amplitudes. These remaining low-amplitude radiation spikes cause high-frequency variations in the time series. Therefore, a low-pass filter can be used to remove them. A low-pass filter is also necessary to remove high-frequency (above ) noise from the signal (see Figure 11 of [15]). In order to account for the slow variation of the intensity signal (which will be discussed shortly), a high-pass filter is also required. Therefore, we use a bandpass filter in the range on the BES-measured time series. The outcome is shown in Figure 16(f). The skewness of the fluctuating-intensity distribution in this figure is nearly zero, demonstrating that this process can successfully account for the skewness caused by high-energy radiation and gives us an estimate of the uncertainty with which we can measure the skewness of the distribution of the turbulent field (in the absence of other effects; see Appendix A.2). Returning to the DLM case, in Figure 16(b), we see that the spike filter successfully removes the low-probability tail.

Now we consider the second spurious source of skewness in the distribution: the slow () time evolution of the intensity signal due to changes in the equilibrium and slow MHD modes (). If, during a given time interval, the (running) mean intensity signal varies unevenly (for example, decreases slightly for the last third of the time interval), then the mean over the entire time interval will be weighted by this variation. The distribution of the fluctuating part of the intensity signal, , will then be skewed. Therefore, a high-pass filter is used to remove all variations of the intensity signal with frequencies below . We demonstrate the effect of using such a filter (in combination with the low-pass filter discussed above) on the DLM case, which has, on inspection of the raw time series, a slow, uneven variation in the intensity signal over a time interval. In Figure 16(c), the bandpass filter reduces the skewness in the distribution significantly, compared to that in Figure 16(b). We have now recovered the distribution that was shown in Figure 5(b).

We finish by noting another possible source of spurious skewness in the distribution of fluctuating intensities, the background emission from the plasma. Intensity distributions for the cases either when there is no neutral beam heating or when only the SW neutral beam is active (the BES images the SS beam), are similar to that of the shutter-closed case in Figure 16(d). Therefore, we do not expect the background plasma emission to produce a skewness larger than that caused by the high-energy radiation, and that we have been able to account and correct for.

Figure 17: Effect of PSFs on the distribution of amplitudes from the GKa5, GKa1, and Marginal simulations. The intensity fields (d-l) are calculated from the density fields (a-c) by evaluating (18). The dashed line in each plot is the normal distribution.

a.2 PSF effects

We now consider PSF effects on the distribution of the fluctuating field. In the first row of Figure 17, panels (a-c), we reproduce the three distributions of the numerically simulated fluctuating density field already plotted in Figure 12 (these are the distributions for the GKa5, GKa1, and Marginal runs). In each of the following rows, we plot the distributions of the fluctuating intensity field calculated from the density field used in panels (a-c) and using the PSFs of the BES diagnostic calculated for each experimental case that we considered in Section 2.4: the second row, panels (d-f), is for the BLM case, the third row, panels (g-i), is for the DLM case, and the fourth row, panels (j-l), is for the IFS case. The relationship between the density field, the intensity field and the PSFs is given by


where is the PSF of channel , calculated according to [46], and