Microturbulence in magnetic confined plasmas contributes to energy exchange between particles of different species as well as the particle and heat fluxes. Although the effect of turbulent energy exchange has not been considered significant in previous studies, it is anticipated to have a greater impact than collisional energy exchange in low collisional plasmas such as those in future fusion reactors. In this study, gyrokinetic simulations are performed to evaluate the energy exchange due to ion temperature gradient (ITG) turbulence in a tokamak configuration. The energy exchange due to the ITG turbulence mainly consists of the cooling of ions in the B-curvature drift motion and the heating of electrons streaming along a field line. It is found that the ITG turbulence transfers energy from ions to electrons regardless of whether the ions or electrons are hotter, which is in marked contrast to the energy transfer by Coulomb collisions. This implies that the ITG turbulence should be suppressed from the viewpoint of sustaining the high ion temperature required for fusion reactions since it prevents energy transfer from alpha-heated electrons to ions as well as enhancing ion heat transport toward the outside of the reactor. Furthermore, linear and nonlinear simulation analyses confirm the feasibility of quasilinear modeling for predicting the turbulent energy exchange in addition to the particle and heat fluxes.

Numerous studies on anomalous transport of particles and heat generated by microscopic turbulence have been done, based on gyrokinetic theory and simulation.1–8 However, there have been fewer theoretical works on energy exchange between different particle species due to turbulence,9–12 and a small number of simulations have been performed to investigate the effect of turbulent energy exchange on the evolution of temperature profiles.13 In Ref. 13, it is shown from gyrokinetic simulations that the effect of turbulent energy exchange is negligibly small under conditions of DIII-D shot 128 913. However, sufficient comparative studies have not been made between collisional and turbulent energy exchanges in a wide range of conditions. In particular, in the case of high temperature plasmas, the impact of collisional energy exchange is expected to be small due to the low collision frequency, while turbulent energy exchange can work actively even in collisionless plasmas. As examples of other works related to turbulent energy exchange, there have been studies on the scaling of the turbulent transport and heating of impurities in magnetized plasmas14 and the thermal equilibration of ions and electrons by turbulence in astrophysical plasmas.15 

Since collisional heat transfer from alpha-heated electrons to ions is expected to play a critical role in sustaining burning plasmas in future reactors, it is an important issue to compare the effects of turbulence and collisions on the energy exchange between ions and electrons in high temperature plasmas. In the present paper, we evaluate the energy exchange in ion temperature gradient (ITG) turbulence by gyrokinetic turbulence simulations and investigate its properties, such as dependence on the ratio between ion and electron temperatures in comparison to those of the collisional energy exchange.

To perform simulations that predict global density and temperature profiles, it is practical to use turbulent transport models, such as quasilinear ones,16–21 rather than running direct turbulence simulations for all cases. In fact, it has been shown that these models can reproduce particle and heat fluxes obtained from gyrokinetic turbulence simulations within acceptable errors. A quasilinear model for turbulent energy exchange is shown in Ref. 22, where an electron drift wave instability is treated in a slab geometry with no temperature gradient. Detailed studies have not been conducted on modeling turbulent energy exchange under more complex conditions such as those for toroidal ITG mode remains. Quasilinear models are based on the assumption that the ratios of turbulent transport fluxes and energy exchange to the squared potential fluctuation amplitude estimated by linear analyses, which are called quasilinear weights,18 take approximately the same values as those ratios in a steady state of turbulence obtained by nonlinear simulations. In this work, linear and nonlinear simulation results are compared to show the validity of this assumption in the quasilinear modeling of energy exchange as well as particle and heat transport fluxes in the ITG turbulence. We here note that this work demonstrates only the feasibility of quasilinear modeling but does not present a saturation rule which is necessary for developing a quasilinear model. We also discuss physical mechanisms and the quasilinear modeling of turbulent energy exchange by wavenumber spectral analyses of entropy balance in microturbulence.9,10,23,24

The rest of this paper is organized as follows. In Secs. II A and II B, gyrokinetic equations and two balance equations related to perturbed entropy density and thermal energy are presented. The turbulent energy exchange is represented by wavenumber spectral functions in Sec. II C. In Sec. II D, the entropy balance in wavenumber space is investigated to consider the conditions for the quasilinear model to correctly predict the turbulent energy exchange and turbulent particle and heat transport fluxes. In Sec. III, results of the ITG turbulence by the GKV code,25 which uses a flux tube domain,28 are shown. Simulation settings are described in Sec. III A, and the turbulent energy exchange and transport fluxes obtained by the GKV simulation are shown as functions of T e / T i in Sec. III B. Comparisons between collisional and turbulent energy exchanges are made in Sec. III C. In addition, the result reported in Ref. 13, where the effect of turbulent energy exchange is shown to be negligible in DIII-D shot 128 913, is verified by the simulation using the same shot conditions. In Sec. III C, spectrum analyses of the turbulent energy exchange are performed to investigate its mechanisms and they are found to be directly connected with those of destabilizing the ITG modes. In Sec. III D, linear and nonlinear simulation results of the entropy balance in wavenumber space are compared to confirm the validity of the assumption of the quasilinear model regarding the quasilinear weights for the turbulent transport fluxes and energy exchange. Finally, conclusions and discussion are given in Sec. IV.

The plasma distribution function Fa of the position vector x, velocity vector v, and time t can be written as F a ( x , v , t ) = f a ( x , v , t ) + f ̂ a ( x , v , t ), where fa and f ̂ a represent the ensemble-averaged and fluctuation parts, respectively, and the subscript a denotes the particle species. The space-time scales of variations in the ensemble-averaged part fa are much larger than those of the fluctuation part f ̂ a so that the ensemble average can also be regarded as the local space-time average. In the gyrokinetic theory, perturbations such as f ̂ a are assumed to satisfy the gyrokinetic ordering
f ̂ a f a e a ϕ ̂ T a | B ̂ | | B | ω Ω a k k ρ t a L δ 1 ,
(1)
where ϕ ̂ and B ̂ are the perturbed electrostatic potential and the perturbed magnetic field, respectively. Here, Ω a = e a B / m a c , ρ t a = v t a / Ω a, ma, ea, Ta, B, c, and v t a T a / m a are the gyrofrequency, thermal gyroradius, mass, charge, temperature, background magnetic field, speed of light, and thermal velocity, respectively. The characteristic wavenumbers (in directions parallel and perpendicular to the background magnetic field), frequency, and equilibrium scale length are represented by k , k , ω, and L, respectively.
It is useful to express any perturbed function W ̂ in the WKB form,29 
W ̂ ( x , t ) = k W k exp ( i S k ( x , t ) ) ,
(2)
where S k is the eikonal whose gradient gives the perpendicular wavenumber vector S k = k . On the zeroth order in δ, the ensemble-averaged distribution function fa is assumed to be the local Maxwellian fMa, which is given in terms of the background density na and temperature Ta by
f M a = n a ( m a 2 π T a ) 3 / 2 exp ( m a v 2 2 T a ) ,
(3)
and the perturbed distribution function f a k can be written as follows:
f a k = e a ϕ k T a f M a + h a k e i k · ρ a ,
(4)
where ρ a = b × v / Ω a with b = B / B and B = | B |. Here, h a k represents the nonadiabatic distribution function which is calculated by the gyrokinetic equation,
( t + v b · + i k · ( v E + v d a ) ) h a k = e a T a f M a ( t + i k · ( v E + v * a ) ) ψ a k + c B k + k = k [ b · ( k × k ) ] ψ a k h a k + C a G K ,
(5)
where v = v · b, v E = c E × b / B, v d a = c b × ( μ B + m a v 2 b · b ) / e a B, v * a = c T a b × { ln n a + ( w 3 / 2 ) ln T a } / e a B , C a G K is the gyrokinetic collision term,10 and E = Φ is the background electric field. Here, the background E × B flow is assumed to be v E δ v t i, and the gradient scale length of v E is estimated as L. Therefore, the effect of electric field shear is neglected. The detail is discussed in  Appendix A. In Eq. (5), h a k is regarded as a function of time t and phase space variables (x, w = m a v 2 / 2 , μ = m a v 2 / 2 B), where v = | v | and v = | v v b |. The gyrophase-averaged perturbed potential function ψ a k is defined in terms of the perturbed electrostatic potential ϕ k and the perturbed vector potential A k as follows:
ψ a k = d ξ 2 π exp ( i k · ρ a ) ( ϕ k v c · A k ) = J 0 ( k ρ a ) ( ϕ k v c A k ) + J 1 ( k ρ a ) v c B k k ,
(6)
where J0 and J1 are the zeroth- and first-order Bessel functions, respectively, A k = b · A k , and B k = i b · ( k × A k ). The perturbed electric and magnetic fields are determined by Poisson's equation and Ampère's law in the gyrokinetic form as follows:
( k 2 + λ D 2 ) ϕ k = 4 π a e a d 3 v J 0 ( k ρ a ) h a k ,
(7)
k 2 A k = 4 π c a e a d 3 v J 0 ( k ρ a ) h a k v ,
(8)
k B k = 4 π c a e a d 3 v J 1 ( k ρ a ) h a k v ,
(9)
where λ D = 1 / a 4 π n a e a 2 / T a is the Debye length. The time evolution of nonadiabatic distribution function h a k can be obtained by solving Eq. (5) combined with Eqs. (7)–(9).
We now consider two balance equations related to the perturbed entropy density and the energy. These variables are associated with the particle and heat fluxes and the energy exchange between electrons and ions. The entropy density due to turbulence δ S a is defined by the difference between the macroscopic entropy density S M a d 3 v f a log f a and the ensemble averaged microscopic entropy density S m a ens d 3 v ( f a + f ̂ a ) log ( f a + f ̂ a ) ens as follows:
δ S a S M a S m a ens = k d 3 v | f a k | 2 2 f M a ens ,
(10)
where ens represents the ensemble average, and the terms of O ( δ 3 ) are neglected. Using Eq. (4), the integral in the angle bracket in Eq. (10) can be rewritten as follows:
d 3 v | f a k | 2 2 f M a = d 3 v | h a k | 2 2 f M a n a e a 2 2 T a 2 | ϕ k | 2 e a T a Re [ ϕ k * n a k ] ,
(11)
where * denotes the complex conjugate and n a k d 3 v f a k is the density perturbation with the perpendicular wavenumber vector k . Multiplying Eq. (5) by h a k * / f M a and taking the ensemble and flux-surface averages, we can derive the entropy balance equation as follows:
t ( δ S h a ) = Γ a turb L p a + q a turb T a L T a + Q a turb T a + D a ,
(12)
where δ S h a = k d 3 v | h a k | 2 / 2 f M a , and denotes a double average over the ensemble and the flux surface. The turbulent particle and heat fluxes Γ a turb , q a turb and the inverse gradient scale lengths L p a 1 , L T a 1 are defined by
[ Γ a turb , q a turb T a ] k [ Γ a k turb , q a k turb T a ] = k Re d 3 v [ 1 , w T a 5 2 ] × h a k * ( i c B ψ a k k × b ) · r ,
(13)
[ L p a 1 , L T a 1 ] = [ ln p a r e a T a Φ r , ln T a r ] ,
(14)
respectively, where the minor radius r is used as a label of flux surfaces of a toroidal plasma and Q a turb represents the turbulent heating of the particles of species a given by
Q a turb k Q a k turb = e a k Re d 3 v h a k * ψ a k t .
(15)
Taking the summation of Eq. (15) over the particle species and using Eqs. (7)–(9), we obtain
a Q a turb = 1 8 π t k ( k 2 + λ D 2 ) | ϕ k | 2 | B k | 2 ,
(16)
where B k = i k × A k . Using the transport timescale ordering ( / t ) / = O ( δ 2 ), we obtain a Q a turb = 0 to the lowest order in δ. Therefore, we can interpret Q a turb as the energy exchange between different plasma species. The last term in Eq. (12) Da represents collisional dissipation written as
D a k D a k = k Re d 3 v h a k * f M a C a G K .
(17)
The left-hand side of Eq. (12) vanishes in the steady state of turbulence where the entropy production due to turbulent transport and heating balance with the collisional dissipation.
The ensemble and flux-surface averaged energy balance equation is given by9 
t ( 3 2 p a ) + 1 V r [ V ( q a + 5 2 T a Γ a ) ] = ( Γ a cl + Γ a ncl ) n a p a r + u a · ( · π a ) + Q a coll e a Γ a turb Φ r + Q a turb ,
(18)
where represents the flux-surface average, and pa and π a are the pressure and viscosity tensor, respectively. Here, V is the derivative of the volume V inside a magnetic surface with respect to the minor radius r. The radial particle and heat fluxes denoted by Γ a and qa, respectively, are written as follows:
Γ a = Γ a cl + Γ a ncl + Γ a turb ,
(19)
q a = q a cl + q a ncl + q a turb ,
(20)
where the superscripts cl , nc, and turb denote classical, neoclassical, and turbulence parts, respectively. Here, Q a coll is the flux-surface average of the collisional heat generation rate. In the case of a plasma consisting of electrons and single-species ions, we have
Q i coll = 3 m e m i n e ν e ( T e T i ) ,
(21)
ν e = 4 2 π n i e e 2 e i 2 ln Λ 3 m e 1 2 T e 3 2 ,
(22)
and
Q e coll = R e · ( u e u i ) a Q i coll ,
(23)
where νe is the electron-ion collision frequency, ln Λ is the Coulomb logarithm, Re is the collisional friction force for electrons, and ( u e u i ) is the difference between the electron and ion flow velocities. We see that the collisional energy transfer Q i coll from electrons to ions is proportional to the product of the collision frequency and the temperature difference between electrons and ions.

 Appendix A shows that in the present model, the effect of the radial electric field Er on turbulent energy exchange appears only in the Doppler shift. The sum of the last two terms on the right-hand side of Eq. (18) is essentially equivalent to turbulent energy exchange Eq. (15) in Er = 0. Therefore, we henceforth discuss turbulent energy exchange with Er = 0.

Electrons and ions exchange energy via collisions and turbulence. The collisional energy exchange decreases for low collision frequency. It always transfers energy from the hotter species to the colder one and vanishes when the two species have the same temperature. As shown later from the gyrokinetic analysis and simulation, the turbulent energy exchange has quite different properties from the collisional one.

We investigate the physical mechanism of the turbulent energy exchange by using Eq. (15). In the steady state of turbulence, the time derivative acting on the perturbed potential Eq. (5) can be transferred to that on the nonadiabatic perturbed distribution function with the help of the Leibniz rule. As shown in  Appendix B, using the gyrokinetic equation, Eq. (5), we can rewrite Eq. (15) as follows:
Q a turb = e a k Re d 3 v ψ a k * h a k t = k ( Q a k turb + Q a B k turb + Q a ψ k turb + Q a C k turb ) ,
(24)
with
Q a k turb = Re d 3 v G a k · j a k * ,
(25)
Q a B k turb = Re d 3 v G a k · j a B k * ,
(26)
Q a ψ k turb = Re d 3 v G a k · j a ψ k * ,
(27)
Q a C k turb = Re d 3 v C a G K ψ a k * ,
(28)
where Er = 0 is used. The perturbed fields G a k and the perturbed currents j a k at the gyrocenter position are defined by
G a k = ψ a k i k ψ a k = G a k + G a k ,
(29)
j a k = j a k + j a k ,
(30)
j a k = e a h a k v b ,
(31)
j a k = j a B k + j a ψ k ,
(32)
j a B k = e a h a k v d a ,
(33)
j a ψ k = i c e a B k + k = k ( b × k ) ψ a k h a k .
(34)
Here, j a ψ k is derived from the nonlinear term of Eq. (5) and represents the current induced by the perturbed potential. As shown in Eq. (24), the turbulent energy exchange is caused by the product of the perturbed field and current due to the nonadiabatic distribution. The field and current can be decomposed into components in directions parallel and perpendicular to the background magnetic field. The perpendicular current component can also be classified as the two parts: the one is produced by the B-curvature drift in the toroidal magnetic field and the other by the drift due to the turbulent potential field.

Furthermore, the effect of collisions is given by the last term at the right-hand side of Eq. (24), and thus, the turbulent energy exchange is represented by the sum of the four parts.

As mentioned earlier, in the steady state of turbulence, the turbulent energy exchange does not cause a net increase or decrease in energy. This property is also valid for each wavenumber as shown by
a Q a k turb = 0 ,
(35)
where Q a k turb = Q a k turb + Q a B k turb + Q a ψ k turb + Q a C k turb. Furthermore, we have
k Q a ψ k turb = 0 ,
(36)
which can be illustrated by using the definition of Q a ψ k turb given in Eqs. (27), (29), and (34). Equation (36) implies that Q a ψ k turb does not contribute to the net heating or cooling of particles of the species a but it represents the energy transfer in the wavenumber space through nonlinear interactions between different modes. Thus, Q a ψ k turb influences the profile of the total wavenumber spectrum Q a k turb. Because of Eq. (36), the total turbulent energy exchange Q a turb is given by taking the sum of the three components Q a k turb , Q a B k turb, and Q a C k turb over the whole wavenumber space. In particular, Q a k turb , Q a B k turb are the contributions to Joule heating (cooling) via currents parallel and perpendicular to the background magnetic field.
It is noted here that, using Eq. (5) and Eqs. (7)–(9), we can obtain
a e a k Re d 3 v ψ a k * h a k t = t k [ 1 8 π | E k | 2 | B k | 2 + a n a e a 2 | ϕ k | 2 2 T a ] ,
(37)
where the steady state of turbulence is not assumed. Under an electrostatic approximation, magnetic fluctuations are neglected, and the right-hand side of Eq. (37) is regarded as the decrease rate of the energy associated with the electrostatic electric field. Then, it is understood from Eqs. (24) and (37) that the turbulent heating of particles represented by a k ( Q a k turb + Q a B k turb + Q a C k turb ) is brought about by consuming the electrostatic electric field energy.
In this section, the predictability of the quasilinear model for the turbulent energy exchange is discussed with the entropy balance in the linear and nonlinear states. When using the solutions h a k and ψ a k of the linearized version of Eq. (5) combined with Eqs. (7)–(9) to evaluate Q a turb defined in Eq. (15), we obtain
Q a turb = e a k Re d 3 v ( i ω r k + γ k ) h a k * ψ a k ,
(38)
where ω r k and γ k are the real and imaginary parts of the complex-valued linear eigenfrequency ω k ω r k + i γ k for the wavenumber vector k . Although a Q a turb = 0 holds in the steady state of turbulence, it does not for the linear solutions. The presence of the finite growth rate γ k causes a Q a turb 0. Then, we drop γ k from Eq. (38) and define
Y a = e a k Re d 3 v ( i ω r k ) h a k * ψ a k ,
(39)
to estimate the turbulent energy exchange from the linear solutions because a Y a = 0 is satisfied. We now note that Ya can be rewritten as follows:
Y a k Y a k = e a 2 k Re d 3 v ( h a k * ψ a k t h a k * t ψ a k ) ,
(40)
which is the same as the turbulent energy exchange introduced by Candy.13 This expression of Ya in Eq. (40) can be used for both linear and nonlinear cases and rigorously satisfies a Y a = 0 even in non-steady states. The relation between Q a turb and Ya is given from Eqs. (15) and (40) as follows:
Q a turb Y a = t { e a 2 k Re d 3 v h a k * ψ a k } ,
(41)
from which we easily see that Q a turb = Y a in the turbulent steady state.
The entropy balance equation for each wavenumber is derived from Eq. (5) as follows:
t d 3 v | h a k | 2 2 f M a = Γ a k turb L p a + q a k turb T a L T a + Q a k turb T a + D a k + N a k ,
(42)
where the subscript k denotes the contribution from each perpendicular wavenumber vector and N a k is defined by
N a k = d 3 v c B f M a k + k = k [ b · ( k × k ) ] ψ a k h a k h a k * .
(43)
Here, N a k represents the entropy which the mode with the wavenumber vector k gains through nonlinear interaction with other modes, and it satisfies
k N a k = 0 ,
(44)
which implies that the nonlinear interaction produces no net entropy.
Substituting Eq. (41) into Eq. (42) and using the following relation:
d 3 v [ | h a k | 2 2 f M a e a 2 T a Re [ h a k * ψ a k ] ] = d 3 v | f a k | 2 2 f M a + e a 2 T a ( n a k * ϕ k + 1 c n a u a k * · A k ) ,
(45)
we obtain
t d 3 v | f a k | 2 2 f M a + e a 2 T a ( n a k * ϕ k + 1 c n a u a k * · A k ) N a k D a k = Γ a k turb L p a + q a k turb T a L T a + Y a k T a ,
(46)
where u a k is the perturbed flow velocity with the perpendicular wavenumber vector k defined by n a u a k d 3 v f a k v. In the double angle brackets on the left-hand side of Eq. (46), we find the entropy due to the perturbed particle distribution function f a k as well as the inverse temperature multiplied by the electrostatic potential energy and the interaction between the magnetic potential and the current of particle species a.

In the quasilinear model, the transport fluxes divided by the squared potential, Γ a k turb / | ϕ k | 2 and q a k turb / | ϕ k | 2 , in the turbulent steady state are approximated by the corresponding values obtained from the linear analysis. We here include the turbulent energy exchange term into the quasilinear model and estimate Y a k / | ϕ k | 2 from the linear calculation. When this model is valid, the values of the terms on the right-hand side of Eq. (46) divided by | ϕ k | 2 should not change whether linear or nonlinear simulations are performed to evaluate them. On the other hand, when divided by | ϕ k | 2 , the time-derivative term and the nonlinear entropy transfer term on the left-hand side of Eq. (46), take different values between the linear and nonlinear cases. The time-derivative term is dominant and the nonlinear entropy transfer vanishes in the linear case, while the former is negligible and the latter is dominant in the nonlinear case. Here, the ratio of the collisional dissipation to the squared potential, D a k / | ϕ k | 2 , is assumed to take the same value. Then, it is concluded from the above-mentioned discussion of Eq. (46) based on the quasilinear model that the ratio of the time-derivative term divided by | ϕ k | 2 in the linear case should be the same as the value of N a k / | ϕ k | 2 in the nonlinear case in order to keep the balance between the left- and right-hand sides of Eq. (46).

We here point out the resemblance of Eq. (46) to the well-known Landau equation of the weakly nonlinear theory for fluid dynamic systems30,31
d | A | 2 d t = 2 γ | A | 2 α | A | 4 ,
(47)
where | A | and γ are the amplitude of the dominant mode and its linear growth rate, respectively, and α | A | 4, with a positive constant α, represents the nonlinear effect which causes the saturation of the mode. In the linearly growing phase of the mode, the time derivative term on the left-hand side of Eq. (47) equals the first term on the right-hand side while in the steady state and the second nonlinear saturation term on the right-hand side balances with the first term. This exchange of the roles between the time-derivative and nonlinear terms in Eq. (47) is common to the process described about the quasilinear model using Eq. (46).
The summation of Eq. (46) over species can be written as follows:
t a T a d 3 v | f a k | 2 2 f M a + e a 2 ( n a k * ϕ k + 1 c n a u a k * · A k ) = t [ a T a d 3 v | f a k | 2 2 f M a + 1 8 π | E k | 2 + | B k | 2 ] = a [ T a Γ a k turb L p a + q a k turb L T a + T a N a k + T a D a k ] ,
(48)
where Eqs. (7)–(9) are used. We can see that terms inside the brackets on the second line of Eq. (48) represent the fluctuation entropy and the electromagnetic energy which are all positive, so that their time-derivative terms are positive in the unstable wavenumber region in the linear phase of the time evolution of the fluctuation. The magnetic fluctuations can be neglected in a low beta plasma, and the nonadiabatic part of the electron distribution function is small for ITG instability. Then, it is considered from Eq. (45) that, in the case of the ITG mode, the contribution from electrons to the time-derivative part in Eq. (48) is negligible and ions contribute dominantly. Therefore, the time-derivative part on the left-hand side of Eq. (46) is expected to be positive for ions in the linear phase, which implies that the nonlinear entropy transfer term for ions in the unstable wavenumber region should be negative, N i k < 0, in the nonlinear steady state according to the quasilinear argument given earlier. On the other hand, the summation of N i k over the linearly stable wavenumber region should be positive because k N i k = 0. Thus, the entropy due to the fluctuations is transferred from the unstable wavenumber region to the stable one.

The discussion of quasilinear weights with entropy balance is expected to apply not only to tokamaks but also to more generic 3D geometries. In Sec. III D, the relative magnitude of each term in the entropy balance equation, Eq. (42), is evaluated by the linear and nonlinear gyrokinetic simulations, and the speculations about the quasilinear model described above are examined. It should be pointed out that in this work, we discuss only quasilinear weights and do not investigate the saturation rule required to predict fluxes.

In this research, microturbulence simulations are performed by the GKV code,25 which solves the gyrokinetic equation for the perturbed distribution function based on the Eulerian scheme in ( k x , k y , z , v , μ ) space. The nonlinear term is evaluated in the real space, and is transformed back to the wavenumber space by means of 2D fast Fourier transform algorithm and the 2/3 rule in (kx, ky). It employs the local flux-tube domain where the background densities, temperatures, and their gradients are fixed. While we do not deal with 3D geometries here, microturbulence in helical systems has also been studied by GKV.26,27

The flux tube coordinates for a low-β, large aspect ratio, axisymmetric torus with concentric circular cross sections, x = r r 0 , y = r 0 / q 0 ( q ( r ) θ ζ ) , and z = θ are used in this work, where r, θ, and ζ are minor radius, poloidal angle, and toroidal angle, respectively. The subscript 0 denotes parameters at the center of flux tube. The perpendicular wavenumber is given by k = ( k x + s ̂ z k y ) e r + k y e θ. Here, kx and ky are wavenumbers in the directions of x and y, respectively, while e r and e θ are unit vectors parallel to r and θ, respectively.

We here focus on the ion temperature gradient (ITG) mode turbulence in tokamak plasmas so that the electron temperature gradient is set to zero, R 0 / L T e = 0, where R0 represent the major radius. Plasma and field parameters used in the simulations are shown in Table I. Most of them are the same values as in the Cyclone DIII-D base case.4 The ion beta value is set to β i = 1 × 10 4 for which the electrostatic approximation is valid. In simulations performed here, B k is neglected, although A k is retained in order to avoid numerical difficulty due to very rapid electrostatic waves called the ωH mode.32 The Lenard–Bernstein collision operator33 is used here because it takes less computation time than more rigorous collision models. However, we expect that the collision model does not influence results from the present study where the normalized collision frequency is set to ν i i * R 0 q ν i i / ( ϵ 3 / 2 v t i ) = 0.068, which is much smaller than the growth rates of the ITG modes in the present study. Since collisional energy exchange is proportional to the temperature difference between electrons and ions, the temperature ratio is set in a range from T e / T i = 0.80 to 1.5 in order to compare turbulent and collisional energy exchanges.

TABLE I.

Plasma and field parameters.

Plasma and field parameters Value
Normalized ion temperature gradient  R 0 / L T i  6.92 
Normalized electron temperature gradient  R 0 / L T e 
Normalized density gradient  R 0 / L n a  2.22 
Mass ratio  m e / m i  5.446 × 10 4 
Charge ratio  e e / e i  −1 
Ion beta value  β i = 4 π n i T i / B 2  1 × 10 4 
Normalized collision frequency  ν i i * R 0 q 0 ν i i / ( ϵ r 3 / 2 v t i )  0.068 
Temperature ratio  T e / T i  0.8 1.5 
Inverse aspect ratio  ε r  0.18 
Safety factor  q0  1.4 
Magnetic shear  s ̂ = ( r / q ) ( d q / d r )  0.78 
Plasma and field parameters Value
Normalized ion temperature gradient  R 0 / L T i  6.92 
Normalized electron temperature gradient  R 0 / L T e 
Normalized density gradient  R 0 / L n a  2.22 
Mass ratio  m e / m i  5.446 × 10 4 
Charge ratio  e e / e i  −1 
Ion beta value  β i = 4 π n i T i / B 2  1 × 10 4 
Normalized collision frequency  ν i i * R 0 q 0 ν i i / ( ϵ r 3 / 2 v t i )  0.068 
Temperature ratio  T e / T i  0.8 1.5 
Inverse aspect ratio  ε r  0.18 
Safety factor  q0  1.4 
Magnetic shear  s ̂ = ( r / q ) ( d q / d r )  0.78 

The resolution settings are shown in Table II. High resolution for the parallel coordinate z is used in Figs. 1(b) and 4 to suppress an error of entropy balances by a hyper diffusion term.34 There is no significant difference in turbulent fluxes and energy exchange when using either low or high resolutions for z. The output data shown in this paper are normalized by following Ref. 25 except for Fig. 2.

TABLE II.

Resolution settings.

Domain sizes and resolved perpendicular wavenumbers 
64.10 ρ t i x 64.10 ρ t i 
62.83 ρ t i y 62.83 ρ t i 
π z π 
4 v t a v 4 v t a 
0 μ B 0 / T a m a v 2 / 2 T a 8 
4.70 ρ t i 1 k x 4.70 ρ t i 1 , ( k x , min = 0.049 ρ t i 1 ) 
1.55 ρ t i 1 k y 1.55 ρ t i 1 , ( k y , min = 0.050 ρ t i 1 ) 
Grid points in ( x , y , z , v , μ ) 
288 × 96 × 64 × 64 × 32 
288 × 96 × 256 × 64 × 32 
Domain sizes and resolved perpendicular wavenumbers 
64.10 ρ t i x 64.10 ρ t i 
62.83 ρ t i y 62.83 ρ t i 
π z π 
4 v t a v 4 v t a 
0 μ B 0 / T a m a v 2 / 2 T a 8 
4.70 ρ t i 1 k x 4.70 ρ t i 1 , ( k x , min = 0.049 ρ t i 1 ) 
1.55 ρ t i 1 k y 1.55 ρ t i 1 , ( k y , min = 0.050 ρ t i 1 ) 
Grid points in ( x , y , z , v , μ ) 
288 × 96 × 64 × 64 × 32 
288 × 96 × 256 × 64 × 32 
FIG. 1.

(a) The linear growth rate γ and frequency ωr at ( k x ρ t i , k y ρ t i ) = ( 0 , 0.30 ), the turbulent ion heating Q i turb, the turbulent particle flux Γ e turb = Γ i turb, and the turbulent heat fluxes q s turb ( a = e , i ) plotted as functions of the temperature ratio T e / T i. (b) Comparison of all terms in the entropy balance equation, Eq. (12), in the saturated state of the ITG turbulence for T e / T i = 1.0. All terms are normalized by the entropy production q i turb / T i L T i due to the turbulent ion heat flux.

FIG. 1.

(a) The linear growth rate γ and frequency ωr at ( k x ρ t i , k y ρ t i ) = ( 0 , 0.30 ), the turbulent ion heating Q i turb, the turbulent particle flux Γ e turb = Γ i turb, and the turbulent heat fluxes q s turb ( a = e , i ) plotted as functions of the temperature ratio T e / T i. (b) Comparison of all terms in the entropy balance equation, Eq. (12), in the saturated state of the ITG turbulence for T e / T i = 1.0. All terms are normalized by the entropy production q i turb / T i L T i due to the turbulent ion heat flux.

Close modal
FIG. 2.

Comparison of energy exchanges due to Coulomb collisions (blue squares) and turbulence (orange circles). The left figure (a) shows the case of δ ρ t i / R 0 = 9.6 × 10 4, n e = 2.0 × 10 19 ( m 3 ), and T i = 0.9 ( keV ), and the right one (b) shows the case of δ = 1.7 × 10 3, n e = 2.25 × 10 20 ( m 3 ), and T i = 3.0 ( keV ). Red triangles and black stars indicate turbulent energy exchanges obtained from simulations for T e / T i = 1.2 and R 0 / L n a = 3.0 using different temperature gradients given by ( R 0 / L T i , R 0 / L T e ) = ( 5.0 , 0 ) and ( R 0 / L T i , R 0 / L T e ) = ( 5.0 , 6.5 ), respectively.

FIG. 2.

Comparison of energy exchanges due to Coulomb collisions (blue squares) and turbulence (orange circles). The left figure (a) shows the case of δ ρ t i / R 0 = 9.6 × 10 4, n e = 2.0 × 10 19 ( m 3 ), and T i = 0.9 ( keV ), and the right one (b) shows the case of δ = 1.7 × 10 3, n e = 2.25 × 10 20 ( m 3 ), and T i = 3.0 ( keV ). Red triangles and black stars indicate turbulent energy exchanges obtained from simulations for T e / T i = 1.2 and R 0 / L n a = 3.0 using different temperature gradients given by ( R 0 / L T i , R 0 / L T e ) = ( 5.0 , 0 ) and ( R 0 / L T i , R 0 / L T e ) = ( 5.0 , 6.5 ), respectively.

Close modal

Here, we perform linear and nonlinear simulations where the temperature ratio is varied as a parameter. Figure 1(a) shows the linear growth rate for ( k x ρ t i , k y ρ t i ) = ( 0.00 , 0.30 ), particle and heat fluxes, and energy exchange as functions of the temperature ratio T e / T i. The results except for the linear growth rate are obtained by taking a time average in a steady state of turbulence. It is seen that, as T e / T i increases, the linear growth rate, the absolute values of particle and heat fluxes and turbulent energy exchange increase.

The ratio of each entropy balance term to the entropy production term caused by the ion heat flux q i turb / T i L T i is shown in Fig. 1(b). The numerical error in the entropy balance defined by the difference between the left- and right-hand sides of Eq. (12) is within 6% of q i turb / T i L T i for ions and 3% for electrons. In the case of ions, the particle and heat fluxes generate entropy while the turbulent energy exchange and collisions reduce entropy, thus maintaining balance. For electrons, on the other hand, the electron heat flux does not appear as an entropy-producing term because of no electron temperature gradient, although the turbulent energy exchange plays a major role in the entropy production in addition to the electron particle flux. The generated entropy is dissipated by collisions, and the entropy of turbulent fluctuations in the electron distribution function is kept in a steady state. The same result as described above is also reported in Ref. 24.

In ITG turbulence, the turbulence entropy of ions is generated primarily by the product of the ion heat flux and the ion temperature gradient and second by that of the particle flux and the ion pressure gradient. The ion turbulence entropy is lost mainly through collisions, although it is partially transferred to the electron turbulence entropy by the turbulent energy exchange. Electrons increase their turbulence entropy by the energy transfer from ions and the particle flux, while they reduce it through collisions. Thus, the total turbulence entropy balance in the steady state of the ITG turbulence is maintained by the turbulent energy transfer from ions to electrons, which carries the excess of the stronger ion entropy production to the weaker electron portion.

Here, the results of the turbulent energy exchange shown in Fig. 1(a) are used for comparison with the collisional energy exchange calculated by Eq. (21). If we treat the Coulomb logarithm as a constant ( ln Λ = 15.5) and vary the density and temperature with keeping n / T 2 fixed, the normalized collision frequency used as an input parameter of the GKV code does not change. Therefore, the results in Fig. 1(a) can be used for the two density and temperature conditions with the same value of n / T 2 shown in Figs. 2(a) and 2(b). The results for the density–temperature condition at r = a / 2 at the DIII-D128913 shot [ 0.9 ( keV ) , 2.0 × 10 19 ( m 3 )] and for the [ 3.0 ( keV ) , 2.25 × 10 20 ( m 3 )] are plotted in Figs. 2(a) and 2(b), respectively.

The collisional and turbulent energy transfers from electrons to ions, Q i coll and Q i turb, are shown as functions of the temperature ratio T e / T i in Fig. 2, where we can identify a difference between the directions of collisional and turbulent energy transfers. In Coulomb collisions, energy is transferred always from higher temperature particles to lower temperature ones. Thus, when T e / T i < 1 ( > 1 ) is less (more) than unity, Q i coll is negative (positive). On the other hand, the ITG turbulence always transfers energy from ions to electrons regardless of the value of T e / T i. We can see that Q i turb < 0 even in the equithermal condition T e / T i = 1, where Q i coll = 0 and that Q i coll and Q i turb take opposite signs to each other when T e / T i > 1.

It can also be confirmed by comparing with Figs. 2(a) and 2(b) that as the temperatures for electrons and ions are increased with their ratio T e / T i fixed, the energy exchange by turbulence becomes more dominant than that by collision. Note that, based on the gyrokinetic ordering, the turbulent ion heating is written as Q i turb = ( ρ t i / R 0 ) 2 ( v t i / R 0 ) ( n T i ) Q *, where Q * is a dimensionless function of dimensionless variables T e / T i , R 0 / L n , R 0 / L T i , R 0 / L T e , q 0 , s ̂ , β , R 0 ν i / v t i , which are used as input parameters of the local flux-tube simulation. This expression is combined with Eq. (21) to obtain
Q i turb Q i coll = m i / m e 3 ( T e / T i 1 ) v t i R 0 ν e ( ρ t i R 0 ) 2 × Q * ( T e T i , R 0 L n , R 0 L T i , R 0 L T e , q 0 , s ̂ , β , R 0 ν i / v t i , ) T i B 2 R 0 3 ( R 0 ν i / v t i ) Q * ( T e / T i 1 ) T i 3 n B 2 R 0 3 Q * ( T e / T i 1 ) .
(49)
When the normalized values T e / T i , R 0 / L n , R 0 / L T i , R 0 / L T e , q 0 , s ̂ , β , R 0 ν i / v t i , are regarded as constants, Q * is also a constant and the ratio Q i turb / Q i coll between the turbulent and collisional energy exchanges is proportional to the temperature. This temperature dependence of Q i turb / Q i coll is seen by comparing the results in Figs. 2(a) and 2(b). If taking account of the temperature dependence of R 0 ν i / v t i and assuming Q * to depends weakly on R 0 ν i / v t i, Q i turb / Q i coll is proportional to the cubic of the temperature. In the case of high plasma temperatures with T e / T i > 1, the net energy transfer from lower-temperature ions to higher-temperature electrons can occur, contrary to conventional thought.

It is reported in Ref. 13 that the turbulent energy exchange has a negligible effect on the simulation for predicting the temperature profile in the case of DIII-D128913. Here, we compare that case with the Cyclone D-III D base case (CBC) used in our simulations. The normalized density and temperature gradients in DIII-D128913 are estimated as R 0 / L n a = 3.0 , R 0 / L T i = 5.0, which is smaller than R 0 / L T i = 6.92 at r / a = 0.5 in CBC and R 0 / L T e = 6.5 at the minor radius r / a = 0.5 from Refs. 13 and 35. In Fig. 2, turbulent energy exchanges obtained from simulations for T e / T i = 1.2 and R 0 / L n a = 3.0 using ( R 0 / L T i , R 0 / L T e ) = ( 5.0 , 0 ) and ( R 0 / L T i , R 0 / L T e ) = ( 5.0 , 6.5 ) are plotted by red triangles and black stars, respectively. Interestingly, these plots indicate that the dependence of the turbulent energy exchange on R 0 / L T e is weak, while the turbulent transport fluxes are found to increase significantly with increasing R 0 / L T e. It is speculated that, even though both ion and electron turbulence entropy production due to transport fluxes increase, their difference stays the same so that the entropy balance for each species is maintained with the turbulent energy exchange unaltered. The black star in the left figure for n e = 2.0 × 10 19 ( m 3 ) and T i = 0.9 ( keV ) corresponds to the conditions at r / a = 0.5 in DIII-D128913, which shows that the magnitude of the turbulent energy exchange becomes significantly smaller than that of the collisional energy exchange. Thus, in this case, the turbulent energy exchange has only a small influence, which is consistent with the result in Ref. 13.

On the other hand, if the ion temperature gradient is so large that the turbulence has a dominant effect on the energy exchange, a net energy transfer can occur from ions to electrons even for Te > Ti. In particular, since the energy exchange due to the ITG turbulence acts to prevent the energy transfer from electrons to ions, the ITG turbulence is undesirable for fusion reactors in that it interferes with the ion heating by the alpha-heated electrons, as well as degrading the ion energy confinement through enhancing the energy transport toward the outside of the device.

Next, we examine each component of the wavenumber spectrum of the turbulent energy exchange shown in Eq. (24). Figure 3(a) shows the wavenumber spectra of the linear growth rate and frequency. The turbulent energy transfer terms in Eqs. (24)–(28) in the case of T e / T i = 1.0 are shown for electrons and ions in Figs. 3(b) and 3(c), respectively. First, the collisional components Q e C k turb and Q i C k turb have little effect on the turbulent energy exchange. In addition, both electrons and ions show positive values of the parallel-heating components Q a k turb > 0 ( a = e , i ), and negative values of the components Q a B k turb < 0 ( a = e , i ) caused by the product of the perpendicular field and the B -curvature drift velocity in the wavenumber region where the ITG mode is linearly unstable. At the bottom of Fig. 3, the roles of turbulent energy transfer terms in increasing or decreasing the thermal electron and ion energies, and the turbulent electrostatic electric field energy are schematically shown, based on Eqs. (18), (24), and (37) with the magnetic fluctuations neglected. We can see that the perpendicular ion cooling represented by Q i B turb < 0 is dominant and overcomes the parallel ion heating Q i turb > 0, which leads to the net turbulent ion cooling shown by Q i turb < 0. On the other hand, the parallel heating Q e turb > 0 is dominant for electrons. Thus, the perpendicular ion cooling and the parallel electron heating are found to be the main mechanisms in the turbulent energy transfer from ions to electrons in a steady state of ITG turbulence.

FIG. 3.

The wavenumber spectra of the linear growth rate and frequency (a) and the turbulent energy transfer terms in Eqs. (24)–(28) for electrons (b) and ions (c) in the case of T e / T i = 1.0. The spectra are given as functions of k y ρ t i obtained by summing over kx. The peaks and valleys of the turbulent energy transfer terms are found at a wavenumber lower than that of the linearly most unstable mode. The directions and magnitudes of the turbulent energy transfer terms are represented by arrows in the bottom figure which schematically shows the role of turbulent energy transfer terms based on Eqs. (18), (24), and (37) with magnetic fluctuations neglected. The electron heating due to the parallel field denoted by Q e turb > 0 and the ion cooling due to the B-curvature drift denoted by Q i B turb < 0 are dominant mechanisms in the turbulent energy exchange between electrons and ions in ITG turbulence.

FIG. 3.

The wavenumber spectra of the linear growth rate and frequency (a) and the turbulent energy transfer terms in Eqs. (24)–(28) for electrons (b) and ions (c) in the case of T e / T i = 1.0. The spectra are given as functions of k y ρ t i obtained by summing over kx. The peaks and valleys of the turbulent energy transfer terms are found at a wavenumber lower than that of the linearly most unstable mode. The directions and magnitudes of the turbulent energy transfer terms are represented by arrows in the bottom figure which schematically shows the role of turbulent energy transfer terms based on Eqs. (18), (24), and (37) with magnetic fluctuations neglected. The electron heating due to the parallel field denoted by Q e turb > 0 and the ion cooling due to the B-curvature drift denoted by Q i B turb < 0 are dominant mechanisms in the turbulent energy exchange between electrons and ions in ITG turbulence.

Close modal
This perpendicular ion cooling is connected with the mechanism of ITG instability. In low beta plasmas, the B-curvature drift can be expressed as v d a = b × ( v 2 + μ B / m a ) B / Ω a B. Then, under the electrostatic approximation, Q i B k turb in Eq. (26) is rewritten as follows:
Q i B k turb = Re 2 P i k * ( c B E k × b ) · ln B ,
(50)
where E k i k ϕ ̂ is the electric field and P i k d 3 v ( m i v 2 +  μ B ) f i k / 2 roughly represents ion pressure perturbation. The ITG mode is destabilized at the outside of the torus (or the bad curvature region) where the ion pressure perturbation is amplified by the outward E × B flow from the inner hot plasma region. Therefore, the phases of P i k and ( c E k × b / B ) · ln B become opposite to each other, and accordingly, Q i B turb < 0 is expected for the ITG instability from Eq. (49). We also see from Eq. (49) that Q i B turb represents the outward energy flow in the bad curvature. Thus, Q i B turb < 0 means the outward energy transport Q i B turb > 0 due to the ITG turbulence.

In Fig. 3, the nonlinear interaction between different wavenumbers is shown by the components Q e ψ k turb and Q i ψ k turb. They show that the energy in the unstable wavenumber region is carried to the zonal flow mode with ky = 0. This implies that the zonal flows play a significant role in the nonlinear saturation of the ITG turbulence. It is also confirmed that k Q a ψ k turb = 0 ( a = e , i ) and the nonlinear interaction causes no net energy production.

In Secs. III B and III C, we investigated the characteristics and the physical mechanism of the energy exchange due to the ITG turbulence. In particular, the comparison with the collisional energy exchange has clarified that the turbulent energy exchange can play a dominant role in the energy exchange between electrons and ions in low collision or high temperature plasmas. Therefore, it is necessary to take account of the effects of the turbulent energy exchange along with those of the particle and heat fluxes for reliable predictions of the global density and temperature profiles in future fusion reactors. In this subsection, the nonlinear simulation results of the ratio of the turbulent energy exchange to the squared amplitude of the electrostatic potential, are compared with the linear simulation results in order to examine the validity of the quasilinear model for predicting turbulent energy exchange. Equation (40) is used here for evaluating the turbulent energy exchange from the linear simulations.

Figures 4(a) and 4(b) show the wavenumber spectra of all terms in the entropy balance equations for electrons and ions, respectively. They are evaluated in the steady state of turbulence obtained by the nonlinear simulation for T e / T i = 1.0. The turbulent ion heat flux under the ion temperature gradient makes the largest contribution to the entropy production in the unstable wavenumber region where the particle flux under the pressure gradient also produces the entropy for both electrons and ions. On the other hand, no contribution is made by the electron heat flux because the electron temperature gradient is set to zero in the present simulation condition. The entropy produced by the unstable modes in the ITG turbulence is transferred to the zonal flow modes around ky = 0 and to the high-wavenumber modes, while the collisional entropy dissipation represented by D e < 0 and D i < 0 occurs in a wide wavenumber range, which maintains the detailed entropy balance in each wavenumber.

FIG. 4.

The wavenumber spectra of terms in the entropy balance equation in Eq. (42) for electrons (a) and ions (b) in the case of T e / T i = 1.0. The spectra are given as functions of k y ρ t i obtained by summing over kx. They are evaluated in the steady state of turbulence obtained by nonlinear simulation. The ratios of the turbulent particle and heat transport fluxes and the turbulent energy exchange to the squared amplitude of the electrostatic potential obtained by linear and nonlinear simulations are shown for electrons (c) and ions (d). The nonlinear entropy transfer terms N e k y and N i k y are negative in the wavenumber regions colored in sky blue. Dashed lines in (c) and (d) represent the ratios obtained by linearly unstable modes with k x = 0.

FIG. 4.

The wavenumber spectra of terms in the entropy balance equation in Eq. (42) for electrons (a) and ions (b) in the case of T e / T i = 1.0. The spectra are given as functions of k y ρ t i obtained by summing over kx. They are evaluated in the steady state of turbulence obtained by nonlinear simulation. The ratios of the turbulent particle and heat transport fluxes and the turbulent energy exchange to the squared amplitude of the electrostatic potential obtained by linear and nonlinear simulations are shown for electrons (c) and ions (d). The nonlinear entropy transfer terms N e k y and N i k y are negative in the wavenumber regions colored in sky blue. Dashed lines in (c) and (d) represent the ratios obtained by linearly unstable modes with k x = 0.

Close modal
The ratios of the turbulent particle and heat transport fluxes and the turbulent energy exchange to the squared amplitude of the electrostatic potential are plotted for electrons and ions in Figs. 4(c) and 4(d), respectively. Here, it should be recalled that these ratios obtained by linear simulations are called quasilinear weights. On the other hand, those obtained by nonlinear simulations are called nonlinear weights here. Dashed and solid curves represent the quasilinear and nonlinear weights, respectively. Solid curves labeled ( W Γ a , N , W q a , N , W Y a , N ) in Figs. 4(c) and 4(d) show the nonlinear weights as functions of k y ρ t i calculated by
W Γ a , N ( k y ) = T i 2 e 2 n e v t i k x Γ ̃ a k , N k x | ϕ k , N | 2 ,
(51)
W q a , N ( k y ) = T i e 2 n e v t i k x q ̃ a k , N k x | ϕ k , N | 2 ,
(52)
W Y a , N ( k y ) = T i R 0 e 2 n e v t i k x Y ̃ a k , N k x | ϕ k , N | 2 ,
(53)
where Γ ̃ a k , q ̃ a k , and Y ̃ a k are defined by the real parts of integrals inside the double average over the ensemble and the flux surface of Eqs. (13) and (40),
Γ ̃ a k Re [ d 3 v h a k * ( i c B ψ a k k × b ) · r ] ,
(54)
q ̃ a k Re [ d 3 v ( w 5 T a 2 ) h a k * ( i c B ψ a k k × b ) · r ] ,
(55)
Y ̃ a k Re [ e a 2 d 3 v ( h a k * ψ a k t h a k * t ψ a k ) ] .
(56)
On the right-hand side of Eqs. (51)–(53), coefficients for normalization of each weights are included and the time average in the steady state of turbulence is used instead of the ensemble average to evaluate . Other curves labeled ( W Γ a , N 0 , W q a , N 0 , W Y a , N 0 ) and ( W Γ a , L 0 , W q a , L 0 , W Y a , L 0 ) represent the nonlinear and quasilinear weights obtained by
W Γ a , N 0 ( k y ) = T i 2 e 2 n e v t i Γ ̃ a k x = 0 , k y , N | ϕ k x = 0 , k y , N | 2 ,
(57)
W q a , N 0 ( k y ) = T i e 2 n e v t i q ̃ a k x = 0 , k y , N | ϕ k x = 0 , k y , N | 2 ,
(58)
W Y a , N 0 ( k y ) = T i R 0 e 2 n e v t i Y ̃ a k x = 0 , k y , N | ϕ k x = 0 , k y , N | 2 ,
(59)
and
W Γ a , L 0 ( k y ) = T i 2 e 2 n e v t i Γ ̃ a k x = 0 , k y , L | ϕ k x = 0 , k y , L | 2 ,
(60)
W q a , L 0 ( k y ) = T i e 2 n e v t i q ̃ a k x = 0 , k y , L | ϕ k x = 0 , k y , L | 2 ,
(61)
W Y a , L 0 ( k y ) = T i R 0 e 2 n e v t i Y ̃ a k x = 0 , k y , L | ϕ k x = 0 , k y , L | 2 ,
(62)
where only the kx = 0 modes are kept instead of summing over kx. Here, the subscripts L and N denote the results from the linear and nonlinear simulations, respectively. On the right-hand side of Eqs. (59)–(61), denotes the surface average and the ratios of Γ ̃ a k x = 0 , k y , L , q ̃ a k x = 0 , k y , L , and Y ̃ a k x = 0 , k y , L to ϕ k x = 0 , k y , L | 2 are evaluated for the linear unstable mode with the wavenumbers kx = 0 and ky.

The quasilinear weights for the kx = 0 modes show a good agreement with those of the nonlinear weights for kx = 0 in the linearly unstable wavenumber region 0.05 k y ρ t i 0.5 [see Fig. 3(a)]. In the areas colored in sky blue, we have N e k y < 0 and N i k y < 0, which indicate that the entropy of the fluctuation at ky is transferred to other wavenumber regions through nonlinear interaction. Both of the nonlinear weights obtained by keeping the only kx = 0 modes and by summing over kx agree well with each other in the colored regions. Thus, in these regions, the nonlinear weights including all kx's are well approximated by the quasilinear weights for kx = 0 within an error margin of 30% or less. We find that more than 80% of the total values of the transport fluxes and the energy exchange over the whole wavenumber space can be accounted for by contributions from the colored wavenumber regions. Therefore, nonlinear simulation results of the transport fluxes and the energy exchange can be effectively predicted from the quasilinear weights for kx = 0 multiplied by the squared potential amplitude. The model for predicting the magnitude of the potential amplitude is not treated in the present work, although there are many analytical and numerical studies to parameterize the saturation amplitude.16–18,20,21 The results shown above indicate that it is possible to construct a quasilinear model which can accurately predict both of the turbulent transport fluxes and the turbulent energy exchange.

In this study, the effect of ITG turbulence on the energy exchange between electrons and ions in tokamak plasmas is investigated. The ITG turbulence is found to be dominant in the energy exchange in equithermal or high-temperature plasmas in which collisional energy exchange is negligibly small. It is also shown that the direction of net energy transfer can be opposite to that of the collisional one from hotter to colder particle species, since ITG turbulence transfers energy from ions to electrons, even when ions are colder than electrons. This result does not contradict with the second law of thermodynamics because the entropy balance is still maintained by the entropy production, mainly due to the ion heat transport from hot to cold regions. Therefore, the ITG turbulence is anticipated to prevent energy transfer from alpha-heated electrons to ions, which is considered a primary ion heating mechanism in future reactors. The wavenumber spectral analysis reveals that the main physical mechanisms of turbulent energy exchange are the cooling of ions in the B-curvature drift motion and the heating of electrons streaming along the field line, which are caused by the perpendicular and parallel components of the turbulent electric field, respectively. In particular, the perpendicular cooling of ions is closely linked to the physical mechanism of ITG instability which drives the ion heat flux.

Since the effect of turbulence on the energy exchange between electrons and ions can possibly overcome that of Coulomb collisions, the turbulent energy exchange as well as the turbulent transport fluxes of particles and heat should be taken into account for predicting global profiles of the density and temperature in future fusion reactors. In order to examine the predictability of the quasilinear model of the energy exchange and transport fluxes, the quasilinear weights of these turbulent quantities normalized by the squared potential are estimated as functions of the wavenumber by linear simulations. It is found that the quasilinear weights of both energy exchange and fluxes agree with the nonlinear simulation results within an error margin of 30% in the wavenumber region, where more than 80% of total energy exchange and fluxes are covered. This indicates that we can construct the quasilinear model which is valid for predicting the energy exchange as well as the transport fluxes. Although this study has not investigated the saturation model, it would not be difficult to incorporate turbulent energy exchange in existing codes that predict fluxes with a quasilinear model.18,21

The entropy balance in the linearly growing and nonlinearly saturated phases of the ITG modes is examined to understand the conditions for the quasilinear weights estimated from the linear analysis to be applicable to describing the steady state of turbulence. The analogy of the entropy balance equation to the Landau equation for weakly nonlinear fluid dynamic systems is noted in that the time-derivative term in the linear phase should be replaced by the nonlinear term in the steady state for the quasilinear weights to keep the same ratios in both phases. Then, we speculate that the nonlinear entropy transfer term in the steady state should be negative in the linearly unstable wavenumber region for the quasilinear model to be valid. This speculation is confirmed from the ITG turbulence simulation showing that the quasilinear weights agree well with the nonlinear results in the wavenumber regions where the nonlinear entropy transfer term is negative.

It is conjectured from results of this work that energy is generally transferred by turbulence from a particle species with larger entropy production due to particle and heat transport driven by the instability to other species regardless of which species is hotter. To verify this conjecture, turbulent energy exchange and its predictability in cases of other instabilities such as those driven by trapped electrons and electron temperature gradient remain subjects of future studies.

Furthermore, in principle, the theory and simulation methods for turbulent energy exchange in this study are expected to be applicable or extendable to non-tokamak systems, such as helical systems (stellarators and heliotrons).36 The study of turbulent energy exchange and quasilinear weights in such general 3D systems is also the subject of future research.

This study was supported in part by the JSPS Grants-in-Aid for Scientific Research (Grant Nos. 19H01879 and 24K07000) and in part by the NIFS Collaborative Research Program (Grant No. NIFS23KIPT009). This work was also supported by JST SPRING (Grant No. JPMJSP2108) and the NINS program of Promoting Research by Networking among Institutions (Grant No. 01422301). Simulations in this work were performed on “Plasma Simulator” (NEC SX-Aurora TSUBASA) of NIFS with the support and under the auspices of the NIFS Collaboration Research program (Grant No. NIFS23KIPT009).

The authors have no conflicts to disclose.

T. Kato: Conceptualization (equal); Data curation (lead); Formal analysis (equal); Investigation (lead); Methodology (equal); Writing—original draft (lead); Writing—review & editing (lead). H. Sugama: Conceptualization (equal); Formal analysis (equal); Investigation (supporting); Methodology (equal); Supervision (lead); Writing—original draft (supporting); Writing—review & editing (supporting). T.-H. Watanabe: Methodology (equal); Writing—review & editing (supporting). M. Nunami: Project administration (lead); Writing—review & editing (supporting).

The data that support the findings of this study are available from the corresponding authors upon reasonable request.

In this work, we follow the local gyrokinetic theory based on the WKB (or ballooning) representation1 and assume that background E × B velocity is on the order of ( ρ t a / L ) v t a = δ v t a (the same order as that of the diamagnetic velocity v * a) and that the gradient scale length of the background electric field is the same as the equilibrium scale length L. Under these assumptions, the Doppler frequency shift due to the background E × B velocity k · v E is considered to be on the same order as the diamagnetic frequency k · v * a although the background E × B velocity shear effect is estimated by Δ r ( v E / r ) ( r × b ) · f ̂ a / | r | ρ t a ( δ v t a / L ) k f ̂ a ( δ v t a / L ) f ̂ a ( δ k v * a ) f ̂ a, the magnitude of which is smaller than terms included in the gyrokinetic equation, Eq. (5), by a factor of δ. Here, Δ r is the radial correlation length of the fluctuation which is estimated as the typical perpendicular wavelength 1 / k ρ t a.

If the background E × B flow is larger and/or the gradient scale length of the background electric field is smaller than assumed above, the background E × B flow shear effect cannot be neglected and both neoclassical and turbulent processes of transport and energy exchange are considered to have Er-dependence different from the results in this paper. Also, the coupling of neoclassical and turbulent processes which influences the Er-dependence may occur in a global model including the sheared Er profile. A nonlinear gyrokinetic equation with large flow velocities on the order of the ion thermal speed is derived in Ref. 37.

Based on the ordering assumption described above, the sum of the last two terms on the right-hand side of Eq. (18) can be rewritten as follows:9 
e a Γ a turb Φ r + Q a turb = e a k Re d 3 v h a k * ( t + i k · v E ) ψ a k = e a k Re d 3 v [ e i k · v E h a k ] * t [ e i k · v E ψ a k ] .
(A1)
As pointed out in Ref. 9, the background radial electric field E r = Φ / r enters the gyrokinetic equation, Eq. (5) only in the form of the Doppler shift ( / t + i k · v E ) and does not appear explicitly in Eqs. (6)–(9). Therefore, for the solutions of Eq. (5) and Eqs. (7)–(9), e i k · v E h a k and e i k · v E ψ a k are independent of Er. Then, even though each of the first and second terms on the left-hand side of Eq. (A1) depends on Er, their sum does not. In the present paper, we evaluate Q a turb by the gyrokinetic simulation for Er = 0 although we should note that the results of Q a turb for Er = 0 are equivalent to those of e a Γ a turb E r + Q a turb for any Er. We also find from Eq. (14) that the right-hand side of the entropy balance equation, Eq. (12), contains this sum of the terms e a Γ a turb E r + Q a turb, which takes a value independent of Er.
In this appendix, the derivation of Eq. (24) is presented. We first note
e a d 3 v h a k * ψ a k t = t e a d 3 v ψ a k h a k * e a d 3 v ψ a k h a k * t ,
(B1)
where denotes a double average over the ensemble and the flux surface. The first term of the right-hand side in Eq. (B1) is negligible because of the time transport scale ordering ( / t ) / = O ( δ 2 ). Substituting Eq. (5) into Eq. (B1), we obtain
e a k d 3 v ψ a k h a k * t = e a k d 3 v ψ a k ( v b · h a k * + i k · v d a h a k * c B k + k = k [ b · ( k × k ) ] ψ a k * h a k * C a G K * ) .
(B2)
On the right-hand side of Eq. (B2), a pure imaginary term disappears and a time derivative term is neglected again because of the transport timescale ordering. We immediately see that Eqs. (26) and (28) corresponds to the second and fourth terms on the right-hand side of Eq. (B2). Equation (27) is also derived from the third term by k + k = k [ b · ( k × k ) ] = k + k = k [ b · ( k × k ) ].
When phase space variables (x, w = m a v 2 / 2 , μ = m a v 2 / 2 B) are used, the integral over velocity space is described as follows:
d 3 v = B m a 2 σ d ξ 0 d w 0 w / B d μ 1 | v | ,
(B3)
where ξ is the gyrophase and σ = v / | v | = ± 1. The first term can be written as follows:
e a d 3 v ψ a k v b · h a k * = e a m a 2 σ d ξ 0 d w 0 w / B d μ ( B · ( σ h a k * ψ a k ) ) e a d 3 v ( h a k * v b · ψ a k ) = e a m a 2 ( B · ) ( d ξ 0 d w 0 w / B d μ ( σ σ h a k * ψ a k ) ) e a d 3 v ( h a k * v b · ψ a k ) .
(B4)
Since the integral range of μ depends on position x, the boundary of the integral range needs to be taken into account when moving outside of the integral on the right-hand side of Eq. (B4). However, when μ = w / B, it is a bounce point ( v = 0) and the integral function σ σ h a k * ψ a k at the point is zero. Therefore, the boundary effect does not affect the integral calculations. The flux surface average of the first term of the right-hand side in Eq. (B4) vanishes due to B · = 0.
Then, we obtain
e a d 3 v ψ a k v b · h a k * = e a d 3 v h a k * v b · ψ a k ,
(B5)
which gives Eq. (25).
1.
T. M.
Antonsen
, Jr.
and
B.
Lane
,
Phys. Fluids
23
,
1205
(
1980
).
2.
P. J.
Catto
,
W. M.
Tang
, and
D. E.
Baldwin
,
Plasma Phys.
23
,
639
(
1981
).
3.
E. A.
Frieman
and
L.
Chen
,
Phys. Fluids
25
,
502
(
1982
).
4.
A. M.
Dimits
,
G.
Bateman
,
M. A.
Beer
,
B. I.
Cohen
,
W.
Dorland
,
G. W.
Hammett
,
C.
Kim
,
J. E.
Kinsey
,
M.
Kotschenreuther
,
A. H.
Kritz
,
L. L.
Lao
,
J.
Mandrekas
,
W. M.
Nevins
,
S. E.
Parker
,
A. J.
Redd
,
D. E.
Shumaker
,
R.
Sydora
, and
J.
Weiland
,
Phys. Plasmas
7
,
969
(
2000
).
5.
J. A.
Krommes
,
Ann. Rev. Fluid Mech.
44
,
175
(
2012
).
6.
X.
Garbet
,
Y.
Idomura
,
L.
Villard
, and
T.-H.
Watanabe
,
Nucl. Fusion
50
,
043002
(
2010
).
7.
Y.
Idomura
,
T.-H.
Watanabe
, and
H.
Sugama
,
C. R. Phys.
7
,
650
(
2006
).
8.
W.
Horton
,
Turbulent Transport in Magnetized Plasmas
, 2nd ed. (
World Scientific
,
Singapore
,
2018
), Chap. 12.
9.
H.
Sugama
,
M.
Okamoto
,
W.
Horton
, and
M.
Wakatani
,
Phys. Plasmas
3
,
2379
(
1996
).
10.
H.
Sugama
,
T.-H.
Watanabe
, and
M.
Nunami
,
Phys. Plasmas
16
,
112503
(
2009
).
11.
F. L.
Hinton
and
R. E.
Waltz
,
Phys. Plasmas
13
,
102301
(
2006
).
12.
R. E.
Waltz
and
G. M.
Staebler
,
Phys. Plasmas
15
,
014505
(
2008
).
13.
J.
Candy
,
Phys. Plasmas
20
,
082503
(
2013
).
14.
M.
Barnes
,
F. I.
Parra
, and
W.
Dorland
,
Phys. Rev. Lett.
109
,
185003
(
2012
).
15.
Y.
Kawazura
,
M.
Barnes
, and
A. A.
Schekochihin
,
Proc. Natl. Acad. Sci.
116
,
771
(
2018
).
16.
A.
Casati
,
C.
Bourdelle
,
X.
Garbet
,
F.
Imbeaux
,
J.
Candy
,
F.
Clairet
,
G.
Dif-Pradalier
,
G.
Falchetto
,
T.
Gerbaud
,
V.
Grandgirard
,
Ö. D.
Gürcan
,
P.
Hennequin
,
J.
Kinsey
,
M.
Ottaviani
,
R.
Sabot
,
Y.
Sarazin
,
L.
Vermare
, and
R. E.
Waltz
,
Nucl. Fusion
49
,
085012
(
2009
).
17.
C.
Bourdelle
,
J.
Citrin
,
B.
Baiocchi
,
A.
Casati
,
P.
Cottier
,
X.
Garbet
, and
F.
Imbeaux
, and
JET Contributors
,
Plasma Phys. Controlled Fusion
58
,
014036
(
2016
).
18.
G. M.
Staebler
,
E. A.
Belli
,
J.
Candy
,
J. E.
Kinsey
,
H.
Dudding
, and
B.
Patel
,
Nucl. Fusion
61
,
116007
(
2021
).
19.
E.
Narita
,
M.
Honda
,
M.
Nakata
,
M.
Yoshida
, and
N.
Hayashi
,
Nucl. Fusion
61
,
116041
(
2021
).
20.
S.
Toda
,
M.
Nunami
, and
H.
Sugama
,
Plasma Phys. Controlled Fusion
64
,
085001
(
2022
).
21.
S. E.
Parker
,
C. S.
Haubrich
,
S.
Tirkas
,
Q.
Cai
, and
Y.
Chen
,
Plasma
6
,
611
(
2023
).
22.
M.
Barnes
,
P.
Abiuso
, and
W.
Dorland
,
J. Plasma Phys.
84
,
905840306
(
2018
).
23.
M.
Nakata
,
T.-H.
Watanabe
, and
H.
Sugama
,
Phys. Plasmas
19
,
022303
(
2012
).
24.
S.
Maeyama
,
A.
Ishizawa
,
T.-H.
Watanabe
,
M.
Nakata
,
N.
Miyato
,
M.
Yagi
, and
Y.
Idomura
,
Phys. Plasmas
21
,
052301
(
2014
).
25.
T.-H.
Watanabe
and
H.
Sugama
,
Nucl. Fusion
46
,
24
(
2006
).
26.
T.-H.
Watanabe
,
H.
Sugama
, and
S.
Ferrando-Margalet
,
Nucl. Fusion
47
,
9
(
2007
).
27.
T.-H.
Watanabe
,
H.
Sugama
, and
S.
Ferrando-Margalet
,
Phys. Rev. Lett.
100
,
195002
(
2007
).
28.
M. A.
Beer
,
S. C.
Cowley
, and
G. W.
Hammett
,
Phys. Plasmas
2
,
2687
(
1995
).
29.
R. D.
Hazeltine
and
J. D.
Meiss
,
Plasma Confinement
(
Addison-Wesley
,
Redwood City, CA
,
1992
), Chap. 7.10.
30.
L. D.
Landau
and
E. M.
Lifshitz
,
Fluid Mechanics
(
Pergamon
,
New York
,
1975
), Chap. 3.
31.
P. G.
Drazin
and
W. H.
Reid
,
Hydrodynamic Stability
(
Cambridge University Press
,
Cambridge
,
1981
), Chap. 7.
32.
W. W.
Lee
,
J. Comput. Phys.
72
,
243
(
1987
).
33.
A.
Lenard
and
I. B.
Bernstein
,
Phys. Rev.
112
,
1456
(
1958
).
34.
S.
Maeyama
,
A.
Ishizawa
,
T.-H.
Watanabe
,
N.
Nakajima
,
S.
Tsuji-Iio
, and
H.
Tsutsui
,
Comput. Phys. Commun.
184
,
11
(
2013
).
35.
A. E.
White
,
L.
Schmitz
,
G. R.
McKee
,
C.
Holland
,
W. A.
Peebles
,
T. A.
Carter
,
M. W.
Shafer
,
M. E.
Austin
,
K. H.
Burrell
,
J.
Candy
,
J. C.
DeBoo
,
E. J.
Doyle
,
M. A.
Makowski
,
R.
Prater
,
T. L.
Rhodes
,
G. M.
Staebler
,
G. R.
Tynan
,
R. E.
Waltz
, and
G.
Wang
,
Phys. Plasmas
15
,
056116
(
2008
).
36.
M.
Wakatani
,
Stellarator and Heliotron Devices
(
Oxford University Press
,
Oxford
,
1998
), Chap. 2.
37.
H.
Sugama
and
W.
Horton
,
Phys. Plasmas
5
,
2560
2573
(
1998
).