Nonlinear Flow Generation By Electrostatic Turbulence In Tokamaks

Nonlinear flow generation by electrostatic turbulence in tokamaks

Princeton Plasma Physics Laboratory report PPPL-4532, June 2010.

W. X. Wang, T. S. Hahm, S. Ethier, G. Rewoldt and W. M. Tang, Princeton University, Plasma Physics Laboratory, P.O. Box 451, Princeton, NJ 08543. P. H. Diamond, University of California, San Diego, La Jolla, California 92093.

Prepared for the U.S. Department of Energy under Contract DE-AC02-09CH11466.

PACS numbers: 52.25Fi, 52.35Ra, 52.65Tt

Abstract

Global gyrokinetic simulations have revealed an important nonlinear flow generation process due to the residual stress produced by electrostatic turbulence of ion temperature gradient (ITG) modes and trapped electron modes (TEM). In collisionless TEM (CTEM) turbulence, nonlinear residual stress generation by both the fluctuation intensity and the intensity gradient in the presence of broken symmetry in the parallel wave number spectrum is identified for the first time. Concerning the origin of the symmetry breaking, turbulence self-generated low frequency zonal flow shear has been identified to be a key, universal mechanism in various turbulence regimes. Simulations reported here also indicate the existence of other mechanisms beyond E cross B shear. The ITG turbulence driven "intrinsic" torque associated with residual stress is shown to increase close to linearly with the ion temperature gradient, in qualitative agreement with experimental observations in various devices. In CTEM dominated regimes, a net toroidal rotation is driven in the cocurrent direction by "intrinsic" torque, consistent with the experimental trend of observed intrinsic rotation. The finding of a "flow pinch" in CTEM turbulence may offer an interesting new insight into the underlying dynamics governing the radial penetration of modulated flows in perturbation experiments. Finally, simulations also reveal highly distinct phase space structures between CTEM and ITG turbulence driven momentum, energy and particle fluxes, elucidating the roles of resonant and non-resonant particles.

I. Introduction

Momentum transport and plasma flow generation are complex transport phenomena of great importance in magnetic confinement fusion. An optimized plasma flow is believed to play a critical role in both controlling large-scale (macroscopic) plasma stability and in reducing energy loss due to plasma microturbulence, and thereby achieving high quality performance in plasma confinement. The toroidal momentum transport has been observed to be highly anomalous in various magnetic fusion experiments, not only for its high level compared to the neoclassical value due to Coulomb collisions, but also for its highly pronounced non-diffusive and non-local nature. A striking finding is the observation of automatic toroidal rotation spin up in nearly all tokamaks, the so called intrinsic or spontaneous rotation, that is, toroidal plasmas can self-organize and develop rotation without an external momentum input. Nondiffusive phenomena can also exist in other transport channels such as energy and particle; however evidences so far in experiments appear not as intriguing as the intrinsic rotation. This phenomenon may play a critical role in determining plasma flows and, consequently, confinement performance, particularly in the International Thermonuclear Experimental Reactor (ITER). Note that intrinsic rotation in tokamaks is an example of a "negative viscosity phenomenon" in which an up-gradient component of the momentum flux organizes a structured mean flow. Negative viscosity phenomena are of broad interest in the context of atmospheres, oceans, stellar interiors, and other rotating fluids.

Understanding the momentum transport and flow generation is one of the highlighted issues of current fusion research. Out of various possible physical mechanisms governing plasma flow dynamics, the strong coupling between toroidal momentum and energy transport universally observed in fusion experiments suggests that micro-turbulence is a key player in determining plasma rotation, as well. The strong momentum-energy transport coupling via microturbulence as a "medium" was predicted by theory and observed in experiments about two decades ago. The strong coupling was also obtained by gyrokinetic simulations of ion temperature gradient turbulence over a wide range of plasma parameters.

For turbulence driven toroidal momentum flux, a generic structure can be expressed as three terms. The flux is proportional to minus the momentum diffusivity times the radial gradient of the toroidal rotation velocity, plus a pinch velocity times the rotation velocity, plus the residual stress.

In addition to diffusion (first term), there are two nondiffusive components, momentum pinch (second term) and residual stress (third term). The three components in the momentum flux are highly distinctive not only formally but also physically. Besides their different physical origins under turbulence circumstances, they have qualitatively distinct effects on the toroidal flow formation. The diffusive transport is well known in the direction opposite to the rotation gradient, leading to the relaxation of the rotation profile and the release of associated free energy. The momentum pinch term is a convective flux which is directly proportional to the rotation velocity, with the pinch velocity as the coefficient. Both momentum diffusion and pinch can move plasma mechanical momentum (that is, toroidal momentum carried by particles), and then rearrange the rotation profile, radially. A qualitative distinction is that momentum pinch can transport momentum in either direction, up-gradient or down-gradient.

The residual stress is defined as a specific part of the Reynolds stress, which depends directly on neither the rotation velocity nor its gradient. The residual stress has a fundamentally distinct effect on rotation profiles, and is shown to drive intrinsic rotation as a type of wave-driven flow phenomenon which operates via wave-particle momentum exchange. Obviously, it has no counterpart in the turbulence driven particle flux which, under the constraint of particle number conservation, consists of only diffusive (including sub- and super-diffusive) and convective components. On the other hand, the energy flux may contain a residual-stress-like component due to energy exchange between particles and waves. The residual stress can be shown in the momentum transport equation to be isomorphic in mathematical form to the integrated external momentum source which acts as a torque to drive the rotation. Thus, the residual stress can act as an internal local torque to spin up a plasma, offering an ideal mechanism to drive intrinsic rotation. For this reason, the divergence of the residual stress is widely referred to as the intrinsic torque in experimental and theoretical investigations. Note that all three components have been observed in tokamak experiments. Searching for nondiffusive elements and understanding underlying mechanisms have been the focus of recent intensive theoretical and experimental efforts.

In this paper, new results of non-diffusive toroidal momentum transport found from our global gyrokinetic simulations are reported. We focus our study on understanding the nonlinear residual stress generation and its effect on toroidal flow formation in electrostatic turbulence regimes of ion temperature gradient (ITG) modes and trapped electron modes (TEM). This study concerns a few critical issues which are highly relevant to experimental observations and theoretical studies. These include: i) mechanisms for turbulence driving residual stress; ii) mechanisms for breaking symmetry in the parallel wavenumber spectrum, where the symmetry breaking is shown to be a critical ingredient for turbulence to generate net acceleration of parallel (and toroidal) flows; iii) impacts of trapped electrons and electron turbulence on residual stress; iv) characteristic dependences of intrinsic rotation on plasma parameters and directional tendency of the rotation; and vi) the fraction of residual stress in the momentum flux. Remarkable results also include the finding of an interesting meso-scale phenomenon, "flow pinch", in collisionless TEM turbulence, which appears to phenomenologically reproduce the radial penetration of modulated flows demonstrated by perturbation experiments. Also presented are highly distinct phase space structures between TEM and ITG turbulence driven fluxes to elucidate the roles of resonant and non-resonant particles.

The remainder of this paper is organized as follows. In Sec. II, gyrokinetic simulation models employed in this work are described, and a benchmark study of a CTEM case is presented. In Sec. III, we discuss generic pictures of turbulence driven toroidal momentum flux obtained in our global simulations. We attempt to partition the momentum flux and calculate the fraction due to residual stress. We also examine the relationship between the momentum and the energy transport, calculating the intrinsic Prandtl number. In Sec. IV, we address the mechanism of nonlinear residual generation in ITG turbulence with focus on the effect of zonal flow shear on parallel wavenumber symmetry breaking. The parametric dependence of ITG driven intrinsic torque on the ion temperature gradient is explored in order to understand empirical trends observed in experiments. The key results of nonlinear residual stress and flow generation in CTEM turbulence and trapped electron effects in the ITG regime are presented in Sec. V. The role of both the turbulence intensity and the intensity gradient in driving residual stress is explored. Also discussed are highlighted meso-scale phenomena, particularly the flow pinch effect, in CTEM dominated regimes. In Sec. VI, the phase structures of momentum, energy and particle fluxes are presented with a lot of interesting details with regard to which and how particles contribute to plasma transport due to turbulence. Section VII presents conclusions.

II. Gyrokinetic simulation models of rotating plasma and treatment of kinetic electrons

In this work, our global turbulence simulation is carried out using the Gyrokinetic Tokamak Simulation (GTS) code. The GTS code is based on a generalized gyrokinetic simulation model using a delta-f particle-in-cell approach, and incorporates the comprehensive influence of non-circular cross section, realistic plasma profiles, plasma rotation, neoclassical (equilibrium) electric field, Coulomb collisions, and other features. It can directly read plasma profiles of temperature, density and toroidal angular velocity from the TRANSP experimental database, and a numerical magnetohydrodynamic (MHD) equilibrium reconstructed by MHD codes using TRANSP radial profiles of the total pressure and the parallel current (or safety factor), along with the plasma boundary shape.

First, we give a brief description of our gyrokinetic simulation model for rotating plasmas in this section. In a delta-f simulation, the turbulence fluctuations are considered as perturbations on top of the neoclassical equilibrium. The gyrokinetic particle distribution function is expressed as an equilibrium part plus a perturbation. The equilibrium distribution function of ions, with magnetic moment and parallel velocity as independent velocity variables, is determined by the neoclassical dynamics and obeys a drift-kinetic equation (Eq. 1 of the report) in which the time derivative of the equilibrium distribution, its advection by parallel motion, the equilibrium E cross B drift and the grad-B drift, and the parallel-gradient force associated with the magnetic moment and the equilibrium potential, are balanced against the ion Coulomb collision operator.

The lowest order solution of that equation is a shifted Maxwellian consistent with (large) plasma rotation (Eq. 2), in which the parallel flow velocity is related to the toroidal rotation by the toroidal angular velocity times the toroidal current divided by the magnetic field strength, and the ion density carries a poloidal variation associated with plasma rotation. The total equilibrium potential consists of two parts, a flux-surface-averaged part and a poloidally varying part. The poloidally varying component can be generated by the centrifugal force which drives charge separation on a magnetic surface in strongly rotating plasmas. Generally the radial potential is dominant. The equilibrium radial electric field can be calculated from a first-principles based particle simulation of neoclassical dynamics with important finite orbit effects, or obtained by direct experimental measurement if available.

Instead of using a true neoclassical equilibrium distribution function, which is unknown analytically, we use this lowest order solution for equilibrium toroidal plasmas in the present simulations. A shifted Maxwellian with either model or experimental profiles of density, ion temperature and toroidal angular velocity is prescribed for the ions. In the electrostatic limit, the ion gyrokinetic equation for the turbulence perturbed distribution of ion guiding centers (Eq. 3) carries, on its right hand side, the drive terms: the temperature-gradient and density-gradient drives advected by the fluctuating E cross B velocity, a Kelvin-Helmholtz-type drive proportional to the gradient of the parallel equilibrium flow, terms containing the parallel flow velocity that matter when the Mach number of the plasma flow is high, a mirror term, and the parallel acceleration by the gyro-averaged fluctuating potential; the linearized Coulomb collision operator closes it.

The GTS code solves the gyrokinetic Poisson equation in configuration space for the turbulence potential at the particle coordinates. Unlike in flux-tube or wedge codes, the real space, global Poisson solver, in principle, retains all toroidal modes from the zero-zero mode all the way to a limit which is set by grid resolution, and therefore retains full-channel nonlinear energy couplings. There are two largely different Poisson solvers implemented in the GTS simulation. In a simple geometry limit, that is, large aspect ratio and circular cross section, turbulence fluctuations on small spatial and fast time scales and the axisymmetric zonal flow on larger (meso-scale) spatial and slow time scales can be decoupled using i) a Pade approximation for the gyroaveraging operator, and ii) the assumption that the gyrokinetic double average of the flux-surface-averaged potential is the potential itself, so that flux-surface averaging and gyrokinetic double averaging commute. This results in two decoupled equations (Eqs. 4 and 5), one for the fluctuating potential in terms of the difference between the gyroaveraged ion and the electron density fluctuations, and one, a radial differential equation weighted by the flux-surface geometry, for the zonal potential.

Because turbulence dynamics on different spatio-temporal scales are separated in solving the Poisson equation, the advantages are apparent. However, the above approximations, particularly the second one, are not well justified in general toroidal geometry. This has motivated us to develop a generalized Poisson solver which solves an integral equation for the total potential, the sum of the fluctuating and zonal parts (Eq. 6).

While the adiabatic electron model has been widely used for simplicity in many earlier numerical and theoretical studies of ITG driven turbulence, non-adiabatic electron physics is in general irreducible in turbulence dynamics of toroidal systems. For ITG and TEM turbulence, where the perpendicular wavenumber times the electron gyroradius is much less than one, we use a drift kinetic description for electrons, neglecting the finite gyroradius effect. However, for electron gyroradius scale turbulence, such as electron temperature gradient driven turbulence, electrons are treated as fully gyrokinetic. Similarly as for ions, the delta-f method can be used to solve for the total perturbed electron guiding center distribution function corresponding to turbulence fluctuations. The equilibrium electron distribution satisfies the electron version of Eq. (1) and can be approximated by a shifted Maxwellian containing a parallel flow similar to that for the ions. Apparently, the perturbed electron distribution contains both adiabatic and non-adiabatic electron response. Another simulation model to treat kinetic electrons is to separate the non-adiabatic electron response explicitly, writing the electron distribution as the equilibrium, plus the adiabatic Boltzmann response, plus a non-adiabatic remainder, and to solve for the non-adiabatic part. In this case Eqs. (4) and (6) are replaced by versions (Eqs. 4b and 6b) in which the polarization term carries the factor one plus the ion-to-electron temperature ratio and the right hand side uses the non-adiabatic electron density fluctuation.

The left-hand side of the non-adiabatic electron equation contains a time derivative of the potential, which can easily give rise to numerical instability if it is calculated using direct finite differences. To avoid the numerical problem, a split-weight scheme was proposed, which uses a separate equation for calculating the time derivative of the potential. That equation, which is not a new equation, is obtained by taking the time derivative of the gyrokinetic Poisson equation and using the ion and electron continuity equations. In a toroidal system, the equation for the time derivative of the potential (Eq. 7) relates it to the difference of the divergences of the electron and ion particle fluxes plus a term involving the magnetic drift acting on the fluctuating potential and the equilibrium potential gradient. The fluctuating part of the time derivative is then obtained directly by subtracting its flux surface average.

It is noticed that many previous simulations include only trapped electrons for the non-adiabatic electron response. Numerically, the fast parallel streaming of passing electrons gives rise to a strict constraint on the time step size, which adds to the computational challenge. While the trapped electrons are the primary origin of non-adiabatic response, some passing electrons can be non-adiabatic too. In fact, dynamical division between trapped and passing electrons is, though not impossible, highly non-trivial during simulations because of the dependence of the trapping-passing boundary on the electric potential which evolves in time, and of the collisional trapping-detrapping process. Nevertheless, thanks to the availability of supercomputing capabilities, we retain full electron dynamics by including both trapped and untrapped electrons in the simulations.

The GTS simulation has been benchmarked against other gyrokinetic codes in the electrostatic regime and the large aspect ratio circular concentric geometry limit. Presented here are benchmark results of the trapped electron mode instability against the FULL code. The FULL code is a linear eigenvalue code, which calculates linear growth rates and real frequencies; it is radially local (corresponding to flux tube geometry), using the so-called ballooning representation. For this benchmark, an analytical equilibrium based on the so-called s-alpha model with alpha equal to zero is used in the FULL local calculations, and a corresponding numerical equilibrium is produced for GTS. The numerical equilibrium includes a small Shafranov shift due to non-zero plasma beta and higher order (in the small inverse aspect ratio) corrections, which are neglected in the analytical equilibrium.

The representative parameters used in the benchmark are: inverse aspect ratio 0.35; density and electron temperature gradient profiles peaked at a normalised value of 6.0 with a flat-topped radial envelope centred at half radius; an ion temperature gradient profile of unit amplitude with a Gaussian envelope; electron-to-ion temperature ratio 3; and a safety factor rising quadratically from 0.854. For local FULL simulations, the corresponding parameters used are: minor-to-major radius ratio 0.175, normalised density and electron temperature gradients 6.0, normalised ion temperature gradient 1.0, safety factor 1.4 and magnetic shear 0.78. For the benchmark, the GTS simulation, which is always global, is carried out in a radial domain from 0.1 to 0.9 in terms of normalized minor radius, and the TEM instabilities are measured at half radius.

Note that the magnetic axis is not included in the simulation domain in this study based on considerations from both the physical and the numerical aspect. First, numerical MHD equilibria expressed in magnetic coordinates, which is currently used by GTS, usually have insufficient resolution near the magnetic axis due to mapping from the original cylindrical coordinates. This may cause numerical problems when simulation particles get into the region. Furthermore, plasma profiles are usually flat in the region near the magnetic axis, which makes the region not essential for turbulence physics.

The linear benchmark results are presented in Fig. 1. The growth rates and the real frequencies from the global GTS calculation are slightly higher than the local eigenvalue FULL calculation. The overall difference is less than 10 percent. There are a few effects which may contribute to the difference. First, in particular, FULL is radially local, whereas GTS is radially global. Further, as a subtle detail in this benchmark simulation, GTS includes multiple (all) toroidal modes which start at very low level initially, so that interactions between the modes are negligible during the linear phase, whereas FULL calculates a single toroidal mode at each time. Finally, as previously reported, differences in magnetic geometry between the s-alpha model and an MHD equilibrium may contribute to discrepancies in gyrokinetic turbulence calculation results. It is also often observed in local linear calculations using the FULL code that the linear frequency and growth rate are rather sensitive to subtle differences of Shafranov shift and finite aspect ratio corrections in equilibrium. Taking into account all these distinctions between the two simulation models, the overall agreement is reasonable.

As a nonlinear benchmark effort, GTS global and GEM local simulations were carried out for electron-temperature-gradient-driven CTEM turbulence in a specific, experimentally relevant parameter regime. The purpose of these simulations is to verify the nonlinear generation of blob-like, large fluctuation structures with toroidal mode number of about 10 or below via dramatic inverse toroidal energy cascades, and they are beyond the scope of this paper and will be discussed elsewhere in a future publication.

Figure 1. Growth rate and real frequency for the TEM instability versus poloidal wave number, compared with the FULL code calculation.

Unlike ITG turbulence with adiabatic electrons, non-adiabatic electron dynamics can drive particle transport for turbulence such as TEM. It is well known that turbulence driven particle transport across the magnetic field lines is ambipolar, that is, flux-surface-averaged radial particle fluxes for electrons and ions are equal, so as to maintain the overall quasineutrality in a toroidal system. The ambipolarity property of TEM driven cross-field particle transport is tested in the GTS simulation. The time history of particle fluxes at 0.54 of the minor radius is plotted in Fig. 2, showing that electron and ion fluxes very closely track with each other all the time during the simulation. Moreover, the ambipolarity of turbulence driven particle transport is obtained locally over the entire radial domain of the global simulation, as is seen in the right panel of Fig. 2 which plots the steady state particle fluxes versus minor radius. This guarantees that quasineutrality is satisfied radially locally.

Figure 2. Time history of ion and electron particle fluxes (left) and steady state ion and electron particle fluxes versus minor radius (right), from the same simulation as Figure 1.

Global gyrokinetic turbulence is characterized by distinguishable dynamical phases in both coordinate space and wavenumber space. Ideally, the dynamics of gyrokinetic turbulence should be robust to numerical techniques. The robustness of turbulence dynamics with respect to different approaches for solving the gyrokinetic Poisson equation, the size of the simulation grids and the number of simulation particles, was carefully examined for ITG simulations previously. A further convergence study for CTEM turbulence is presented in Fig. 3. Two simulations using 50 and 100 particles per cell per species, respectively, are shown to produce well converged results for electron particle transport that displays no noticeable difference in a statistical sense (left panel). In other words, the difference in the simulated fluxes between the two cases is within the same range of statistical error as different simulation runs with the same number of particles but with different initial conditions. At the same time, the time evolution of the corresponding average electron weight square of the two simulations is shown to be almost identical (right panel).

Figure 3. Time history of electron particle fluxes (left) and average weight squares of electrons (right) from two simulations using different numbers of simulation particles. The major parameters are normalised electron temperature and density gradients of 6.5 and a normalised ion temperature gradient of 2.4, with a shaped DIII-D type MHD equilibrium.

This result indicates that the observed weight growth does not depend on whether 50 or 100 particles per cell per species are used in these simulations, and is driven by physics, corresponding to the increase of the amplitude of the perturbed distribution associated with plasma profile evolution induced by turbulence-driven fluxes during the transport time scale. Furthermore, while the particle weight is physically growing during the simulations, there is no observable correlation between the weight evolution and the dynamics of electron particle flux, as is shown in Fig. 3. The particle weight remains at a low level, with an average weight square below 0.09, at the end of the simulations, which does not impact the results of simulated transport. These convergence studies clearly indicate that the noise-induced transport in our simulations is negligible with respect to turbulence driven transport. For most of the simulations in this study, we use 100 particles per cell per species.

III. Momentum flux partition and Prandtl number

In this section, we first present generic pictures of turbulence driven toroidal momentum transport based on global gyrokinetic simulation results. We then discuss the relationship between the momentum and the energy transport. The critical quantity in the discussion is the Prandtl number, that is, the ratio of ion momentum and thermal diffusivities, which has attracted a lot of attention in experimental and theoretical studies. It is also highly interesting to examine the partition of the turbulence driven momentum flux, particularly the percentage of the non-diffusive component.

As we mentioned before, a generic structure for turbulence-driven toroidal momentum flux (Eq. 8) is the sum of three terms: minus the momentum diffusivity times the radial gradient of the rotation velocity, plus the pinch velocity times the rotation velocity, plus the residual stress. The three components, diffusion, momentum pinch and the residual stress, are highly distinctive not only formally but also physically, and have different effects on toroidal flow formation. The distinction among the three components, however, is a highly nontrivial task in practice, particularly in experiments. This formulaic difference can be used to design simulations and experiments for the identification and partition of the three components, as in our following numerical studies.

The relationship between turbulence-driven toroidal momentum and energy transport has long been an issue of interest in both experimental and theoretical investigations. Experiments on various machines have established a fairly comprehensive database over various regimes, including L-mode and H-mode plasmas, for the ratio of effective momentum and thermal diffusivity, also referred to as the raw Prandtl number. As a general validation against the experimental database, systematic simulations have been carried out, over a wide range of experimentally relevant plasma parameters, to investigate this issue with a focus on the Prandtl number. The results of this simulation study are summarized below.

At first, it is helpful to write down precisely the definitions used to calculate the relevant quantities in the simulations. The effective momentum diffusivity and the associated total (net) radial flux of toroidal momentum are defined (Eq. 9) by the velocity-space integral of the ion mass times the major radius times the toroidal velocity, weighted by the radial E cross B velocity and the perturbed ion distribution, set equal to minus the ion mass density times the effective momentum diffusivity times a geometric factor times the radial gradient of the toroidal angular velocity. This form is suitable for general geometry, with the radial coordinate denoting the magnetic flux surface.

In ion-dynamics-dominated regimes, our simulations verify that there exists strong coupling between ion momentum and heat transport for ITG driven turbulence, and the effective Prandtl number is on the order of unity. This is in broad agreement with a theoretically predicted trend and experimental observations in conventional tokamaks where low-wavenumber fluctuations are believed to be responsible for a high level of plasma transport. A typical simulation result is shown in Fig. 4 where an effective Prandtl number near unity is obtained in the long-time steady-state.

On the other hand, global gyrokinetic turbulence is characterized by distinguishable dynamical phases in both configuration space and wave-number space, and correspondingly the turbulence driven momentum transport can display different behavior over different dynamical phases of turbulence evolution and the Prandtl number varies. Particularly in the ITG turbulence regime with adiabatic electrons, a significant inward, non-diffusive momentum flux associated with residual stress is robustly observed in the post saturation phase, which is after the nonlinear saturation of the ITG instability but before a long term steady state. As is seen in Fig. 4, there is an up-gradient momentum flux generated in the post saturation phase. This non-diffusive flux can result in a great departure of the Prandtl number from unity. Moreover, the Prandtl number in the long-time steady state is also shown to vary over a certain range around unity, showing fairly sensitive dependence on plasma parameters. Finally, as a further validation effort, GTS simulation predictions of toroidal momentum and ion thermal transport have been directly compared with experimental measurements on DIII-D. Results for an ion-transport-dominated DIII-D discharge with relatively high toroidal rotation are presented in Fig. 5. Reasonably good agreement between the simulation of ITG turbulence and the experiment is obtained not only for the Prandtl number, but also for the individual values of the momentum and thermal diffusivities.

Figure 4. Time evolution of effective toroidal momentum and heat diffusivity (left) and the initial toroidal rotation profile (right).

Figure 5. Time history of ITG turbulence driven effective momentum and ion thermal diffusivity from a simulation of DIII-D discharge 125236 at 3.6 seconds, compared with experimental results from TRANSP analysis.

The calculation of the raw Prandtl number in experiments is relatively straightforward, and many experimental results reported are raw Prandtl numbers. While it is useful to look at the raw Prandtl number, a more meaningful physics quantity to examine is the ratio of pure momentum diffusivity and thermal diffusivity, referred to as the intrinsic Prandtl number, though it is harder to calculate. Obviously, the intrinsic and raw Prandtl numbers are equal if there is no non-diffusive contribution to the momentum flux, and the difference between the two reflects the fraction of non-diffusive contributions. Specifically, the ratio of non-diffusive to total momentum flux equals one minus the ratio of the intrinsic to the raw Prandtl number.

An intrinsic Prandtl number of unity was theoretically predicted for drift wave turbulence. Recent studies indicate departures from unity. To obtain the intrinsic Prandtl number, one has to separate out the nondiffusive components from the total momentum flux, which is a highly non-trivial task, particularly in experimental measurements. To this end, a set of numerical experiments have been carefully designed and carried out. In these simulations, toroidal rotation profiles are set to have zero rotation velocity at the normalized minor radius of 0.5, and to have a step-type profile for the rotation gradient — a flat-topped exponential envelope of index 6 and width 0.28 centred there. Here we use ITG turbulence with adiabatic electrons.

We first look at the marginally unstable ITG regime using two simulations with a normalised ion temperature gradient of 5.5 and equal electron and ion temperatures. For both simulations, only one parameter varies, the amplitude of the rotation gradient, which corresponds to the use of different initial rotation gradients. The simulation domain runs from 0.1 to 0.9 in normalised radius, and we focus on a narrow radial annulus centered at the point where the rotation velocity and the momentum pinch vanish. The remaining nondiffusive component of momentum flux in the region is the residual stress, which is independent of the variation of the rotation gradient. From Eqs. (8) and (9), the raw Prandtl number equals the intrinsic Prandtl number minus the residual stress divided by the product of a geometry factor, the ion thermal diffusivity and the rotation gradient.

The key idea for these simulations to allow for separation of the diffusive and non-diffusive components is based on the fact that the shear flow instability is very hard to drive unstable in a toroidal system because of the strong stabilization effect of the magnetic shear. However, the equilibrium (mean) E cross B shear flow, which is determined by neoclassical dynamics to relate to the toroidal rotation via the radial force balance relation, can influence both turbulence and residual stress generation. For simplicity, we exclude the equilibrium electric field in these simulations. In this case, we expect that rotation, particularly with relatively low gradient, has negligible effect on ITG driven turbulence, that is, fairly similar turbulence fields can be produced in the two cases with different rotations. Hence, it is reasonable to argue that the intrinsic Prandtl number and the ratio of residual stress to ion thermal diffusivity are held constant in these simulations. Note that turbulence intensity may slightly vary from one simulation to another; while both the momentum and thermal diffusivities are each roughly proportional to the intensity, their ratio is less sensitive to the intensity variation, as is the ratio of residual stress to thermal diffusivity. Therefore, from the two aforementioned simulations I and II, we can estimate the intrinsic Prandtl number (Eq. 10) as the difference of the two raw Prandtl numbers each weighted by its own rotation gradient, divided by the difference of the two rotation gradients.

Figure 6 (top panel) shows the time history of the raw Prandtl number, from linear phase to saturated steady state, of two simulations in which the second rotation gradient is twice the first. Note that a semi-quantitative definition should be used for "steady state", relative to linear growth and post-saturation phases, particularly in a global simulation. Specifically, the time averaged growth rate of the turbulence intensity vanishes, relative to the linear growth rate, though the instantaneous growth rate fluctuates. On the other hand, the toroidal spectra of fluctuations in the "steady state" turbulence regime, which can be used for comparison with experimental measurements in validation studies, are typically characterized by a significant down-shift from linearly unstable modes due to nonlinear toroidal energy cascades. However, a "steady state" may still exhibit considerable temporal variations in various transport quantities, particularly in a global simulation, due to meso-scale dynamics such as turbulence spreading, self-consistent plasma profile evolution, turbulence avalanches and so on.

As generally done for calculating transport fluxes in this type of turbulence simulation study, an averaged raw Prandtl number in the saturated turbulence steady state is calculated by time average. Generally, the interval of time averaging should be some time scale between the correlation time and the profile evolution time. In this case, an averaged raw Prandtl number is calculated over a period from t equal to 1000 to 2400, which spans many turbulence growth times. The obtained values are 0.961 and 0.703 in the two cases, with standard deviations of 0.19 and 0.14 respectively. These results, in terms of Eq. (10), give an estimate of 0.445 for the intrinsic Prandtl number in the marginal ITG regime. This result appears to be consistent with a recent theoretical prediction of an intrinsic Prandtl number of 0.2 to 0.5 in stiff profile regimes.

Plotted in the middle of Fig. 6 is a scan of the intrinsic Prandtl number versus the normalised ion temperature gradient, showing that the intrinsic Prandtl number increases with the temperature gradient. The ratio between the residual stress component and the total momentum flux is plotted at the bottom of Fig. 6, which shows that the residual stress contribution to the total momentum flux is significant — more than 50 percent for case I — and is increased with the decrease of the rotation gradient. This is readily understandable because the residual stress, unlike the diffusive component, is independent of rotation gradient. On the other hand, the fraction of residual stress does not show a clear, conclusive scaling trend with the ion temperature gradient.

Figure 6. Time evolution of the raw Prandtl number from two simulations with a normalised ion temperature gradient of 5.5 and different rotation gradients (top); intrinsic Prandtl number versus normalised ion temperature gradient (middle); and ratio of residual stress to total momentum flux versus rotation gradient and normalised ion temperature gradient (bottom).

IV. Nonlinear residual stress generation and the scaling of intrinsic torque

We have shown that the residual stress component can be a substantial portion of the total momentum flux driven by ITG turbulence. As discussed in the previous section, residual stress, acting as an internal torque, may play a critical role in driving intrinsic rotation. The micro-picture of this mechanism is that the net parallel (toroidal) flow is accelerated by turbulence, which critically depends on the details of the parallel wavenumber spectrum. For most drift wave instabilities, both signs of the parallel wavenumber are equally excited, resulting in a reflection symmetry in the spectrum. Perfect local symmetry means perfectly balanced population density between co- and counter-propagating acoustic waves and thus a vanishing net local momentum torque. Therefore, a critical, generic piece of physics behind the residual stress spinning up the plasma is the breaking of that symmetry.

Out of various theoretical possibilities, one of the leading candidates is a mean E cross B flow shear, which shifts the eigenmode to one side radially, and thus produces a non-vanishing spectrum-averaged parallel wavenumber. Another symmetry breaking mechanism, which leads to an inward pinch, can come from the interplay of magnetic field curvature and ballooning mode structure in toroidal geometry.

Recently, using global gyrokinetic simulation, a universal mechanism for parallel wavenumber symmetry breaking has been identified due to turbulence self-generated zonal flow shear, and an associated residual stress has been robustly observed in ITG simulations over a wide range of experimentally relevant parameters. From the viewpoint of local analysis and simulation, the turbulence self-generated zonal flow shear has no preferred direction in a long time statistical sense, and therefore was expected to have little direct effect on the parallel wavenumber spectrum. However, for global simulations, the zonal flow dynamics is found to be significantly different from the local picture. Specifically, zonal flow is shown to be slowly varying in time and of large scale in space. This is also an indication of the existence of toroidal zonal flow. A slowly varying large scale zonal flow structure has been clearly identified recently in drift wave turbulence in a linear machine. The observed low frequency, large scale zonal flow structure is shown to have a remarkable effect on the parallel spectrum of potential and density fluctuations.

The most straightforward way to examine turbulence driven residual stress is to set the initial rotation to be zero in simulations. Thus, only the residual stress component remains in the momentum flux. To elucidate the critical role of zonal flows in the nonlinear flow generation, unless explicitly specified, equilibrium E cross B flows are excluded in the following simulations. This may correspond to typical core turbulence apart from internal transport barriers and in L-mode plasmas, where equilibrium shear is not dominant.

The quantity used to characterize the symmetry breaking in our study is the average parallel wavenumber of the turbulence spectrum (Eq. 11), defined as the intensity-weighted average over poloidal and toroidal mode numbers of the quantity formed from the toroidal mode number times the safety factor minus the poloidal mode number, normalised by the safety factor and the major radius.

Figure 7 illustrates the simulation results of ITG turbulence with adiabatic electrons. For this case, the ITG instability is quite marginal, with a normalised ion temperature gradient of 4.9 and equal ion and electron temperatures. First, the upper-left panel of Fig. 7 shows that significant inward flux of toroidal momentum is driven in the whole radial range with ITG turbulence present. Particularly, a large inward momentum flux emerges in the post saturation phase, which is after the nonlinear saturation of the ITG instability but before a long term steady state. Because of the zero initial toroidal rotation used, the momentum flux is, by definition, essentially residual stress.

Plotted in the lower-left panel is the intensity-weighted parallel wavenumber sum, a quantity resembling the residual stress expression. One can see that it indeed reproduces a similar spatio-temporal behavior to the directly calculated momentum flux. Further, in the upper-right panel, the spectrum-averaged parallel wavenumber shows an apparent spatio-temporal correlation with the momentum flux, indicating the importance of the non-vanishing wavenumber. The whole picture for the residual stress generation is completed by finding out what causes the symmetry breaking, giving rise to the non-zero parallel wavenumber. This is in the lower-right panel, which plots the shearing rate of turbulence self-generated zonal flows according to a formula appropriate to shaped tokamak geometry (Eq. 12), built from the poloidal magnetic field, the total field, and the derivative with respect to poloidal flux of the zonal radial electric field divided by the major radius times the poloidal field.

A clear correlation between the zonal flow shearing rate and the average parallel wavenumber indicates that the breaking of symmetry and the yielding of a non-vanishing wavenumber are caused by the zonal flow shear. Since zonal flows are turbulence self-generated, this process represents a universal, nonlinear mechanism for residual stress generation. It is expected to play an important role in flow generation, particularly in L-mode plasmas where the E cross B shear of the equilibrium electric field is weak.

Figure 7. Spatio-temporal evolution of the radial flux of toroidal momentum (upper left), the intensity-weighted parallel wavenumber sum (lower left), the spectrum-averaged parallel wavenumber (upper right) and the zonal flow shearing rate (lower right).

As discussed previously, residual stress, acting like an intrinsic (internal) torque, is the only way to spin up a plasma from rest. The observation that an external torque in the counter-current direction is required to hold the plasma from rotating is another direct evidence of the existence of intrinsic torque. Recently, intensive experimental studies carried out on various machines attempted to identify the role of residual stress, and to characterize the dependence of the intrinsic rotation and intrinsic torque on plasma parameters. The empirical tendency obtained in H-mode plasmas shows that the offset value of the toroidal rotation typically scales with the increment in stored energy, and the rotation is usually in the co-current direction. Intrinsic rotations are observed to increase with increasing pressure gradient in various JT-60 plasmas.

The characteristic dependence of intrinsic torque driven by ITG turbulence is investigated using a set of systematic simulations. The radial profiles of ion temperature are given by specifying a temperature gradient profile as a flat-topped exponential envelope of index 6 centred at half radius, along with a fixed temperature of 1 keV there. This gives a fairly uniform ITG drive in a region centered at half radius, as shown in the top panel of Fig. 8. The temperature gradient amplitude varies from 4.9 to 8.2 for these simulations, covering a wide range from near to well beyond ITG marginality. The rest of the input parameters are the same for all six simulations.

The intrinsic momentum torque appearing in the toroidal momentum balance (transport) equation takes the form of the divergence of the residual stress. Instead of calculating that local torque, we examine the rate of toroidal momentum generation due to ITG turbulence. The residual stress is the only quantity responsible for the momentum build-up in this case. The mid-panel of Fig. 8 illustrates the spatio-temporal evolution of the flux-surface averaged toroidal momentum density, defined as the velocity-space integral of the ion mass times the major radius times the toroidal velocity weighted by the perturbed distribution. The quantity calculated here is the rate of total toroidal momentum generation, the time derivative of the volume integral of the momentum density magnitude. Apparently, this quantity is a measure of total (or spatially averaged) torque driven by turbulence, which has better correspondence to the intrinsic torque inferred from experiments.

As illustrated in Fig. 8 (bottom), the ITG driven intrinsic torque is shown to increase with the temperature gradient. A slightly stronger than linear scan of torque versus the normalised ion temperature gradient — and equivalently, torque versus the normalised pressure gradient, because of the fixed density profile used in all these simulations — is obtained. This result is consistent with experimental trends observed in various devices, including Alcator C-Mod where the central flow velocity scales linearly with the edge pressure gradient.

Figure 8. Radial profiles of the normalised ion temperature gradient, illustrated for the cases with amplitudes 5.5 and 7.6 (top); spatio-temporal evolution of the flux-surface averaged toroidal momentum density (middle); and total spatially averaged intrinsic torque versus the normalised ion temperature gradient (bottom).

The dominant underlying physics governing this characteristic dependence is that the residual stress is proportional to the turbulence intensity which, in turn, increases with the strength of the ITG drive. However, this does not explicitly give a linear dependence from a simple argument. The zonal flows and their effect on symmetry breaking, on the other hand, are also expected to increase with the increase of turbulence intensity. Therefore, we may expect a stronger than linear scan of torque versus the temperature gradient for the nonlinearly driven residual stress.

It is noticed that the dependence on pressure gradient can also be introduced via the equilibrium radial electric field (not included in these simulations) which relates to the pressure gradient through the well known radial force balance relation. However, this connection is less transparent. First of all, the equilibrium E cross B flow shear effects are twofold: reducing fluctuation intensity and breaking up the parallel wavenumber symmetry. Its overall effect on residual stress generation depends on the balance between the two. Secondly, the mean E cross B shear is proportional to both the second radial derivative of the pressure and the product of the density and pressure gradients. It should be pointed out that the scaling of the torque against the pressure gradient does not hold locally. For instance, one can have zero local torque, that is, the divergence of the residual stress is zero, at a location of strong temperature gradient and maximal residual stress.

The nonlinear residual stress generation is also observed for electron driven turbulence such as CTEM. We will present these results in the next section. The characteristic dependence of the associated residual-stress-driven torque on gradients of electron profiles, such as the electron temperature, density and pressure gradients, can be established through electron driven turbulence. Even more complex connections are expected between the intrinsic torque and plasma profiles in the presence of hybrid ITG and TEM turbulence, which is more likely to be the case in experiments. However, all these are beyond the scope of this paper and will be discussed elsewhere in future publications.

V. Nonlinear residual stress generation in CTEM turbulence and trapped electron effects in the ITG regime

In this section, we discuss the effects of trapped electrons, focusing on nonlinear residual stress generation by TEM turbulence and ITG turbulence with non-adiabatic electrons. For simulations with kinetic electrons presented in this section and hereafter, unless explicitly specified, the working gas is hydrogen, that is, the ion-to-electron mass ratio is 1836.

We first examine ITG turbulence. The major parameters used are: normalised ion temperature gradient 5.3, normalised electron temperature gradient 1.6, normalised density gradient of about 1 or below, and initial rotation zero. For these parameters, the ITG modes are marginally unstable, as found in many experiments, and TEM modes are stable. This allows us to investigate the same turbulence, that is ITG, when we switch the electron response in the simulations from adiabatic to non-adiabatic. Figure 9 shows the results of ITG turbulence with adiabatic electrons. Similar to Fig. 7, close spatio-temporal correlations among the momentum flux, the intensity-weighted parallel wavenumber sum, the spectrum-averaged parallel wavenumber and the zonal flow shearing rate illustrated in Fig. 9 clearly demonstrate that the residual stress is nonlinearly driven by the fluctuation intensity, acting with the zonal flow shear which induces symmetry breaking in the parallel wavenumber spectrum.

Figure 9. Spatio-temporal evolution of the radial flux of toroidal momentum (upper left), the intensity-weighted parallel wavenumber sum (lower left), the spectrum-averaged parallel wavenumber (upper right) and the zonal flow shearing rate (lower right), from an ITG simulation with adiabatic electrons.

The results for ITG turbulence with non-adiabatic electrons are presented in Fig. 10, which uses exactly the same set of simulation parameters as in Fig. 9. An immediate observation is that the spatio-temporal correlations among the plotted four quantities become obviously less clear, compared to the ITG case with adiabatic electrons. First, the lesser similarity in spatio-temporal structures between the momentum flux (residual stress, upper-left panel) directly calculated from Eq. (9) and the estimate from the intensity-weighted parallel wavenumber sum (lower left) indicates that the turbulence intensity driven residual stress does not fully account for the residual stress produced by the turbulence. Further, less correlation in spatio-temporal structures between the spectrum-averaged parallel wavenumber (upper right) and the zonal flow shearing rate (lower right) indicates that the zonal flow shear does not fully account for the origin of the non-vanishing wavenumber.

The non-adiabatic electrons are shown to introduce finer radial scales into the zonal flows. In configuration space, this appears as small wiggles — finer structures with small amplitude — sitting on a large scale structure with large amplitude. The corresponding E cross B shear at small scales and low frequencies, however, appears too weak to have a visible impact on the parallel wavenumber spectrum. The key points made by these results clearly indicate: i) the existence of other possibilities for driving residual stress, and ii) the existence of other mechanisms beyond E cross B shear for parallel wavenumber symmetry breaking. For the former, one interesting candidate is the turbulence intensity gradient, whose important role will be elucidated in our CTEM simulations to be presented later. For the latter, the possible mechanisms include the effects of magnetic shear, nonlinear mode couplings, and turbulent radial current, which will be addressed in a future publication.

Figure 10. Spatio-temporal evolution of the radial flux of toroidal momentum (upper left), the intensity-weighted parallel wavenumber sum (lower left), the spectrum-averaged parallel wavenumber (upper right) and the zonal flow shearing rate (lower right), from an ITG simulation with non-adiabatic electrons, using the same parameters as Figure 9.

Now we present a simulation of an experimental case, which shows trapped electron effects on ion turbulence and transport. Simulation results presented in Fig. 11 are for a DIII-D experiment. This is an ion transport dominated DIII-D discharge with low toroidal rotation. A relatively large ion temperature gradient exists in the range of normalised radius 0.2 to 0.5, which makes ITG modes unstable, while TEMs are stable for most minor radii. The real mass ratio of 3672 for a deuterium plasma is used in this simulation.

As is illustrated in the left panel of Fig. 11, the ITG turbulence with adiabatic electrons is shown to produce ion heat transport at a much higher level than the neoclassical in the inner core area, which matches the experimental level in the region. However, the ITG simulation with adiabatic electrons fails to account for the observed high level ion transport in the outer core region beyond about 0.45 of the minor radius, where the ITG instability is marginal or even stable. Trapped electron physics is found to play a critical role in this region. When trapped electrons are included, they substantially destabilize the ITG mode due to the resonance between the mode frequency and the toroidal precession frequency, with a corresponding change in the electron response as compared to the adiabatic case. This resonance occurs for precession drift-reversed trapped electrons, and has a dependence on the magnetic shear. The net effect of trapped electrons is to increase the linear growth rate and the nonlinear saturation level. Consequently, ITG driven fluctuation intensity is substantially enhanced, particularly in the outer core region where the pure ITG modes are marginal or stable. As a result, the simulated ion heat flux is increased to be closer to the experimental observations in the region, while the ion transport is not considerably affected in the inner core region where the ITGs are strongly unstable.

In the two cases, the core ITG turbulence cannot reproduce the experimental level of ion heat transport in the further outer core region beyond 0.6 of the minor radius, where ion transport may be largely controlled by edge-core coupling. On the other hand, the ITG driven toroidal momentum flux is also substantially increased by non-adiabatic electrons (the right panel of Fig. 11). For this DIII-D shot, the toroidal rotation is small and flat in the region beyond 0.4 of the minor radius due to the use of counter neutral beam injection, which balances the intrinsic torque. This implies that the momentum flux observed in the simulation mostly comes from the residual stress.

Figure 11. Ion heat fluxes versus minor radius from ITG simulations with adiabatic and non-adiabatic electrons, compared with the experimental result from TRANSP and the neoclassical level from GTC-Neo, for DIII-D discharge 129533 at 4.075 seconds (left); and time history of ITG driven toroidal momentum fluxes, showing the enhancement due to trapped electrons (right).

Now we turn to discussing one of the key results of our simulation study, that is, nonlinear residual stress and flow generation in CTEM turbulence. The CTEM simulation presented below employs typical parameters of DIII-D plasmas. The major parameters used here are: normalised electron temperature and density gradients both 6.0, normalised ion temperature gradient 2.4, electron temperature 4.8 keV and ion temperature 3.5 keV at half radius, and an initial rotation of zero. A numerical MHD equilibrium corresponding to a real DIII-D discharge is used. An equilibrium electric field, which satisfies the radial force balance relation, is also included in this CTEM simulation.

As a key player in drift wave turbulence in general and in residual stress generation specifically, zonal flows generated by global CTEM turbulence display distinct characteristics compared to the large scale, stationary ones typically observed in global ITG turbulence. As illustrated in Fig. 12, in addition to radially global structures with near-zero radial wavenumber which are dominant, there are significant shorter scale structures at a normalised radial wavenumber of about 0.2 and even fine but weak structures at about 0.6. Recent theoretical calculations of zonal flow growth rate indicate zonal flow generation at fine scales, with radial wavenumber times ion gyroradius of order unity, in CTEM turbulence. The zonal flows are also shown to be less stationary, exhibiting significant temporal variation in amplitude. In the frequency domain, the zonal flows peak at zero frequency, but with a certain extension to the low frequency range. At the same time, the zonal flow component at the geodesic acoustic frequency is very weak. The zonal flow shearing rate, however, is high, which is shown to have a strong effect on the turbulence parallel wavenumber spectrum.

Figure 12. Spatio-temporal evolution of zonal flows generated in CTEM turbulence (left) and the corresponding spectra in frequency and radial wavenumber space (right).

The nonlinear residual stress generation by CTEM turbulence, for the first time, is clearly observed in global simulations, as illustrated in Fig. 13. First, the CTEM-driven residual stress exhibits coherent spatio-temporal bursting behavior with momentum flux pulses propagating both inward and outward in the radial direction, as shown in the top-left panel of Fig. 13. The residual stress at steady state changes direction from outward in the inner core region to inward in the outer core region. The mid-left panel of Fig. 13 is the intensity-weighted parallel wavenumber sum, which represents the component of the residual stress driven by the turbulence intensity in the presence of a non-vanishing parallel wavenumber. The observation of a clear correlation between the residual stress and that quantity indicates that the turbulence intensity plays a major role in driving the residual stress, particularly in the inner core region inside 0.55 of the minor radius. In the outer core region, however, the turbulence intensity effect appears not to account well for the residual stress generation.

Again, the strong correlation between the spectrum-averaged parallel wavenumber (top right) and the zonal flow shearing rate (mid right) shown in Fig. 13 elucidates that the CTEM self-generated zonal flow shear plays a key role in breaking the parallel wavenumber symmetry. It is interesting to compare the effect of equilibrium shear, which is included in this simulation via the radial force balance relation, and the effect of zonal flow shear. The total E cross B shear rate, zonal flows plus equilibrium flow, is plotted in the bottom panel, showing very similar structures to those of pure zonal flows. This is because the equilibrium shear is much weaker than the zonal flows. In this particular case, the equilibrium shear is shown to have a minor effect on symmetry breaking, and thus on residual stress generation. However, the equilibrium flow shear is expected to have a significant effect on turbulence driven residual stress in regions of transport barriers, both in the core and at the edge.

Figure 13. Spatio-temporal evolution of the radial flux of toroidal momentum (top left), the intensity-weighted parallel wavenumber sum (mid left), the spectrum-averaged parallel wavenumber (top right), the zonal flow shearing rate (mid right) and the total E cross B shearing rate (bottom), from the simulation of CTEM turbulence.

It has been remarked that the mechanism of turbulence intensity does not fully explain the residual stress generation, particularly in the outer core region beyond about 0.55 of the minor radius. It has been suggested in theory that turbulence intensity gradients can also contribute to driving residual stress, in addition to the turbulence intensity. The intensity gradients in both the radial direction and the radial wavenumber direction may act to drive a residual stress. Here, we focus our discussion on the role of the intensity gradient in the radial direction. Figure 14 shows the spatio-temporal evolution of the quantity formed from minus the average parallel wavenumber times the radial derivative of the fluctuation intensity, which can be used to approximately represent the intensity gradient driven residual stress. Its apparent correlation with the directly calculated residual stress (top-left panel of Fig. 13) is noted particularly in the outer core region where a significant effect of the turbulence intensity is not observed. This simulation result is a clear identification of the turbulence intensity gradient driving residual stress in the presence of zonal flow shear induced symmetry breaking.

Figure 14. Spatio-temporal evolution of the average parallel wavenumber times the radial derivative of the fluctuation intensity, which represents the part of the residual stress driven by the turbulence intensity gradient.

A few highly remarkable, interesting features revealed in the CTEM simulation are worth further discussion. First, the CTEM-driven residual stress changes sign from outward in the inner core region to inward in the outer core region, as is more clearly seen in Fig. 15 (upper-left panel) which shows the radial profile of the residual stress, time averaged over the steady state. What determines the sign of the residual stress, particularly its relation with plasma parameters, remains an important issue.

The residual stress is shown to act as an intrinsic torque effectively. The resultant parallel flow (or toroidal flow) generation process is demonstrated in the lower panels of Fig. 15. In this case, a parallel flow is driven mostly in the counter-magnetic-field direction in the whole region from 0.25 to 0.8 of the minor radius where CTEM turbulence is excited. The corresponding toroidal rotation is in the co-current direction for this DIII-D geometry case, which is in agreement with the experimental trend of intrinsic rotation observed in various tokamaks. The maximum parallel flow velocity generated at the end of the simulation reaches about 5 percent of the local ion thermal velocity.

Further, the CTEM turbulence and transport are characterized by burstings which emerge regularly in time and propagate radially. The coherent spatio-temporal bursting phenomenon was observed in ITG simulations, but appears more pronounced in the TEM turbulence regime. The bursting generation frequency and the radial propagation velocity are estimated to be about a tenth of the sound speed divided by the minor radius, and about five to ten thousandths of the sound speed, respectively. More interestingly, it is found that the temporal burstings and radial propagation are also directly displayed in the parallel flow during its generation process, as is clearly seen in the lower panels of Fig. 15. Particularly, it is shown that small parallel flow perturbations are generated locally, in the center of the plasma in the simulation case, by the turbulence, and then propagate radially. The measured propagation velocity is about seven thousandths of the sound speed.

This "flow pinch" phenomenon revealed in the simulations may have analogues in experiments. The radially inward propagation of toroidal flow perturbations generated by modulated beams in the peripheral region near the plasma edge was demonstrated in JT-60U experiments, which was attributed to the turbulence driven momentum pinch. Nevertheless, our simulation results of flow pinch may offer a new insight into the underlying dynamics governing the radial penetration of localized modulated flows in perturbation experiments. We should particularly point out that the meso-scale phenomena and their critical role in determining plasma transport and its radially nonlocal nature are highly pronounced in the TEM turbulence regime.

Figure 15. Averaged momentum flux (residual stress) at steady state versus minor radius (upper left); radial profile of the ion parallel flow, in the counter-field and therefore co-current direction, at the end of the simulation (upper right); spatio-temporal evolution of the parallel flow (lower right); and time history of the parallel flow at three radial locations (lower left), which more clearly illustrates the flow pinch phenomenon.

VI. Phase space structures of turbulence driven fluxes

It is highly interesting and instructive to examine the phase space structures of various turbulence driven fluxes. This type of study can provide physical pictures at a very fundamental level with regard to which and how particles contribute to plasma transport due to turbulence. Particularly, this can help elucidate the roles of resonant and non-resonant particles.

In Fig. 4, it is observed that the ITG driven toroidal momentum flux reverses its sign from inward in the post saturation phase to outward in the long-time steady state. The phase space structures of the momentum flux at the two stages are presented in Fig. 16. To be precise, illustrated in Fig. 16 is the momentum flux resolved in radius, parallel velocity and perpendicular velocity, which integrates over the two velocity variables to give the momentum flux defined in Eq. (9) without a flux surface average, and is calculated at the midplane at 0.54 of the minor radius. A similar definition applies to the particle and energy fluxes whose structures are also discussed in this section.

First, it is observed that the momentum flux is carried mostly by passing and barely trapped ions. There are four signed peaks which are regularly located in the velocity plane, indicating dominant contributions from four different particle groups. The four groups of particles are distinguished by the amplitude of energy, low or high, and the sign of the parallel velocity, positive or negative, and contribute to the momentum flux in different ways. Specifically, one higher energy group with positive parallel velocity and one lower energy group with negative parallel velocity make positive contributions to the momentum flux, that is, outward momentum flux; another higher energy group with negative parallel velocity and another lower energy group with positive parallel velocity make negative contributions, that is, inward momentum flux. On the other hand, contributions from high energy ions with velocities beyond three thermal velocities are small.

Note that these characteristic phase space structures persistently appear in both the post-saturation stage and the long-time steady state with no considerable difference, while the net momentum fluxes at the two stages are in opposite directions. As for what makes the total momentum flux inward or outward, the simulation results suggest that it depends on the relative amount of each particle group's contribution, which may have to do with the details of the turbulence spectrum.

Figure 16. Phase space dependence of the ITG turbulence driven ion toroidal momentum flux at the long-time steady state (top) and at the post saturation phase (bottom), corresponding to Figure 4. The straight lines denote the boundaries between trapped and passing particles.

It is highly instructive to compare phase space structures in different transport channels. As illustrated in Fig. 17, the ITG driven heat flux is carried by different ions in a different way. First, the heat flux is carried mostly by trapped and barely trapped ions. Higher energy, mostly trapped, ions make the largest contribution, which is positive and peaked at about 2.5 thermal velocities; lower energy ions around the trapped-passing boundaries, centered at about one thermal velocity, make a negative contribution. These features share similarity to some extent with neoclassical transport in the collisionless regime, which may imply that the ion transport driven by the fluid-type ITG turbulence is non-resonance dominated.

Figure 17. Phase space dependence of the ITG turbulence driven ion heat flux at the steady state, corresponding to Figure 4.

We have shown in Sec. V that the trapped electron physics has a strong impact on ITG turbulence, particularly in a regime close to or below the ITG marginality. Particularly, the non-adiabatic electrons are shown to substantially enhance residual stress generation. It is interesting to understand how the inclusion of trapped electrons could change the way ions interact with turbulence and thus carry the momentum flux. The phase space dependence of the momentum flux is compared between the ITG turbulence with adiabatic and non-adiabatic electrons in Fig. 18, which is obtained from the same simulation of the DIII-D discharge as in Fig. 11. It is shown that the ITG driven momentum flux in experimental conditions possesses all the characteristics described previously for Fig. 16 for the case of large aspect ratio circular geometry and model plasma profiles. More importantly, trapped electrons are found not to change the qualitative phase space structure of ITG driven momentum and heat fluxes. The enhancement in the momentum flux production due to trapped electrons in ITG turbulence is mainly related to the increase of turbulence intensity.

Figure 18. Phase space structures of the toroidal momentum flux driven by ITG turbulence with kinetic electrons (top) and adiabatic electrons (bottom), from the same simulation as Figure 11.

For typical plasma parameters of fusion experiments, collisionless TEM turbulence can be a major source to drive multiple-channel transport. Now we turn to discuss the phase space characteristics of plasma transport produced by CTEM turbulence. First, the TEM turbulence driven momentum flux shows a highly distinct topology in phase space structures compared to that of ITG turbulence. Figure 19 shows the CTEM simulation results of the residual stress component. For both electron-temperature-gradient-driven and density-gradient-driven CTEM turbulence, the ion momentum flux of residual stress is carried by two groups of ions which are divided by the sign of the parallel velocity: ions with positive parallel velocity carry outward flux and ions with negative parallel velocity carry inward flux. The dominant contributions come from passing ions at around one perpendicular thermal velocity and about one and a half to two parallel thermal velocities of either sign. While the phase space structures between the two drive cases look qualitatively similar, it is also possible to notice a subtle difference: in electron-temperature-gradient-driven CTEM turbulence, the trapped ion region is basically a null region for the momentum flux.

Figure 19. Phase space structures of the momentum flux in electron-temperature-gradient-driven (top) and density-gradient-driven (bottom) CTEM turbulence.

Compared to the ITG case, the CTEM turbulence driven ion energy transport is also caused by ions from different regions and in a different way. As shown in the lower panel of Fig. 20, the dominant contributions come from two groups of passing ions, both carrying positive, outward ion flux. On the other hand, the electron responses in CTEM turbulence are shown to concentrate sharply in the trapped region, clearly spelling out the effect of trapped electron modes. As illustrated in the upper panel of Fig. 20, an outward flux of energy is carried by trapped electrons with higher energy, with little contribution from passing electrons. At the same time, deeply trapped, low energy electrons are shown to carry an inward, but small, flux for energy. Moreover, the electron phase space structures appear symmetric around zero parallel velocity. Apparently, the electron transport is dominated by the precession drift resonance of trapped electrons.

Figure 20. Phase space structures of the electron (top) and ion (bottom) energy flux in CTEM turbulence.

Remarkably, our simulation results clearly reveal that the particle flux is carried by the same particles in the phase space as the energy flux. This result holds true for both ions and electrons at different turbulence regimes. The results of CTEM driven ion and electron flux are presented in Fig. 21, which displays high similarity to the energy fluxes in Fig. 20.

Figure 21. Phase space structures of the electron (top) and ion (bottom) particle flux in CTEM turbulence, from the same simulation as Figure 20.

In the regime of ITG turbulence with adiabatic electrons, Fig. 22 illustrates how a net, outward heat flux and a vanishing ion particle flux can be produced by the same groups of ions. Higher energy ions, mostly trapped ones, carry outward fluxes for both heat and particles; and lower energy ions, mostly barely trapped ones, carry inward fluxes for both heat and particles. For the heat flux, the outward flux exceeds the inward, and a net outward flux remains. For the particle flux, however, the inward and outward components are balanced with each other, resulting in a vanishing net flux.

Figure 22. Phase space structures of the heat flux (top) and particle flux (bottom) in ITG turbulence with adiabatic electrons.

VII. Conclusions

Global gyrokinetic simulations using experimentally relevant parameters have revealed an important nonlinear flow generation process due to the residual stress produced by electrostatic turbulence of ion temperature gradient modes and trapped electron modes. Turbulence self-generated low frequency zonal flow shear has been identified to be a key, universal mechanism in various turbulence regimes for parallel wavenumber symmetry breaking, which is a critical ingredient for parallel (and toroidal) flow generation by turbulence. The principal results of this study are summarized as follows.

i) The nonlinear residual stress generation has been clearly observed, for the first time, in CTEM turbulence. Particularly, in addition to turbulence fluctuation intensity driving residual stress, which was also reported previously for ITG turbulence, the intensity gradient is also identified to drive significant residual stress, by acting with the CTEM self-generated zonal flow shear which induces symmetry breaking in the parallel wavenumber spectrum.

ii) The residual stress, acting as an intrinsic torque, is shown to spin up toroidal rotation effectively. In the simulated CTEM case with typical DIII-D parameters, where the plasma is initially rotation-free and momentum-source-free, a net toroidal rotation is produced in the co-current direction in the whole turbulence region. This is consistent with the experimental trend of observed intrinsic rotation. The total toroidal momentum is generated at an approximately constant rate, with the peak of the corresponding parallel flow profile at the end of the simulation reaching about 5 percent of the local ion thermal velocity. This net, directional mechanical flow generation phenomenon is an indication of momentum transfer from turbulence waves to particles via residual stress.

iii) The CTEM turbulence and transport including the momentum flux are characterized by burstings which emerge regularly in time and propagate radially both inward and outward. The meso-scale phenomena appear more pronounced than in the ITG turbulence regime, and are found to play a critical role in determining plasma transport and its radially nonlocal nature. One highly remarkable result is the observation of the "flow pinch" phenomenon. Specifically, toroidal flow perturbations, which are generated locally by the turbulence in the center of the plasma in the simulation case, are found to propagate radially. This result may offer an interesting new insight into the experimental phenomenon of radially inward penetration of perturbed flows created by modulated beams in peripheral regions.

iv) In the ITG turbulence regime, the intrinsic torque associated with residual stress is predicted to increase close to linearly with the value of the temperature gradient offset from the ITG critical gradient. The dominant underlying physics governing this scaling is that both the residual stress and the zonal flow shear increase with the turbulence intensity which, in turn, increases with the strength of the ITG drive. This simulation result is in qualitative agreement with experimental trends observed in various devices, such as the Rice scaling in which the increment of central toroidal flow velocity for H-mode plasmas scales linearly with the increment in the plasma stored energy divided by the plasma current.

v) For typical tokamak parameters, the nonlinearly generated residual stress is found to contribute up to more than 50 percent of the total momentum flux produced by ITG turbulence. It is plausible that the portion of the residual stress increases with the decrease of the rotation gradient. The intrinsic Prandtl number is shown to increase with the ion temperature gradient, specifically ranging from about 0.4 to 0.7 for normalised ion temperature gradients from 5.5 to 8.5. This result is in general agreement with observations in NSTX, where an estimated intrinsic Prandtl number of 0.5 to 0.8 was reported from the experimental database of various shots.

vi) While the critical effect of zonal flow shear on the parallel wavenumber spectrum is clarified, our simulations, particularly with electron physics included, also indicate the existence of other mechanisms beyond E cross B shear for symmetry breaking. The possible mechanisms include the effects of magnetic shear, nonlinear mode couplings, and turbulent radial current, which will be addressed in a future publication.

vii) Our simulations reveal highly distinct phase space structures between ITG and TEM turbulence for momentum, energy and particle fluxes, with a lot of interesting details with regard to which and how particles contribute to ion and electron transport in different channels. This study can ultimately help elucidate the roles of resonant and non-resonant particles in plasma transport in different turbulence regimes, which is a highly non-trivial issue under turbulence circumstances with many modes nonlinearly coupled together.

viii) In the ITG marginality regime, trapped electron physics is shown to play a critical role in determining plasma transport, not only producing the proper ion heat flux in experiments but also largely enhancing the residual stress generation. However, trapped electrons do not change the qualitative phase space structure of ITG driven momentum and heat fluxes.

Acknowledgments

We would like to acknowledge useful discussions with Drs. S. M. Kaye, J. Lang, W. W. Lee, C. J. McDevitt and W. Solomon. This work was supported by U.S. DOE Contract No. DE-AC02-09CH11466 and the SciDAC project for Gyrokinetic Particle Simulation of Turbulent Transport in Burning Plasmas. Simulations were performed at the National Energy Research Scientific Computing Center (NERSC) and the National Center for Computational Sciences (NCCS).

(The reference list, the twenty-two figures and the original equations are in the complete report at osti.gov/servlets/purl/984349. On this site, the alpha-particle-driven Alfvén eigenmodes measured in the same laboratory's TFTR tokamak are at /library/stm-acf328b26d, the field-reversed configuration — a plasma that generates and holds its own confining field — is at /library/stm-1aa14f2192, and the dense plasma focus, the laboratory's most heavily instrumented self-organising pinch, is at /library/stm-05100e66da.)