Comparison of BES measurements of ion-scale turbulence with direct gyrokinetic simulations of MAST L-mode plasmas
Observations of ion-scale () density turbulence of relative amplitude are available on the Mega Amp Spherical Tokamak (MAST) using a 2D (8 radial 4 poloidal channel) imaging Beam Emission Spectroscopy (BES) diagnostic. Spatial and temporal characteristics of this turbulence, i.e., amplitudes, correlation times, radial and perpendicular correlation lengths and apparent phase velocities of the density contours, are determined by means of correlation analysis. For a low-density, L-mode discharge with strong equilibrium flow shear exhibiting an internal transport barrier (ITB) in the ion channel, the observed turbulence characteristics are compared with synthetic density turbulence data generated from global, non-linear, gyro-kinetic simulations using the particle-in-cell (PIC) code NEMORB. This validation exercise highlights the need to include increasingly sophisticated physics, e.g., kinetic treatment of trapped electrons, equilibrium flow shear and collisions, to reproduce most of the characteristics of the observed turbulence. Even so, significant discrepancies remain: an underprediction by the simulations of the turbulence amplituide and heat flux at plasma periphery and the finding that the correlation times of the numerically simulated turbulence are typically two orders of magnitude longer than those measured in MAST. Comparison of these correlation times with various linear timescales suggests that, while the measured turbulence is strong and may be ‘critically balanced’, the simulated turbulence is weak.
A ‘grand challenge’ of fusion research is to develop reliable, first-principles predictive capability of plasma confinement, e.g., of the ITER device. An essential part of this process is to perform quantitative comparisons of simulation results with experimental observations and thereby assess their validity . Prediction of heat and particle fluxes require non-linear simulations of the saturated state of plasma turbulence, which can be validated at one level by comparing with the results of transport analysis and, at a deeper level, with measured turbulence characteristics. An excellent overview of the methodology of such studies can be found in Ref. .
The micro-instabilities responsible for anomalous transport exist over a wide range of spatial scales: from electro-static ion-temperature-gradient (ITG), trapped-electron (TEM) and parallel-velocity gradient (PVG) modes (or electro-magnetic micro-tearing modes at high-) at scales larger than the ion Larmor radius (where is the wavenumber perpendicular to the magnetic field, within a flux surface and is the ion Larmor radius), down to electron-temperature-gradient (ETG) modes with . A wide range of diagnostic techniques is required to detect the resulting fluctuations in plasma parameters. Such verification programmes have been the focus of much research activity on conventional tokamaks, e.g., on DIII-D, which has a comprehensive set of turbulence diagnostics [4, 2, 5] and on Tore Supra [6, 7]. Studies on the NSTX spherical torus have focused on electron-scale turbulence, comparing data from a micro-wave scattering system with global, non-linear gyro-kinetic simulations using GYRO  (earlier studies focussed on linear gyro-kinetic calculations [9, 10]). The availability of ion-scale density fluctuation measurements from BES  will facilitate future multi-scale validation studies on the NSTX device.
On DIII-D, a unique array of multi-scale, multi-field turbulence diagnostics has facilitated detailed validation studies comparing gyro-kinetic turbulence simulations with experimental data from L- and H-mode plasmas [13, 14, 15]. A variety of studies have been undertaken, e.g., of the dependence of the turbulence characteristics on  and elongation  and of the cross-phase between and fluctuations , thus probing the underlying physics, especially in cases where disagreement was found. Meaningful comparison of simulation results with experiment is only possible by implementing synthetic diagnostics that mimic the instrumental characteristics of the various fluctuation measurements, as was done, e.g., on DIII-D for BES and Correlation Electron Cyclotron Emission (CECE) diagnostics for measurements of and , respectively [18, 17, 14]. In the L-mode studies, although good agreement was found between simulated and measured heat flux and fluctuation characteristics in the mid-core region (, where denotes the normalised radius), simulations generally underpredict the heat fluxes and amplitudes by almost an order of magnitude in the peripheral region () [17, 14, 4, 15, 2]. This discrepancy has been attributed to increasing importance of edge-core coupling  or to an inward propagation of turbulence (or ‘avalanche’) from the plasma edge . More recently, there have been concerted efforts to resolve this discrepancy, which now appears not to be ubiquitous.
No such validation excercise using global, non-linear simulations of ion-scale turbulence has as yet been performed for a spherical-tokamak plasma, whose particular characteristics make this especially interesting. Heating of the low-aspect-ratio plasma with tangentially directed neutral-beam-injection (NBI) heating results in strong equilibrium rotation in the core (toroidal Mach number , where is the toroidal rotation rate and is the ion thermal velocity) and hence strong shear , where and the safety factor, which can be strong enough to stabilise ion-scale turbulence. The low aspect ratio also results in a large trapped particle fraction (), enhancing the drive for Trapped-Electron Mode (TEM) micro-instabilities, although, as will be seen, collisions can reduce the trapped particle drive substantially in MAST plasmas. Finally, the low toroidal magnetic field results in larger values of the ion Larmor radius normalised to the plasma radius, . The assumption of constant gradient scale lengths in local gyro-kinetic simulations is questionable when gradients vary appreciably over the simulation domain, which typically has to be several in radial extent to ensure convergence. Therefore, as discussed in Ref. , for current-generation ST plasmas in which , global codes may be required to perform meaningful non-linear simulations of the ion-scale turbulence.
Here we present the first validation exercise for a MAST plasma comparing the characteristics of ion-scale turbulence measured using a 2D BES turbulence imaging system  with synthetic data produced from non-linear simulations performed using the global particle-in-cell (PIC) code NEMORB . The comparison is performed for a low-density, L-mode plasma exhibiting an internal transport barrier (ITB), similar to those used for earlier studies of ITB formation and dynamics [24, 25], for which linear gyro-kinetic stability calculations have been performed using the local code GS2 . This equilibrium configuration was also used as the basis for the earlier global, non-linear simulations presented in Ref. . An L-mode plasma is ideal for these studies because, in the peripheral region, the equilibrium flow shear is too weak to stabilise the ITG turbulence fully.
The MAST device is a medium-sized, low-aspect-ratio tokamak (aspect ratio , plasma current , toroidal field at 0.7 m), which is equipped with tangentially directed NBI heating (injected power at injection energy). Studies of transport in MAST are facilitated by the availability of high-resolution kinetic profile data from a suite of advanced diagnostics (see §3) and an integrated analysis chain () to prepare this data for transport analysis using TRANSP . The BES turbulence imaging system on MAST has sufficient signal-to-noise ratio (SNR) to detect density fluctuations with relative magnitude if appropriate correlation analysis is used. Calculations of the three-dimensional spatial response and sensitivity of the BES system, accounting for relevant physical effects were performed in Ref.  as described in §3.1. This is quantified in terms of two-dimensional point-spread functions (PSFs), which are used here to generate synthetic BES data from the gyro-kinetic simulation data. Previous results of an analysis of the BES data  to elicit the structure and dynamics of the ion-scale turbulence  and the dependence of the ion temperature gradient on the magnetic configuration, rotational shear and heat flux  are pertinent to our discussion of the results presented here.
The remainder of the paper in structured as follows. The results of gyro-kinetic simulations of MAST L-mode equilibria are presented in §2, with a summary of results of previous linear calculations in §2.1 and the non-linear simulations in §2.2. Our experimental observations of the ion-scale turbulence are explained in §3, with a brief description of the capabilities of the BES system in §3.1 and an explanation of the correlation analysis used to determine the turbulence characteristics in §3.2. The method used for the generation of synthetic data is then explained in §3.3. The observed and simulated turbulence characteristics are then compared in §4. Finally, our conclusions are presented in §5.
2 Gyro-kinetic simulations of an L-mode discharge
The validation study presented here is for an L-mode discharge formed using a scenario with early NBI heating during the current ramp, which slows current penetration, giving rise to strong toroidal rotation and negative magnetic shear () in the plasma core. Previous linear-stability calculations for such a discharge  revealed the peripheral region, where the flow shear is much weaker, to be unstable to ITG modes. Therefore, such discharges hosting ITG turbulence are ideal candidates on which to validate non-linear, ion-scale turbulence simulations.
Kinetic profiles of the equilibrium used for this study are shown in Fig. 1. They are taken from MAST L-mode discharge #27268 at 0.25 s, the parameters of which are: plasma current , toroidal field at , co-injected NBI heating of , line-average density , with , and , corresponding to Mach number . The foot of the ITB in the ion thermal and momentum channels is located at the normalised radius (where is the normalised poloidal flux) in the region where . Note that the temperature and rotation profiles of the bulk ions are calculated from measurements on the impurities using the NCLASS  neo-classical package within the TRANSP code.
Linear-stability calculations were performed using both GS2  and NEMORB  for a similar equilibrium from an earlier MAST discharge (#22087 at 0.25 s), albeit with somewhat higher plasma current of , which was used for a study of ITB formation and evolution . The kinetic profiles from this equilibrium are very similar to those of the equilibrium studied here, so the results of these linear calculations are relevant and briefly summarised in §2.1. This equilibrium was also used as the basis for first global, non-linear simulations  using NEMORB and some of the relevant results are also summarised in §2.2.
2.1 Linear stability calculations of #22807 at 0.25 s.
As presented in Ref. , local linear-stability calculations with kinetic electrons were performed using GS2 for the equilibrium from discharge #22807 at the time of peak ITB strength at at three locations: inside the ITB, just outside and in the plasma periphery. At the innermost of these surfaces (at = 0.3, where is the normalised toroidal flux, corresponding to ) all modes at electron and ion scales were found to be stable in both electrostatic and electromagnetic calculations. At mid-radius (, corresponding to ), in the absence of flow shear, the relatively low collisionality results in appreciable trapped-electron drive and TEM modes are unstable in the intermediate wave-number range () – these are not completely stabilised by the flow shear. (In contrast, in calculations with adiabatic electrons, ITG modes that are unstable without flow shear are completely stabilised when the flow shear is included). At the outermost surface (, corresponding to ), the TEM modes are stable because the plasma is more collisional, while strongly growing ITG modes are not stabilised by the flow shear, which is weak at this radius.
To assess the importance of non-local effects, the results of linear stability calculations with NEMORB were compared with those from local, linear GS2 calculations (with and without collisions) . Linear, electrostatic calculations were made with NEMORB (for the equilibrium from #22807 at 0.25 s) with adiabatic and kinetic electrons, both with and without flow shear , but without collisions. As in the GS2 calculations, the core was found to be stable to ITG modes, while the most unstable region was found to be . By varying the normalised ion Larmor radius in the NEMORB simulations, it was found that global effects significantly reduce the growth rates, starting at just below the experimental value of 0.018.
2.2 Global, non-linear simulations with NEMORB
Previous simulations of #22807 at 0.25 s.
The results of global, electrostatic, non-linear simulations with NEMORB of the equilibrium from discharge #22807 at were presented in Ref.  and the reader is referred to this article for details. In the simulations with adiabatic electrons but without sheared flow, the turbulence was found to spread from the linearly unstable region at into the stable region down to during the non-linear phase. The ion heat flux was found to be very low , which is below the neo-classical level. Increasing the normalised, inverse ion-temperature-gradient scale length by 30% increased the heat flux considerably, indicating that the gradient is close to the non-linear threshold for ITG turbulence. In simulations with kinetic electrons but without sheared flow, the heat flux increased to , which is far above the experimental level, . Again, the turbulence spread into the linearly stable core region. Including equilibrium flow suppressed the turbulence in the region with strong shear () and modestly reduced the peak heat flux. Results of runs incorporating electron collisions, which had been shown to reduce the linear drive due to trapped electrons , were not reported in Ref.  (simulations with collisions have however been peformed for the more recent equilibrium discussed in §2.2.2).
Non-linear simulations of #27268 at 0.25 s.
In order to be able to compare our BES measurements with the non-linear simulations, further simulations have been carried out for an equilibrium from which the measurements are available, i.e., that corresponding to the kinetic profiles shown in Fig. 1. Electrostatic simulations are available for five cases, distinguished by whether they were run with adiabatic (AE) or kinetic electrons (KE), with or without ion-electron and electron-electron collisions and with or without equilibrium toroidal flow. The various combinations for the five runs (I-V) are summarised in Table 1, where the start time and duration of the non-linear phase used to generate the synthetic data are also stated. As discussed in §3.2 below, these periods are considered the minimum required to extract reasonable estimates of the turbulence characteristics. The evolution of the maximum heat flux over the radial profile during each of the non-linear runs is shown in Fig. 2.
|Case||KE||Flow||Coll.||Start ||Duration ||Symbol|
Details of numerical set-up.
In NEMORB, the heat source is adjusted until the kinetic profiles match those prescribed and the resulting heat flux is a prediction which can be compared with the experimental value, although, as is discussed in §4.1, a precise match cannot be expected. In order to maintain the signal-to-noise ratio for the duration of the simulation, markers were used on an grid, for the simulations with collisions the number of markers was doubled. For the simulations with kinetic electrons, the trapped-electron response is treated kinetically while that of the passing electrons is treated adiabatically. The reader is referred to Ref.  for further numerical details.
In order to ensure numerical stability, a boundary condition is adopted in NEMORB which forces both the density and potential perturbations at to be zero, . Markers near the plasma boundary carry some density which is zeroed just before the Poisson equation is solved. Moreover, if a marker or one of its associated gyro-points lies outside the plasma, it is not taken into account. Therefore, the quasi-neutrality condition is violated near the boundary and charge accumulation can occur, leading to a numerical sheath region with spurious electric fields. In order to circumvent this problem, the Poisson equation is solved with an artificial ‘shielding’ density in addition to the prescribed density near the boundary. The shielding density, as shown in Fig. 1(b), peaks at the plasma boundary, decaying exponentially over a scale length . Sensitivity studies have been carried out to determine the effect of changing the scale length of the shielding density, whereby it was found that increasing this from to 0.03 has little effect on the amplitude profile, which in all simulated cases peaks deeper inside the plasma (at or inside ).
As discussed in §3.3, generation of the synthetic BES data requires output of the simulated density fluctuations as a function of time over a 2D (radial
3 Experimental measurements
The ion-scale density fluctuations in MAST plasmas are measured using a 2D imaging BES diagnostic , which is briefly described in §3.1. The characteristics of the turbulence (amplitude, correlation times and lengths and apparent poloidal phase velocity of the density contours) are determined using the correlation analysis described in §3.2, applied to short () time series corresponding to the equilibrium used for the gyro-kinetic simulations. Care was taken to choose a quiescent time period free of the fast-ion-driven MHD events, with frequencies in the range , which are frequent in NBI heated discharges on MAST. Other data required for our analysis are obtained from the high-resolution profile diagnostics available on MAST: the magnetic pitch angle (, where and are the poloidal and toroidal components of the magnetic field) is measured by a 32-channel Motional Stark Effect (MSE) diagnostic ; ion temperature and toroidal flow velocity are obtained from multi-channel, Charge-Exchange-Recombination Spectroscopy (CXRS) measurements on impurity ions with spatial resolution of ; and electron temperatures and densities from a 132-channel NdYAG Thomson Scattering system with comparable spatial resolution . The components of the magnetic field vector are obtained from MSE-constrained EFIT equilibrium reconstructions .
3.1 BES turbulence measurements
The ion-scale density fluctuations are measured using a 2D imaging BES diagnostic . This system has a 2D Avalanche Photo-diode Detector (APD) in the form of an 8 radial 4 poloidal channel array, which is directly imaged (rather than using optical fibres) onto the heating neutral beam with a resulting nominal spatial resolution of radially and poloidally. The Doppler-shifted emission from the beam is selected using a band-pass interference filter. The incident light on the APD sensors with a photon flux of is detected with a SNR at 2 MHz digitisation rate simultaneously for all channels. If appropriate correlation analysis is used , the system is able to detect density fluctuations with relative amplitude at wave numbers and at frequencies up to the Nyquist frequency of 1 MHz.
The actual spatial resolution is degraded below that determined by the imaging properties of the optical system, both due to geometrical effects caused by the finite depth of the line of sight through the beam (), together with field-line curvature and any mismatch between the -field direction and the line of sight, which is optimal when viewing at mid-radius (). The finite lifetime of the excited deuterium atoms () also causes the emission due to instantaneous excitation to be spatially de-localised by in the direction of beam propagation. The spatial response of the BES system is determined using a numerical simulation code , which takes account of these effects in terms of 2D point-spread functions, as discussed in §3.3. This spatial response depends on the particular magnetic equilibrium, beam parameters, viewing location and plasma profiles and so has to be calculated explicitly for each measurement.
The emissivity of the beam is proportional to the beam density , the plasma density at the observed location and to the relative population of the excited state, which is determined by a collisional-radiative balance between excitation and de-excitation of the beam atoms by collisions with the plasma ions and electrons and by radiative decay. Neglecting any low-frequency (a few kHz) fluctuations in the beam density, the relative density fluctuation level can be determined from the intensity fluctuations using the relation , where () and () are the fluctuating and mean components of the plasma density (intensity of emission), respectively . The differential excitation rate is obtained from the data in Ref. , where it is calculated using an appropriate collisional-radiative model. This parameter is a weak function of density (decreasing from 0.93 to 0.27 over the density ranging from to ) and a very weak function of temperature.
The ability of the BES diagnostic to detect density fluctuations due to low-frequency, ion-scale turbulence - the phenomenon of interest in this study - can be compromised by density fluctuations due to fast-ion-induced MHD activity superimposed on the data. Examples of cross-power and cross-phase spectra (between two poloidally adjacent channels at the same nominal radius corresponding to ) of the BES fluctuation data are shown in Fig. 3 for two periods corresponding to the MHD-quiescent phase used for the correlation analysis and a period just after during which there is fast-ion-induced MHD activity. During the MHD-free phase, the turbulent component of the signal decreases to the level of the broad-band noise () at frequencies . The linear increase in the cross-phase with frequency up to observed during the first period is expected for broadband turbulence propagating poloidally at constant phase velocity. During the second period, the cross-power increases by an order of magnitude and is dominated by the MHD component, which exhibits a wide range of frequencies due to the sweeping nature of the fast-ion driven ‘fish-bone’ instabilities . The spatially coherent nature of the MHD-induced fluctuations is manifest in the zero cross-phase over all frequencies below . In order to avoid the complications imposed by the MHD, we have applied the correlation analysis described in §3.2 only to short time periods relatively free of this activity. (In previous work , we applied data selection cuts which reject data from periods with strong, coherent MHD activity.)
Another general issue with core fluctuation diagnostics based on BES, which was first raised with reference to the measurements on TFTR  is that strong density fluctuations at the plasma periphery can be imprinted on the beam density, hence superimposing a spatially coherent component onto the observed fluctuations in the plasma core, which is in anti-phase with those at the edge. Methods exist to account for this effect by performing an inversion of an integral transformation, involving the so-called ’beam-transfer function’, over the beam path from the edge to the observed location in the plasma, e.g. as described in Ref. . Such methods can, however, only be applied if the fluctuation measurements are available simultaneously over the beam path from the edge to the deepest observed location. As the MAST BES system has a fixed detector geometry, it is impossible to observe the edge and core plasma simultaneously, so such techniques cannot be applied to our data. However, the correlation analysis described in §3.2 accounts for spatially constant contributions to the measured correlation functions. Let us nevertheless give a simple estimate that shows that the beam-imprinting effect on our measurements is small.
The detection limit for core fluctuations in terms of imposed by this beam imprinting effect due to a given level of edge fluctuations can be estimated as follows. The emissivity from the beam is proportional to the beam density , hence fractional changes in the beam density will produce a similar fluctuation in the observed intensity and a factor larger fractional change in the inferred density fluctuation level. An approximate upper limit to the apparent core fluctuation level due to an edge density fluctuation in an edge-localized ‘shell’ of line-integral density , where is the shell’s thickness, is , where is the beam stopping rate coefficient and is the beam-atom velocity (for a beam at ). Estimating the average fluctuation level as over a peripheral shell with (equivalent to the region ) yields , which is below the observed fluctuation level in the plasma core inferred from our measurements (see Fig. 6)
3.2 Correlation analysis
The statistical characteristics of the fluctuations are determined using correlation analysis techniques very similar to those used in Ref. , to which the reader is referred for further details. The data is first band-pass filtered over the frequency interval to reject high-frequency noise and the beam fluctuations at a few kHz, which otherwise disturb the correlation analysis . A matrix of spatio-temporal correlation functions is then calculated according to:
where , and denote the radial and vertical (poloidal) positions and time, respectively, and are the radial and poloidal channel separations, is the time lag and denotes a time average over a period of for the experimental data and over the periods given in Table 1 for the simulated data. Note that we use the notation () for the radial and perpendicular directions relative to the flux surfaces, which for the case of the up/down symmetric, double-null diverted (DND) equilibria considered here are aligned with the () directions at the horizontal mid-plane.
The auto-covariances at , contain not only a component due to the true density fluctuations but also photon and electronic noise. The density fluctuation level at each radial location can be obtained from the (noise-subtracted) auto-covariance functions according to: , where is the auto-covariance of a signal from a calibration source containing only noise, determined at all 32 locations, then averaged (as denoted by ) over the four poloidally separated channels at the same radial location.
The poloidal correlation length is estimated using data from the four poloidally spaced channels at each radial location by fitting the averaged cross-correlations (over channels with the same at each radial location ), to the function , where and are the two fit parameters. The parameter effectively subtracts any spatially constant component, usually due to large-scale, global MHD modes. The field-perpendicular correlation length on a given flux surface is then , where is the pitch angle of the magnetic field. An example of a fitted poloidal correlation function is shown in Fig. 4(a) for channels at a nominal location of for which the polodially constant offset term was . The relatively large value of this offset is likely to be due to residual MHD activity, which is never entirely absent during the beam-heated phase. Note that fitting this offset serves to subtract both any component due to residual, global MHD activity or from edge-induced beam density fluctuations as discussed in §3.1 above. This technique is essentially equivalent to that discussed in Ref. , where the contribution due to edge-induced beam density fluctuations was accounted for by subtracting a long-range, spatially constant component from the local cross-correlation functions.
In contrast to the poloidal correlation functions, the radial correlation functions have a monotonically decaying rather than a wave-like structure. Radial correlation lengths are therefore obtained by fitting the averaged radial cross-correlations (over the four poloidal channels at each radial location ), to the function , where the fit parameters are and the constant , again serves to subtract any spatially constant component. An example of a fitted radial correlation function is shown in Fig. 4(b). The radially constant offset term in this case was , which is quite close to obtained above from fitting the poloidal correlation function.
The correlation time is estimated from the temporal cross-correlation functions between the four poloidally separated channels at each radial location, by finding for a particular the time delay at which the correlation function has its maximum. The peak values of the correlation functions are then fitted with the function to determine , where it is implicitly assumed that any influence of the finite parallel correlation length on can be ignored.
The apparent phase velocity of the fluctuations is determined at each radial location using a cross-correlation time delay (CCTD) technique , namely it is calculated by making a linear fit to . It is found that the observed velocity is dominated by the apparent poloidal motion of field-aligned, elongated () eddies through a poloidal plane due to the dominant toroidal rotation of the plasma: as discussed in Ref. , it can be shown that , where the toroidal velocity . This is because any poloidal flows are strongly damped , leaving of the order of the diamagnetic velocity .
This type of analysis, characterising the turbulence in terms of a few scalar parameters ( and ) determined from the ensemble-averaged correlation functions, effectively yields these parameters averaged over all density fluctuations at spatial scales larger than the instrumental resolution (), weighted by the RMS value of the density fluctuations . These parameters are therefore representative of the fluctuations at the dominant ‘outer’, energy-containing scale of the turbulence . In order to calculate meaningful ensemble averages, the time averaging must be done over many correlation times and wave periods, i.e., over a sample of duration , where is the ion-diamagnetic frequency and is the scale length of the ion temperature gradient. Taking the fiducial value , sample periods of order 1 ms are required.
3.3 Synthetic BES data
For the comparisons between simulations and measurements to be meaningful, ‘synthetic’ data, analogous to the real BES data, is generated from the simulated density fluctuations and analysed using the same correlation techniques as those used for the experimental measurements. The synthetic data is generated using a synthetic ‘diagnostic’, which mimics the physical characteristics of the measurement technique. The effect of applying this processing to the simulated data is firstly to spatially average over fluctuations at smaller scales than that of the instrumental resolution and secondly to impose a lower limit to the fluctuation amplitude below which the signal is dominated by noise.
In order to generate the synthetic data, the characteristics of the BES system have to be quantified in terms of its spatial resolution, sensitivity and noise properties. These are determined using a numerical simulation, described in more detail in Ref. , which accounts for the beam absorption and emission, magnetic field and viewing geometry. Assuming that the turbulent eddies are highly elongated along the field, , the spatial response of each channel can be quantified in terms of 2D point-spread functions (PSFs), , which are effectively cross-sections for the excitation of the beam emission per unit area of the field of view at the focal plane at the beam.
The PSFs depend on the beam properties, magnetic equilibrium and plasma profiles, as well as on the diagnostic parameters (e.g., on the APD bias voltage and radial viewing location) and so have to be calculated explicitly for each simulated equilibrium and diagnostic setting. The validity of this procedure of utilising 2D PSFs to represent the spatial response relies on the assumption that the turbulence is to a good approximation two-dimensional, specifically that its parallel correlation length , where is the line-of-sight path length through the beam. If the parallel correlation length , where is the connection length at the low-field side of the plasma [30, 47], then this approximation is well satisfied.
Generation of synthetic BES data from the non-linear simulations requires a 2D map of the relative density fluctuations over a poloidal cross section (i.e. at fixed toroidal angle) as a function of time on a regular () grid covering the region of the plasma viewed by the BES diagnostic. With three viewing locations used to measure a full radial profile, the BES measurements cover a radial range of , while the vertical extent of the measurements is . Relative density fluctuations from a NEMORB simulation over such a radial stripe are shown in Fig. 5 with the PSFs of the BES system superimposed. Because the spatial extent of each channel is of order of a few cm, a resolution of is sufficient for the simulated data. The sampling rate of the BES diagnostic is 2 MHz, hence the sample period for the simulated data is also chosen to be .
As well as the relative density fluctuation data, some other data is required for the synthetic data generation. Calculation of the rate of photons incident on each sensor channel requires the mean density (which does not need to be a function of time as the duration of the simulation is much shorter than the timescale of the profile evolution) and the simulated relative density fluctuation data. The photon rate can then be calculated from the following expression:
where the differential excitation rate was discussed in §3.1 and the integral is performed over the region covered by the simulated data.
The statistical photon noise within the digitisation period is simulated as follows. The number of detected photons in a sample period , where is the bandwidth of the pre-amplifiers, is , where is the quantum efficiency of the detectors and is a factor to take account of the transmission of the optics. Pseudo-photon noise is then added according to , where is a random number from a normal distribution of unit standard deviation. The coefficient is the excess noise factor due to the amplification process in the APD, which is given by , where is the excess noise exponent due to the amplification processes in the APD, and is the gain of the APD ( at the bias voltage used for the measurements of 310 V). The output voltage is then calculated from , where is the electron charge and is the trans-impedance of the amplifiers. The second term is the electronic noise of the amplifier, also simulated by adding a normally distributed random noise voltage with standard deviation . Typically, the output voltage is , hence the SNR due to the amplifier noise alone is and the overall SNR including the photon noise is .
The procedure followed to generate the turbulence characteristics for the simulated turbulence data can be summarised as follows: PSFs for each of the three BES view radii are generated for the specific equilibrium and beam parameters used in the experiment; data is written out from the simulation over the required domain; the synthetic BES data are generated for each available PSF; finally, the synthetic data is analysed using the same correlation analysis as for the real BES data, which was described in §3.2.
4 Comparison of non-linear simulations with measurements
Results of applying the above procedure to the five non-linear NEMORB simulations of the MAST L-mode discharge #27268 discussed in §2.2.2 are presented below and the simulated turbulence characteristics are then compared with experimental measurements. Profiles of the experimental and predicted fluctuation amplitudes and ion heat fluxes are shown in Fig. 6 as a function of normalised radius . Because the radial extent of the BES array covers only one third of the outboard plasma radius of , the experimental data in this figure are taken from three similar discharges, with nominally the same equilibrium and kinetic profiles, but with the BES viewing position set to three different radial locations (#27272, 1.05 m, #27268, 1.2 m and #27274, 1.35 m). The correlation analysis is applied to a data sample at the time , corresponding to the equilibrium used for the simulations (note that this is before the time of the beam cut-off at 0.3 s).
4.1 Turbulence amplitudes and ion heat fluxes
As shown in Fig. 6 (a), the experimental fluctuation level in the core is , which is significantly above that due to background light emission of and increases strongly to reach towards the periphery. The background fluctuation level is estimated by normalising the fluctuation amplitude just after the beam cutoff to the signal level just before the beam cut. Hence, should be representative of the background fluctuation level during the beam phase provided this does not change significantly between the time of observation and the beam cut-off. Spectral observations show that the background emission is predominantly so-called FIDA (fast-ion emission) from re-neutralised beam ions, which increases towards the plasma periphery where the neutral density is highest. Note that, for , the observed fluctuation level in the core is larger than the detection limit estimated in §3.1 as due to edge fluctuation induced beam density fluctuations.
The experimental heat flux, determined from transport analysis  performed using the TRANSP code , is shown in Fig. 6 (c). The uncertainties in the absolute levels of are quite large due to the difficulty in quantifying the level of anomalous fast-ion redistribution and losses arising from fast-ion driven MHD activity, which is particularly strong during the phase with both on-axis directed NBI heating beams for which this comparison is made. Outside , the heat flux is approximately constant at and the value normalised to the gyro-Bohm level , where , increases by over three orders of magnitude from in the core to in the plasma periphery. In the outer region of the plasma, is up to a factor larger than the neo-classical heat flux determined from NCLASS , but it is below the neo-classical level in the core, where the power balance is particularly sensitive to assumptions made about the anomalous fast-ion diffusion.
Determining the turbulent heat flux in principle requires knowledge of the density, temperature and potential fluctuations and cross-phases between them . Measurements of potential fluctuations in the core plasma are challenging and rarely available  and those of ion-temperature fluctuations even more so . In the absence of such measurements on MAST, the gyro-Bohm-normalised turbulent ion heat flux can be estimated from the simplified relation: , where the Boltzmann relation between and is assumed and any possible reductions due to the cross-phase are ignored. As can be seen in Fig. 6 (c), this simplified estimate agrees well with that obtained from power balance at the periphery , but underestimates by a factor over the rest of the profile. Interestingly, measurements on TFTR  found that in plasmas with dominant ITG turbulence, over the outer half of the plasma radius, so an appreciable fraction of the ion heat flux might arise from enhanced ion temperature fluctuations in regions with dominant ITG turbulence.
The normalised shearing rate, , is shown in Fig. 7, together with maximum growth rates from linear GS2 calculations (in this case from the similar equilibrium of discharge #22807 at 0.25 s) for two cases (I and II) without sheared flow or collisions. While for case I (adiabatic electrons), over the whole profile, for case II (kinetic electrons), only in the region with strong flow shear inside , which is consistent with the turbulence being suppressed by the flow shear there.
As shown in Fig. 6 (d), the simulation with adiabatic electrons but without sheared flow (I) results in an ion heat flux up to an order of magnitude below and a correspondingly small fluctuation amplitude. Kinetic treatment of the trapped electrons only, again for a case without flow (II), results in values of up to an order of magnitude larger than at mid-radius, while in the core and peripheral regions is underpredicted. The fluctuation level exhibits similar behaviour. The underprediction in the core is consistent with the previous linear studies , which showed the region where to be linearly stable. The effect of including equilibrium flow shear (III) in the simulations with kinetic electrons is not dramatic, with a slight decrease in compared to case II in the core where the flow shear is strongest, but an increase in the peripheral region, where approaches the experimental level – the corresponding changes in are rather slight.
Introducing collisions (IV and V) reduces the drive from trapped electrons, hence decreasing the level of turbulence and the heat flux. The effect of this is to reduce the radial extent over which there is significant turbulence, in the case IV (without flow) to a limited region around , while introducing sheared flow (V) broadens and shifts the unstable region inwards (compared to case IV) to . In case V, with the most complete physics, the peak value of exceeds somewhat in the mid-radius region, while the fluctuation level is close to that measured. In the peripheral region, both of these cases underpredict and (this is not surprising because in NEMORB the fluctuation level is forced to be zero at the plasma boundary).
It is perhaps surprising that introducing sheared flow into case V actually increases the predicted heat flux by almost an order of magnitude at mid-radius, when sheared flow is normally expected to suppress the turbulence. Currently, we do not have an explanation for this result, which would necessiate at least performing further global, linear calculations for the specific equilibrium.
The fact that, particularly in the cases with collisions, the simulated turbulence has significant amplitude in relatively restricted radial regions indicates that the gradient is marginally close to the non-linear critical gradient to excite turbulence. This was borne out by the non-linear studies with NEMORB  carried out for the similar equilibrium of discharge #22807 (see §2.2), which showed that relatively small increases in above the experimental value could produce large increases in the turbulent heat flux. Although such calculations have not been carried out for the equilibrium used for these studies, comparison of similar non-linear simulations for both equilibria show to be even closer to marginality in the equilibrium from #27268 due to the relatively smaller predicted heat fluxes.
Using our database of equilibrium and turbulence data, primarily from MAST L-mode discharges, it is found [29, 31] that exhibits a dependence on flow shear and magnetic geometry, specifically on the parameters and . This dependence shows remarkable similarity to a numerically predicted ‘manifold’ of (where ), required for the excitation of marginally unstable ITG- or PVG-driven turbulence . This observation is further evidence that the is generally close to marginal stability in MAST L-mode plasmas. A consequence of this closeness of to marginality is that the use of the heating operator in NEMORB to match the predicted profile to that prescribed will yield heat fluxes which are very sensitive to experimental uncertainties. Hence, it is not surprising that there is not better agreement between predicted profile of and experiment.
The shortfall in the predicted level of turbulence and heat flux in the peripheral region is a clear deficiency of these non-linear NEMORB simulations, which probably arises from the necessary boundary conditions (see §2.2.3) that suppress any turbulence at the plasma edge. As mentioned in the introduction, such a shortfall has also been found in local, non-linear simulations of DIII-D L-mode discharges using the global code GYRO [13, 14, 2]. In these studies, although good agreement was found in the plasma core between simulated and experimental ion and electron heat fluxes and turbulence characteristics, a systematic shortfall in the heat flux and fluctuation levels was observed at . More recently, concerted efforts have been made by several groups to either understand or resolve this discrepancy [52, 53] and currently this apparent shortfall appears not to be ubiquitous.
4.2 Perpendicular and radial correlation lengths
As shown in Fig. 8, the measured perpendicular correlation length is approximately constant, whereas the radial correlation length decreases somewhat with radius. The observed mean anisotropy is consistent with the findings of Ref. , where it was calculated over a much larger database and was . The corresponding values of (where ) are and , generally decreasing with increasing radius.
In all of the simulations, the perpendicular correlation length is comparable to that observed over most of the profile except in the periphery (), where the simulated turbulence may be impacted by boundary conditions. In the simulations, the radial correlation lengths agree within a factor with the observed values, although the degree of anisotropy in the simulations is somewhat larger. Note in this context that a related NEMORB study  of the dependence of the linear stability of ITG modes on sheared flows showed the anisptropy of the unstable modes generally increases with the shearing rate.
4.3 Correlation times
The measured correlation times, shown in Fig. 9, are rather short and in the range . One of the most striking of our observations is that in all of the simulations, the correlation times are almost two orders of magnitude longer than those measured. Note that the values of for the simulation cases in Fig. 9 had to be calculated using synthetic data generated without adding noise. This is because, when the turbulence amplitude is weak, the determination of is particularly subject to the effect of noise on the correlation functions and not many valid data points remained for cases IV and V. Longer-duration simulations would help overcome this problem but prohibitively long computational times would be required. Note, however, that the derived turbulence characteristics were largely insensitive to the addition of noise.
It was shown in Ref.  that the measured are comparable with various linear time scales. These time scales are: the drift-wave time , the parallel streaming time (thermal ion transit time over the parallel connection length ) and the magnetic drift time (perpendicular magnetic drift time over a radial correlation length ). This observation is consistent with the turbulence being ‘critically balanced’  and implies that the turbulence is anisotropic in the poloidal plane with , where is taken as the shorter of and . It was also found that the observed correlation times are always shorter than or comparable to the shearing time , i.e., , which implies that equilibrium shear does not always determine in the experiment.
In Fig. 10, the correlation times for the L-mode experimental data considered here, together with those from the five NEMORB simulations, are compared with these time scales. For this experimental data, although , the observed are usually somewhat shorter: the mean logarithms of ratios of these timescales for the experimental data are: , , and . This is consistent with strong, critically balanced turbulence. In the case of the simulation results, the derived are substantially longer than any of the linear time scales. This indicates that in these simulations, the turbulence is weak (whereby ), in contrast to the turbulence measured in the experiment, which appears to be strong ().
The de-correlation of the turbulence by the fluctuating potential from the drift waves is characterised by the non-linear time , which can be estimated from the turbulence characteristics by assuming and are related through the Boltzmann response. Under this assumption, does not include any contribution from the trapped-electron response or toroidally and poloidally symmetric zonal flows, so this timescale is denoted and given by . In Ref. , it was found that the measured correlation times were generally much shorter than and it was conjectured that this could be consistent with a strong zonal component to the turbulence contributing dominantly to its de-correlation (the ratio was also found to increase with collisionality, suggesting that the ratio of the zonal to the drift-wave component of the turbulent amplitude increases with decreasing collisionality).
In the experiment considered here, we also find that the non-linear time is always substantially longer than (see Fig. 10(e)), the mean logarithmic ratio being . In contrast, in our simulations, the non-linear time is mostly comparable to (e.g. for case V, ), which indicates that the simulated turbulence is de-correlated by the drift-wave potential fluctuations and a dominant component from zonal flows need not be invoked to explain the correlation times. Note that with the damping rate for ion-scale zonal flows scaling as (perpendicular ion viscosity), non-linear simulations of long duration () would be required to capture correctly the zonal-flow dynamics.
4.4 Poloidal propagation velocities
The interpretation of the poloidal propagation velocity of the density patterns observed with the BES system is discussed in some detail in Ref. . This motion is not that of the poloidal plasma flow but arises largely from the dominant toroidal flow of the field-aligned, elongated ‘eddies’ moving through the radial-poloidal focal plane. In this case, it is evident from simple geometry that this apparent velocity , where is the pitch angle of the local magnetic field line. This expression can also be derived from the continuity equation as in Ref. , where various terms were neglected, including the net, temperature-gradient-driven poloidal flow of the bulk ions, . Taking to be the velocity of the ions measured by the CXRS system , we can also expect differences between this velocity and that of the ions of .
The observed velocity can also be considered as the poloidal projection of the field-aligned density patterns propagating with velocity perpendicular to , i.e. , the parallel velocity component not contributing to any apparent poloidal motion. Neglecting the diamagnetically small net poloidal plasma flow, this can be related to as , which is the same relation as above. The perpendicular velocity , where is the radial electric field component, vanishes in the moving frame of the plasma, so any residual perpendicular flow relative to represents the phase velocity of the density patterns in the plasma frame due to all the effects, as discussed in Ref. .
The observed propagation velocity is shown in Fig. 11(a) together with that calculated from the expression using the velocity from CXRS, where it can be seen that the fluctuations propagate mostly in the positive, ion-diamagnetic direction relative to . The poloidal projection of the velocity, is also shown in Fig. 11(a). This is determined using the profile from TRANSP, which is calculated using the radial force balance . The components of the velocity are calculated in TRANSP from the measured velocity using neo-classical expressions in the NCLASS package .
It can be seen from Fig. 11 (a) that is close to . This can be explained by the approximate cancellation of the diamagnetic contribution to and the difference velocity between the and the ions, which can be approximated as . Considering both of these effects it can be shown that . In Fig. 11, we also plot the sum , where the poloidal projections of the ion and electron diamagnetic velocities are and the pressure-gradient scale lengths are estimated using the TS and CXRS measurements. It can be seen that the observed density patterns propagate in the ion-diamagnetic direction relative to at a velocity .
Poloidal velocity profiles from the simulations with and without equilibrium flow are shown in Fig. 11(a, b), respectively. In case I with adiabatic electrons, the drift velocity is small and in the ion-diamagentic direction. In simulations II and III (kinetic electrons, no collisions), without and with equilibrium flow respectively, the density fluctuations propagate at approximately the electron-diamagnetic drift velocity with respect to , which is consistent with the turbulence being driven by the trapped electrons, and in the opposite relative direction to that observed. In the cases with collisions, the trapped electron drive is reduced and consequently the fluctuations drift more towards the ion-diamagnetic direction than in the cases without collisions, more strongly in case IV (no flow) than in case V (with flow). Some of the cases clearly exhibit radially sheared zonal flows, which are not damped in these simulations due to the absence of collisions. As mentioned above, there are various effects that can affect the apparent velocities, therefore, while the relative differences found between the simulation cases and observations can be considered as an indication of the underlying physics, close quantitative agreement is not expected.
This validation exercise, comparing synthetic BES data generated from global, non-linear gyro-kinetic simulations with observations, has revealed that the simulations with the most complete physics, i.e., including kinetic electrons, sheared equilibrium flow and collisions, exhibit a degree of agreement in terms of the radial and poloidal correlation lengths, the ion heat flux and fluctuation levels (at mid-radius) and the poloidal propagation velocity of the turbulence. There are, however, notable discrepancies that suggest that the simulations have a long way to go before they can be viewed as reliably reproducing reality. Firstly, the measured correlation times are very short and always comparable to the linear timescales (drift-wave drive, parallel streaming and magnetic drift), which is consistent with strong (and probably critically balanced) turbulence, whereas the simulated turbulence exhibits much longer correlation times and appears to be only weakly non-linear. Secondly, the predicted profiles of fluctuation level and turbulent heat flux exhibit some degree of agreement with experiment only in the mid-radius region. This is partly because neo-classical physics is not included in the NEMORB simulations but also because, under conditions where the gradients are close to marginal stability, the prescribed ion-temperature gradient would have to match very precisely in order to yield the experimental heat flux. Finally, the simulations exhibit a marked shortfall in the predicted fluctuation level and ion heat flux in the peripheral region, which is likely both to be a result of the ‘zero-turbulence’ boundary condition required to ensure numerical stability and quite plausibly some missing physics related to non-local interaction of peripheral turbulence with edge/SOL instabilities. Further work is planned to understand these result and to improve the simulations, e.g. to investigate the sensitivity of the predicted heat flux to local changes in the ion-temperature gradient, to resolve the components of the heat flux due to density and ion temperature fluctuations and to examine the possible role of parallel dynamics.
We are grateful to I. Abel, M. Barnes, G. Colyer, S. Cowley, M. Fox, E. Highcock and F. Parra for discussions. This work was funded by the RCUK Energy Programme under grant EP/I501045 and by the European Communities under the Contract of Association between EURATOM and CCFE. To obtain further information on the data and models underlying this paper please contact PublicationsManager@ccfe.ac.uk. The views and opinions expressed herein do not necessarily reflect those of the European Commission. The gyro-kinetic simulations were performed on the HECTOR supercomputer (EPSRC Grant EP/H002081/1) and on HELIOS.
- : Plasma Phys. Control. Fusion
- See for example Ref. , where study of inter-ELM turbulence in the pedestal region of H-mode plasmas in NSTX was carried out using this system and the parametric dependences of the spatial and temporal characeristics were compared with expectations for different types of turbulence.
- Here the radial coordinate is defined relative to the symmetry axis of the tokamak.
- Note that density turbulence has been observed using BES system on the NSTX device in the pedestal region of H-mode plasmas . Although the amplitudes reported there of are similar to those we observe in the peripheral region of L-mode plasmas, this turbulence is highly localised to the steep gradient region of the pedestal.
- A possible explanation is suggested by the the results of the linear stability calculations presented in Fig. 4 of Ref. , where it was found that introducing moderate levels of co-current flow shear into the global simulations modestly increased the maximum growth rates relative to the case without shear flow. It can be seen from Fig. 7 that at mid radius, so such an increase in growth rates there might explain the increased heat flux compared to case IV. The explanation of this effect suggested in Ref.  is that the sign of the toroidal flow determines whether the diamagnetic contribution to the net perpendicular flow shear either enhances or reduces that due to the sheared equilibrium flow, hence either stabilising or destabilising the turbulence .
- Similar values of were already reported in the study of inter-ELM, pedestal turbulence in NSTX  but in a very different turbulence regime to the L-mode plasmas studied here.
- Note, however, that the large spread of values of these timescales manifest in Fig. 10 means that, while they clearly are all of the same order, this limited dataset does not provide strong evidence that they are balanced in the same way as was found in Ref.  - there a much larger dataset was used and it was possible to infer that held even though local equilibrium parameters changed over a relatively broad range. What we see from the present dataset is in fact just sufficient to confirm that all the linear timescales remain roughly comparable to the general gyrokinetic sound time . The key finding is that the measured nonlinear timescales associated with the density fluctuations are much longer than - and so are the found in the simulations.
- Tang W M 2002 Phys. Plasmas 1856
- Rhodes T L et al2011 Nucl. Fusion, 063022
- Connor J W and Wilson H R 1994 Plasma Phys. Control. Fusion 719
- White A E et al2010 Phys. Plasmas, 056103
- Shafer M W et alPhys. Plasmas 2012 032504
- Casati A et al2009 Phys. Rev. Lett. 165005
- Bourdelle C et al2011 Nucl. Fusion 063037
- Ren Y et al2012 Phys. Plasmas 056125
- Mazzucato E et al2008 Phys. Rev. Lett. 075001
- Yuh H Y et al2011 Phys. Rev. Lett. 055003
- Smith D R et al2010 Rev. Sci. Inst. 10D717
- Smith D R et al2013 Phys. Plasmas 055903
- Holland C et al2008 J. Phys. Conf. Series 012043
- Holland C et al2009 Phys. Plasmas 052301
- Holland C et al2011 Phys. Plasmas 056113
- Petty C C et al2008 Phys. Plasmas 080501
- White A E et al2008 Phys. Plasmas 056116
- Shafer M W, et al2006 Rev. Sci. Inst. 10F110
- Lin Z and Hahm T S 2004 Phys. Plasmas 1099
- Diamond P H et al2001 Nucl. Fusion 1067
- Saarelma S et al2012 Plas. Phys. Control. Fusion 085012
- Field A R et al2012 Rev. Sci. Inst. 013508
- Joliet S, Bottino A and Angelino P 2007 Comp. Phys. Comm. 409
- Field A R et al2004 in Proc. 20th IAEA Fusion Energy Conf. (Vilamoura, Portugal, 2004) EXC/P2-11
- Field A R et al2011 Nucl. Fusion 063006
- Kotschenreuther M, Rewoldt G and Tang W M 1995 Comp. Phys. Comm. 128
- Hawryluk R J 1980 in An Empirical Approach to Tokamak Transport, ed. Coppi B, Vol. 1 CEC Brussels
- Ghim Y-c, Field A R Zoletnik S and Dunai D 2010 Rev. Sci. Inst. 10D713
- Field A R et al2012 in Proc. IAEA Fusion Energy Conf. (IAEA CN-197) San Diego USA EXC/P7-01
- Ghim Y-c et al2013 Phys. Rev. Lett. 145002
- Ghim Y c et al2013 Nucl. Fusion Lett. submitted [arXiv:1211.2883]
- Houlberg W A et al1997 Phys. Plasmas 3230
- McMillan B F et al2011 Phys. Plasmas 112503
- Roach C M et al2009 Plas. Phys. Control. Fusion 124020
- De Bock M F M et al2008 Rev. Sci. Inst. 10F524
- Conway N J et al2006 Rev. Sci. Inst. 10F131
- Scannell R et al2010 Rev. Sci. Inst. 10D520
- Lao L L et al1985 Nucl. Fusion 1611
- Fonck R J, Duperrex P A and Paul S F 1990 Rev. Sci. Inst. 3487
- Hutchinson I H 2002 Plas. Phys. Control. Fusion 71
- Chen L, White R B and Rosenbluth M N 1984 Phys. Rev. Lett. 1122-1125
- Zoletnik S et al1998 Plas. Phys. Control. Fusion 1399
- Smith D et al2013 Nucl. Fusion 113029
- Ghim Y-c et al2012 Plas. Phys. Control. Fusion 095012
- Durst R D et al1992 Rev. Sci. Instrum. 4907
- Hinton F L and Wong S K 1985 Phys. Fluids B 3082
- Barnes M, Parra F I and Schekochihin A A et al2011 Phys. Rev. Lett. 115003
- Takeshi I et al2008 Rev. Sci. Inst. 10F318
- Evensen H T et al1998 Nucl. Fusion 237
- Kishimoto Y et al1999 Plas. Phys. Control. Fusion A663
- Highcock E G et al2012 Phys. Rev. Lett. 265001
- Howard N T et alPhys. Plasmas 032510
- Jenko F April 2013 presentation at 10th ITPA Transport Topical Group Meeting, IPP Garching, Germany
- Hill P et al2012 Plas. Phys. Control. Fusion 065011
- Kim Y B, Diamond P H and Groebner R J 1991 Phys. Fluids B 8 2050