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 $ \u2207 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.

## I. INTRODUCTION

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 plasmas^{14} 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.

## II. THEORETICAL MODEL

### A. Gyrokinetic equations

*F*of the position vector

_{a}**, velocity vector**

*x***, and time**

*v**t*can be written as $ F a ( x , v , t ) = f a ( x , v , t ) + f \u0302 a ( x , v , t )$, where

*f*and $ f \u0302 a$ represent the ensemble-averaged and fluctuation parts, respectively, and the subscript

_{a}*a*denotes the particle species. The space-time scales of variations in the ensemble-averaged part

*f*are much larger than those of the fluctuation part $ f \u0302 a$ so that the ensemble average can also be regarded as the local space-time average. In the gyrokinetic theory, perturbations such as $ f \u0302 a$ are assumed to satisfy the gyrokinetic ordering

_{a}*m*,

_{a}*e*,

_{a}*T*,

_{a}**,**

*B**c*, and $ v t a \u2261 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 \u2225 , \u2009 k \u22a5$,

*ω*, and

*L*, respectively.

^{29}

*δ*, the ensemble-averaged distribution function

*f*is assumed to be the local Maxwellian

_{a}*f*, which is given in terms of the background density

_{Ma}*n*and temperature

_{a}*T*by

_{a}^{10}and $ E = \u2212 \u2207 \Phi $ is the background electric field. Here, the background

*E*×

*B*flow is assumed to be $ v E \u223c \delta 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 \u22a5$ is regarded as a function of time

*t*and phase space variables (

**, $ w = m a v 2 / 2 , \u2009 \mu = m a v \u22a5 2 / 2 B$), where $ v = | v |$ and $ v \u22a5 = | v \u2212 v \u2225 b |$. The gyrophase-averaged perturbed potential function $ \psi a k \u22a5$ is defined in terms of the perturbed electrostatic potential $ \varphi k \u22a5$ and the perturbed vector potential $ A k \u22a5$ as follows:**

*x**J*

_{0}and

*J*

_{1}are the zeroth- and first-order Bessel functions, respectively, $ A \u2225 k \u22a5 = b \xb7 A k \u22a5$, and $ B \u2225 k \u22a5 = i b \xb7 ( k \u22a5 \xd7 A k \u22a5$). The perturbed electric and magnetic fields are determined by Poisson's equation and Ampère's law in the gyrokinetic form as follows:

### B. Energy and entropy balance equations

^{*}denotes the complex conjugate and $ n a k \u22a5 \u2261 \u222b d 3 v f a k \u22a5$ is the density perturbation with the perpendicular wavenumber vector $ k \u22a5$. Multiplying Eq. (5) by $ h a k \u22a5 * / f M a$ and taking the ensemble and flux-surface averages, we can derive the entropy balance equation as follows:

*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

*δ*. Therefore, we can interpret $ Q a turb$ as the energy exchange between different plasma species. The last term in Eq. (12)

*D*represents collisional dissipation written as

_{a}^{9}

*p*and $ \pi a$ are the pressure and viscosity tensor, respectively. Here, $ V \u2032$ is the derivative of the volume

_{a}*V*inside a magnetic surface with respect to the minor radius

*r*. The radial particle and heat fluxes denoted by $ \Gamma a$ and

*q*, respectively, are written as follows:

_{a}*ν*is the electron-ion collision frequency, $ ln \u2009 \Lambda $ is the Coulomb logarithm,

_{e}**R**

_{e}is the collisional friction force for electrons, and $ ( u e \u2212 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 *E _{r}* 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

*E*= 0. Therefore, we henceforth discuss turbulent energy exchange with

_{r}*E*= 0.

_{r}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.

### C. Spectral analysis of turbulent energy exchange in wavenumber space

*E*= 0 is used. The perturbed fields $ G a k \u22a5$ and the perturbed currents $ j a k \u22a5$ at the gyrocenter position are defined by

_{r}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.

*a*but it represents the energy transfer in the wavenumber space through nonlinear interactions between different modes. Thus, $ Q a \psi k \u22a5 turb$ influences the profile of the total wavenumber spectrum $ Q a k \u22a5 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 \u2225 k \u22a5 turb , \u2009 Q a B k \u22a5 turb$, and $ Q a C k \u22a5 turb$ over the whole wavenumber space. In particular, $ Q a \u2225 k \u22a5 turb , \u2009 Q a B k \u22a5 turb$ are the contributions to Joule heating (cooling) via currents parallel and perpendicular to the background magnetic field.

### D. Quasilinear model and entropy balance

*Y*can be rewritten as follows:

_{a}^{13}This expression of

*Y*in Eq. (40) can be used for both linear and nonlinear cases and rigorously satisfies $ \u2211 a Y a = 0$ even in non-steady states. The relation between $ Q a turb$ and

_{a}*Y*is given from Eqs. (15) and (40) as follows:

_{a}*a*.

In the quasilinear model, the transport fluxes divided by the squared potential, $ \Gamma a k \u22a5 turb / \u27e8 \u27e8 | \varphi k \u22a5 | 2 \u27e9 \u27e9$ and $ q a k \u22a5 turb / \u27e8 \u27e8 | \varphi k \u22a5 | 2 \u27e9 \u27e9$, 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 \u22a5 / \u27e8 \u27e8 | \varphi k \u22a5 | 2 \u27e9 \u27e9$ from the linear calculation. When this model is valid, the values of the terms on the right-hand side of Eq. (46) divided by $ \u27e8 \u27e8 | \varphi k \u22a5 | 2 \u27e9 \u27e9$ should not change whether linear or nonlinear simulations are performed to evaluate them. On the other hand, when divided by $ \u27e8 \u27e8 | \varphi k \u22a5 | 2 \u27e9 \u27e9$, 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 \u22a5 / \u27e8 \u27e8 | \varphi k \u22a5 | 2 \u27e9 \u27e9$, 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 $ \u27e8 \u27e8 | \varphi k \u22a5 | 2 \u27e9 \u27e9$ in the linear case should be the same as the value of $ \u2212 N a k \u22a5 / \u27e8 \u27e8 | \varphi k \u22a5 | 2 \u27e9 \u27e9$ in the nonlinear case in order to keep the balance between the left- and right-hand sides of Eq. (46).

^{30,31}

*γ*are the amplitude of the dominant mode and its linear growth rate, respectively, and $ \alpha | 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 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.

## III. SIMULATION RESULTS

### A. Simulation settings

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 \u2225 , \mu )$ 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 (*k _{x}*,

*k*). 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.

_{y}^{26,27}

The flux tube coordinates for a low-*β*, large aspect ratio, axisymmetric torus with concentric circular cross sections, $ x = r \u2212 r 0 , y = r 0 / q 0 ( q ( r ) \theta \u2212 \zeta ) ,$ and $ z = \theta $ 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 \u22a5 = ( k x + s \u0302 z k y ) e r + k y e \theta $. Here, *k _{x}* and

*k*are wavenumbers in the directions of $ \u2207 x$ and $ \u2207 y$, respectively, while $ e r$ and $ e \theta $ are unit vectors parallel to $ \u2207 r$ and $ \u2207 \theta $, respectively.

_{y}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 *R*_{0} 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 $ \beta i = 1 \xd7 10 \u2212 4$ for which the electrostatic approximation is valid. In simulations performed here, $ B \u2225 k \u22a5$ is neglected, although $ A \u2225 k \u22a5$ is retained in order to avoid numerical difficulty due to very rapid electrostatic waves called the *ω _{H}* mode.

^{32}The Lenard–Bernstein collision operator

^{33}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 $ \nu i i * \u2261 R 0 q \nu i i / ( \u03f5 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.

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$ | 0 |

Normalized density gradient | $ R 0 / L n a$ | 2.22 |

Mass ratio | $ m e / m i$ | $ 5.446 \xd7 10 \u2212 4$ |

Charge ratio | $ e e / e i$ | −1 |

Ion beta value | $ \beta i = 4 \pi n i T i / B 2$ | $ 1 \xd7 10 \u2212 4$ |

Normalized collision frequency | $ \nu i i * \u2261 R 0 q 0 \nu i i / ( \u03f5 r 3 / 2 v t i )$ | 0.068 |

Temperature ratio | $ T e / T i$ | $ 0.8 \u2013 1.5$ |

Inverse aspect ratio | $ \epsilon r$ | 0.18 |

Safety factor | q_{0} | 1.4 |

Magnetic shear | $ s \u0302 = ( 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$ | 0 |

Normalized density gradient | $ R 0 / L n a$ | 2.22 |

Mass ratio | $ m e / m i$ | $ 5.446 \xd7 10 \u2212 4$ |

Charge ratio | $ e e / e i$ | −1 |

Ion beta value | $ \beta i = 4 \pi n i T i / B 2$ | $ 1 \xd7 10 \u2212 4$ |

Normalized collision frequency | $ \nu i i * \u2261 R 0 q 0 \nu i i / ( \u03f5 r 3 / 2 v t i )$ | 0.068 |

Temperature ratio | $ T e / T i$ | $ 0.8 \u2013 1.5$ |

Inverse aspect ratio | $ \epsilon r$ | 0.18 |

Safety factor | q_{0} | 1.4 |

Magnetic shear | $ s \u0302 = ( 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.

Domain sizes and resolved perpendicular wavenumbers |

$ \u2212 64.10 \rho t i \u2264 x \u2264 64.10 \rho t i$ |

$ \u2212 62.83 \rho t i \u2264 y \u2264 62.83 \rho t i$ |

$ \u2212 \pi \u2264 z \u2264 \pi $ |

$ \u2212 4 v t a \u2264 v \u2225 \u2264 4 v t a$ |

$ 0 \u2264 \mu B 0 / T a \u2261 m a v \u22a5 2 / 2 T a \u2264 8$ |

$ \u2212 4.70 \rho t i \u2212 1 \u2264 k x \u2264 4.70 \rho t i \u2212 1 , ( k x , min = 0.049 \rho t i \u2212 1 )$ |

$ \u2212 1.55 \rho t i \u2212 1 \u2264 k y \u2264 1.55 \rho t i \u2212 1 , ( k y , min = 0.050 \rho t i \u2212 1 )$ |

Grid points in $ ( x , y , z , v \u2225 , \mu )$ |

$ 288 \xd7 96 \xd7 64 \xd7 64 \xd7 32$ |

$ 288 \xd7 96 \xd7 256 \xd7 64 \xd7 32$ |

Domain sizes and resolved perpendicular wavenumbers |

$ \u2212 64.10 \rho t i \u2264 x \u2264 64.10 \rho t i$ |

$ \u2212 62.83 \rho t i \u2264 y \u2264 62.83 \rho t i$ |

$ \u2212 \pi \u2264 z \u2264 \pi $ |

$ \u2212 4 v t a \u2264 v \u2225 \u2264 4 v t a$ |

$ 0 \u2264 \mu B 0 / T a \u2261 m a v \u22a5 2 / 2 T a \u2264 8$ |

$ \u2212 4.70 \rho t i \u2212 1 \u2264 k x \u2264 4.70 \rho t i \u2212 1 , ( k x , min = 0.049 \rho t i \u2212 1 )$ |

$ \u2212 1.55 \rho t i \u2212 1 \u2264 k y \u2264 1.55 \rho t i \u2212 1 , ( k y , min = 0.050 \rho t i \u2212 1 )$ |

Grid points in $ ( x , y , z , v \u2225 , \mu )$ |

$ 288 \xd7 96 \xd7 64 \xd7 64 \xd7 32$ |

$ 288 \xd7 96 \xd7 256 \xd7 64 \xd7 32$ |

### B. Heat flux, turbulent energy exchange, and entropy balance

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 \rho t i , k y \rho 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.

### C. Comparison between collisional and turbulent energy exchanges

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 \u2009 \Lambda = 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 \u2009 ( keV ) , \u2009 2.0 \xd7 10 19 \u2009 ( m \u2212 3 )$] and for the [ $ 3.0 \u2009 ( keV ) , \u2009 2.25 \xd7 10 20 \u2009 ( 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 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 , \u2009 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 \xd7 10 19 \u2009 ( m \u2212 3 )$ and $ T i = 0.9 \u2009 ( 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 *T _{e}* >

*T*. 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.

_{i}### D. Spectral analysis of the turbulent energy exchange

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 \u22a5 turb$ and $ Q i C k \u22a5 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 \u2225 k \u22a5 turb > 0 ( a = e , i )$, and negative values of the components $ Q a B k \u22a5 turb < 0$ $ ( a = e , i )$ caused by the product of the perpendicular field and the $ \u2207 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 \u2225 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 \u2225 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.

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

### E. Correlation between results from linear and nonlinear simulations

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 *k _{y}* = 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.

*k*= 0 modes are kept instead of summing over

_{x}*k*. Here, the subscripts

_{x}*L*and

*N*denote the results from the linear and nonlinear simulations, respectively. On the right-hand side of Eqs. (59)–(61), $ \u27e8 \cdots \u27e9$ denotes the surface average and the ratios of $ \u27e8 \Gamma \u0303 a k x = 0 , k y , L \u27e9 , \u2009 \u27e8 q \u0303 a k x = 0 , k y , L \u27e9$, and $ \u27e8 Y \u0303 a k x = 0 , k y , L \u27e9$ to $ \u27e8 \varphi k x = 0 , k y , L | 2 \u27e9$ are evaluated for the linear unstable mode with the wavenumbers

*k*= 0 and

_{x}*k*.

_{y}The quasilinear weights for the *k _{x}* = 0 modes show a good agreement with those of the nonlinear weights for

*k*= 0 in the linearly unstable wavenumber region $ 0.05 \u2264 k y \rho t i \u2264 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

_{x}*k*is transferred to other wavenumber regions through nonlinear interaction. Both of the nonlinear weights obtained by keeping the only

_{y}*k*= 0 modes and by summing over

_{x}*k*agree well with each other in the colored regions. Thus, in these regions, the nonlinear weights including all

_{x}*k*'s are well approximated by the quasilinear weights for

_{x}*k*= 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

_{x}*k*= 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.

_{x}^{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.

## IV. CONCLUSIONS AND DISCUSSION

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 $ \u2207 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.

## ACKNOWLEDGMENTS

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).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

### APPENDIX A: DISCUSSION OF THE EFFECT OF ELECTRIC FIELD

In this work, we follow the local gyrokinetic theory based on the WKB (or ballooning) representation^{1} and assume that background *E* × *B* velocity is on the order of $ ( \rho t a / L ) v t a = \delta 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 \u22a5 \xb7 v E$ is considered to be on the same order as the diamagnetic frequency $ k \u22a5 \xb7 v * a$ although the background *E* × *B* velocity shear effect is estimated by $ \Delta r ( \u2202 v E / \u2202 r ) ( \u2207 r \xd7 b ) \xb7 \u2207 f \u0302 a / | \u2207 r | \u223c \rho t a ( \delta v t a / L ) k \u22a5 f \u0302 a \u223c ( \delta v t a / L ) f \u0302 a \u223c ( \delta k \u22a5 v * a ) f \u0302 a$, the magnitude of which is smaller than terms included in the gyrokinetic equation, Eq. (5), by a factor of *δ*. Here, $ \Delta r$ is the radial correlation length of the fluctuation which is estimated as the typical perpendicular wavelength $ 1 / k \u22a5 \u223c \rho 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 *E _{r}*-dependence different from the results in this paper. Also, the coupling of neoclassical and turbulent processes which influences the

*E*-dependence may occur in a global model including the sheared

_{r}*E*profile. A nonlinear gyrokinetic equation with large flow velocities on the order of the ion thermal speed is derived in Ref. 37.

_{r}^{9}

*E*. Then, even though each of the first and second terms on the left-hand side of Eq. (A1) depends on

_{r}*E*, their sum does not. In the present paper, we evaluate $ Q a turb$ by the gyrokinetic simulation for

_{r}*E*= 0 although we should note that the results of $ Q a turb$ for

_{r}*E*= 0 are equivalent to those of $ e a \Gamma a turb E r + Q a turb$ for any

_{r}*E*. 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 \Gamma a turb E r + Q a turb$, which takes a value independent of

_{r}*E*.

_{r}### APPENDIX B: DERIVATION OF EQ. (24)

**, $ w = m a v 2 / 2 , \u2009 \mu = m a v \u22a5 2 / 2 B$) are used, the integral over velocity space is described as follows:**

*x**ξ*is the gyrophase and $ \sigma = v \u2225 / | v \u2225 | = \xb1 1$. The first term can be written as follows:

*μ*depends on position

**, the boundary of the integral range needs to be taken into account when moving $\u2207$ outside of the integral on the right-hand side of Eq. (B4). However, when $ \mu = w / B$, it is a bounce point ( $ v \u2225 = 0$) and the integral function $ \u2211 \sigma \sigma h a k \u22a5 * \psi a k \u22a5$ 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 $ \u27e8 B \xb7 \u2207 \cdots \u27e9 = 0$.**

*x*## REFERENCES

*Turbulent Transport in Magnetized Plasmas*

*Plasma Confinement*

*Fluid Mechanics*

*Hydrodynamic Stability*

*Stellarator and Heliotron Devices*