Collisionless kinetic-fluid closure and its application to the three-mode ion temperature gradient driven system

A novel closure model is presented to give a set of fluid equations which describe a collisionless kinetic system. In order to take account of the time reversal symmetry of the collisionless kinetic equation, the new closure model relates the parallel heat flux to the temperature and the parallel flow in terms of the real-valued coefficients in the unstable wave number space. Effects of the closure model on turbulence saturation and anomalous transport are investigated based on kinetic and fluid entropy balances. When the closure model is applied to the three-mode ion temperature gradient ~ITG! driven system, the fluid system of equations reproduces the exact nonlinear kinetic solution found by Watanabe, Sugama, and Sato @Phys. Plasmas7, 984 ~2000!#. Oscillatory behaviors and initial amplitude dependence of other numerical kinetic solutions of the three-mode ITG problem can also be accurately described by the fluid system. © 2001 American Institute of Physics. @DOI: 10.1063/1.1367319 #


I. INTRODUCTION
Many theoretical studies have been done in order to understand anomalous transport processes observed in magnetically confined plasmas.Recently, predictions of the anomalous transport coefficients have come to be made by firstprinciple simulations of the ion temperature gradient ͑ITG͒ turbulence 1 based on gyrokinetic 2-5 and gyrofluid [6][7][8][9] ͑or gyro-Landau-fluid͒ models.Gyrokinetic simulations directly solve gyrokinetic equations for gyrocenter distribution functions and the Maxwell equations for electromagnetic fields while gyrofluid simulations solve a truncated system of fluid equations instead of the gyrokinetic equation.The gyrofluid simulations are attractive in that less computer memory and time are consumed although they depend on a closure relation between the highest-order fluid variable and lower-order variables which is assumed to be valid even for collisionless kinetic systems.Since conventional gyrofluid closure models 10,11 are derived so as to accurately reproduce kinetic dispersion relations for linear modes, the gyrokinetic and gyrofluid simulations show good agreements in their linear results.However, there exist disagreements in their nonlinear results such as the saturated fluctuation levels and the turbulent transport coefficients. 12,13he three-mode ITG system 2 is one of the simplest examples to show these agreements and disagreements between the kinetic and Landau-fluid models, in which effects of finite gyroradii are neglected by taking the long-wavelength limit.Recently, an exact nonlinear solution of the threemode ITG problem was found by Watanabe, Sugama, and Sato ͑WSS͒. 14The nonlinear solution of the three-mode ITG problem shows periodic oscillation of the amplitude of the electrostatic potential.On the other hand, the solution of the kinetic-fluid equations using the conventional linear closure models fail to reproduce this oscillation and instead give the amplitude saturation at a certain value. 15Also, it is important to note that the exact kinetic solution is symmetric in time reversal while the solution of the conventional kinetic-fluid equations is not symmetric due to the dissipative closure term.
In this paper, a novel collisionless kinetic-fluid closure model is presented, which takes into account the time reversal symmetry of the collisionless kinetic equation by including nondissipative closure terms.It is shown that, when the new closure model is applied to the three-mode ITG problem, the exact kinetic solution by WSS is reproduced as a solution of the fluid equations.Recently, Mattor and Parker presented a nonlinear closure model, 15 which also describes the behavior of the WSS solution better than the conventional linear closure models.However, their model contains complicated procedures for analytical continuation of the plasma dispersion function with matrix arguments including nonlinear effects and it cannot be easily extended to systems with more than three modes.Our closure model does not contain such fundamental difficulties for application to systems with a large number of degrees of freedom.Another nonlinear closure model based on the phase velocity transform ͑PVT͒ was presented by Mattor. 16As written in Ref. 17, one difficulty with the original PVT closure is that calculating the transform into the phase velocity space can be nearly as time-consuming as the original kinetic equation.Therefore, Mattor proposed a Wentzel-Kramers-Brillouin ͑WKB͒-type local approximation to make the PVT closure model more tractable. 17,18In this local PVT model as well as in all other closure models except for the original PVT model and ours, it is assumed that, for the wave number vector k, the distribution function is written as f k (v ʈ ) ϫexp(Ϫi k tϩik•x) with a single dominant complex-valued frequency k .However, we show that, in the nonlinear stage with the time reversal symmetry, a complex-conjugate pair of functions f k (v ʈ ) and f k *(v ʈ ) corresponding to k and k * are required to describe the distribution function for each k, which is what our closure model attempts to describe.The collisionless turbulent system is also an interesting subject from the nonequilibrium thermodynamical point of view.In the collisionless system, an increase rate of the entropy functional defined in terms of the kinetic distribution function is equal to the product of transport fluxes and thermodynamic forces.Thus, no turbulent transport occurs in the collisionless steady state, or the entropy grows indefinitely in time for the transport fluxes to take finite nonzero values.This is called the entropy paradox by Krommes and Hu. 190][21] This thermodynamic aspect of the collisionless closure model is also discussed in the present work.
The rest of this work is organized as follows.In Sec.II, the new closure model is presented to give a closed set of fluid equations, which describe the collisionless system governed by the drift kinetic equation for ions, the adiabatic condition for electrons, and quasineutrality.When there exists a linearly unstable normal mode solution of the collisionless system, its complex conjugate also becomes a linear solution due to the time reversal symmetry.Our closure model is derived in order to express the parallel heat flux in terms of the temperature and the parallel flow for the normal mode as well as for the complex-conjugate mode.Also, the role of the closure model in the turbulence saturation and the anomalous transport is investigated based on the kinetic and fluid entropy balance equations.Application of the closure model to the three-mode ITG problem is given in Sec.III.The fluid system of equations with our model used can reproduce the exact nonlinear kinetic solution by WSS.It is also shown that oscillatory behaviors and initial amplitude dependence of other numerical kinetic solutions can be accurately described by the fluid system.Finally, conclusions are given in Sec.IV.

A. Basic equations
Let us start from the drift kinetic equation for ions with the mass m i and the charge e in a slab geometry where f ϭ f (x,v ʈ ,t) is the ion distribution function integrated over the velocity v Ќ perpendicular to the magnetic-field B ϵBb, v ʈ ϵv•b is the parallel velocity, E ʈ ϵE•b is the parallel electric field, and v E ϭcEϫB/B 2 is the EϫB drift velocity.The long perpendicular wavelength limit (k Ќ i Ӷ1) is assumed in Eq. ͑1͒.We assume that the magnetic-field B is uniform and constant.The electric field is written in terms of the electrostatic potential as EϭϪٌ.The distribution function is given by the sum of a background Maxwellian part F M ϭn 0 (m i /2T i ) 1/2 exp(Ϫm i v 2 /2T i ) and a perturbation part f ˜, where n 0 and T i denote the background density and temperature, respectively, which depend only on a coordinate x in the perpendicular direction.
We assume the adiabatic electron response.Then, the quasineutrality condition is written as Equations ͑1͒ and ͑2͒ form a closed system of equations to determine f ˜and .
Now, we represent f ˜and in terms of the Fourier expansion as

͑3͒
The drift kinetic equation ͑1͒ is rewritten in the wave number representation as where k ʈ ϵk•b is the parallel wave number and the parallel nonlinearity (e/m i )E ʈ ‫ץ‬ f ˜/‫ץ‬v ʈ in Eq. ͑1͒ is neglected.In Eq. ͑4͒, inhomogeneities in the background density and temperature are taken into account only through * i ϵ(cT i /eB)k Ќ •bϫٌ ln n 0 and i ϵd ln T i /d ln n 0 while n 0 and T i in other places as well as * i and i are regarded as constants.Taking the velocity moments of Eq. ͑4͒, we obtain fluid equations where Here all nonlinear terms result from the EϫB drift.From Eq. ͑2͒, we also have which makes the nonlinear terms in Eq. ͑5͒ vanish.A system of the fluid equations ͑5͒-͑7͒ and the quasineutrality condi-tion ͑8͒ are not closed because the parallel heat flow q k is included in Eq. ͑7͒.A new collisionless closure model for q k is given in Sec.II C. Before presenting the closure model, it is useful to review the linear kinetic solution of Eqs.͑4͒ and ͑8͒ in the next subsection.

B. Linear kinetic solution
Here, we consider the linear solution of Eqs.͑4͒ and ͑8͒ for a specified wave number vector k.The initial-value problem of the linearized version of Eqs.͑4͒ and ͑8͒ is easily solved by using the Laplace transform.The linear solutions f k (v ʈ ,t) and k (t) are written by where C is a contour in the complex -plane which lies above all of the singular points of the integrand.Here, f k (v ʈ ,) and k () are given by and respectively, where the dispersion function D k () is defined by and In Eqs.͑12͒ and ͑13͒, L is a Landau contour in the complex v ʈ -plane which pass below the singular point v ʈ ϭ/k ʈ for k ʈ Ͼ0 ͑above the singular point for k ʈ Ͻ0).
We consider the normal-mode solution which takes the form of Here, Lk is a constant ͑the normalization e Lk /T i ϭ1 is used later͒, Lk is determined by the dispersion relation D k ( Lk )ϭ0, and f Lk (v ʈ ) is given in terms of Lk as Here, the system is assumed to be linearly unstable and the most unstable mode which has the largest growth rate Im( Lk )(Ͼ0) is considered.We can directly confirm that the normal-mode solution in Eq. ͑14͒ satisfies the linearized version of Eqs.͑4͒ and ͑8͒.Also, as well known, asymptotic behaviors of general linear solutions are dominantly described by the normal-mode solution which has the largest growth rate.When we use the Fourier transform in time instead of the Laplace transform, we find that D 0k ( Lk )ϭ0 as well as D 0k ( Lk * )ϭ0 are satisfied for the linearly-unstable case Im( Lk )Ͼ0, where Lk * denotes the complex conjugate of Lk and D 0k is defined by Eq. ͑12͒ with the real v ʈ -axis used instead of the Landau contour L. Corresponding to this complex conjugate Lk * , let us take the complex conjugate of the normal-mode solution

͑16͒
Substituting Eq. ͑16͒ into Eq.͑13͒ and using Eqs.͑9͒-͑15͒, e Lk /T i ϭ1 and D k ( Lk )ϭD 0k ( Lk )ϭ0, we find that, for the initial condition f k (v ʈ ,tϭ0)ϭ f Lk * (v ʈ ) which show that the complex conjugate of the normal-mode solution in Eq. ( 16) is a linear solution of Eqs. ( 4) and (8) as well as the normal-mode solution.Since Im( Lk )Ͼ0 is assumed, the complex-conjugate solution exponentially damps in the course of time.The existence of the complexconjugate solution is related to the time reversal symmetry of the original kinetic equation and it should not be confused with what results from the reality condition in the Fourier representation.We should note that, for the case of Im( Lk )Ͻ0, the complex-conjugate of the normal-mode solution is not a linear solution of Eqs.͑4͒ and ͑8͒ as shown from D k ( Lk ) D 0k ( Lk ).

C. Closure model
Let us consider a closure problem for the parallel heat flow term q k in Eq. ͑7͒.Using Eqs.͑8͒ and ͑9͒ in Ref. 15 to give q k (nm) and T k (nm) , respectively, with the phase velocity W→ Lk /k ʈ , we find that, for the normal-mode solution given in the previous subsection, the parallel heat flow is related to the temperature perturbation as where

͑19͒
where B(w) is defined by 15 B͑w ͒ϵϪ 1 for k ʈ Ͼ0.Here, Z (n) () denotes the nth derivative of the plasma dispersion function Z() defined by for Im()Ͼ0 and by its analytic continuation for Im()Ͻ0.
For k ʈ Ͻ0, B(w) is determined by the relation B ϩ (w) ϭ͓B Ϫ (w*)͔* where B ϩ (w) and B Ϫ (w) denote B(w) for k ʈ Ͼ0 and B(w) for k ʈ Ͻ0, respectively, and * represents the complex conjugate.The closure relation given by Eqs.͑18͒-͑21͒, which is exactly valid for the linear normal-mode solution, is the one which is well reproduced by conventional linear kinetic-fluid closure models by Hammett and Perkins, 10 and by Chang and Callen. 11In the limit of Lk /k ʈ v t →0, we obtain C Lk ϭϪ2(2/) 1/2 ik ʈ /͉k ʈ ͉.This case corresponds to the result by Hammett and Perkins ͑hp͒, in which Eq. ͑18͒ leads to the dissipative parallel heat flow q k ϭϪn 0 ʈ hp ik ʈ T k with the parallel heat diffusivity ʈ hp ϭ2(2/) 1/2 v t /͉k ʈ ͉.Chang and Callen suggested in Eqs.͑54͒-͑61͒ of Ref. 11 that simplified versions of normal-mode closure relations, which are obtained by replacing Lk with i‫ץ‬ t under the conditions such as the adiabatic region and the one-pole approximation, should be used for numerical simulations.These conventional linear closure models cause the dissipation or time irreversibility for both stable and unstable modes.Hammett and Perkins 10 showed that there exists a non-Maxwellian equilibrium function f 0 (v ʈ ) for which their closure gives an exact linear response to the electrostatic potential and accordingly an exact dispersion relation for the normal mode.However, this does not imply that their closure exactly describes any linear solution of the initial value problem with an arbitrary initial perturbation f k (v ʈ ,tϭ0) added to that equilibrium.For example, the complex-valued frequency Lk * of the complex- conjugate solution derives not from the dispersion function D k () but from the initial condition term I k () ͓see Eq. ͑17͔͒.In this case of the complex-conjugate solution, their closure cannot correctly give the parallel heat flow q k as shown below.
It is easily found that, for the complex-conjugate solution in Eq. ͑16͒, the parallel heat flow is related to the temperature perturbation in terms of C Lk * ͑complex conjugate of where For the complexconjugate solution, the phase difference between q k and T k , which is an important element to consider by closure models, is opposite to that for the normal mode as seen from Eqs. ͑18͒ and ͑22͒.We should note that the conventional linear closure models can not describe the relation in Eq. ͑22͒ for the complex-conjugate solution and that this is the reason why they fail to reproduce the amplitude oscillations of the nonlinear solution of the three-mode ITG problem as shown in Sec.III.In order to keep the time reversal symmetry of the original kinetic equation, closure models need to be valid for both cases of the normal-mode solution and the complexconjugate solution.In our new collisionless fluid closure model, the parallel heat flux is related to the lower-order moments ͑i.e., density, parallel flow, and temperature pertur-bations͒ such that this relation is valid both for the normalmode solution and the complex-conjugate solution as well as for any linear combination of these solutions.Let us express such a closure relation by q k ϭF k (T k ,u k ,n k ).For distribution functions given by linear combinations of the normalmode solution and the complex-conjugate solution the density, parallel flow, temperature, and parallel heat flow are written as where ␣ k Ͼ and ␣ k Ͻ are arbitrary constants.Substituting these expressions into the closure relation, we obtain Therefore, we write the closure relation Furthermore, from comparison between the cases of (␣ k Ͼ ,␣ k Ͻ )ϭ(1,0) and (0,1) in Eq. ͑24͒, we find that F k (T Lk * ,u Lk * ,n Lk * )ϭ͓F k (T Lk ,u Lk ,n Lk )͔* and, therefore, the dimensionless coefficients C Tk , C uk , and C nk in Eq. ͑25͒ are real-valued constants determined so as to satisfy Equation ͑26͒ is required from Eqs. ͑18͒ and ͑25͒.Equation ͑25͒ is called the nondissipative closure model ͑NCM͒ hereafter.Here, C Lk is given by Eq. ͑19͒, and where * e ϵϪ(T e /T i ) * i .We obtain from Eqs. ͑26͒ and , and k ϵ rk ϩi ik ϵ Lk ( Lk Ϫ * e )/(k ʈ 2 v t 2 )Ϫ1 ϪT e /T i are used.For arbitrarily given C nk , the closure relation in Eq. ͑25͒ with Eq. ͑28͒ satisfies the requirement described before Eq.͑23͒.In order to uniquely determine (C Tk ,C uk ,C nk ), we impose an additional condition that Eq. ͑25͒ should be valid when the perturbation distribution function takes a Maxwellian form, f k ϰexp(Ϫm i v 2 /2T i ), for which u k ϭT k ϭq k ϭ0.Then, we obtain from which with Eq. ͑28͒, C Tk and C uk are also uniquely determined.
So far, we have considered only linearly unstable modes ͓Im( Lk )Ͼ0͔.As mentioned at the end of Sec.II B, the complex-conjugate of the normal mode is not a linear solution for Im( Lk )Ͻ0.Thus, for linearly stable modes, we assume the normal mode to be a dominant contribution to the perturbation, and use the same closure relation as in Eq. ͑18͒ where C Lk is defined by Eq. ͑19͒.It may appear that different closure relations for unstable and stable modes give discontinuity in the parallel heat flux q k across the point of marginal stability ͓Im( Lk ϭ0)͔.However, as seen from Eq. ͑26͒, we find that, even if apparent forms are different, the NCM in Eq. ͑25͒ give the same q k as Eq.͑30͒ when the fluid variables (T k ,u k ,n k ) are described by the unstable normal mode, (T Lk ,u Lk ,n Lk ).Thus, the discontinuity in q k across the marginal point does not occur if the fluid variables change continuously and are dominated by the normal mode near the marginal point.

D. Entropy balances
It is useful to consider entropy balances in nonequilibrium systems for understanding the relation between entropy production, dissipation, thermodynamic forces, and transport fluxes.Here, kinetic and fluid entropy balances are examined for the turbulent system governed by the equations in Sec.
II A in order to elucidate the role of the closure model for the turbulence saturation and the anomalous transport.
We define the microscopic entropy S m and macroscopic entropy S M for ions per unit volume by 20

͑32͒
where f ϭF M ϩ f ˜.Retaining terms up to O( f ˜2), the relation between S m and S M is given by where ͗¯͘ represents the ensemble average.As shown from the basic kinetic equation ͑1͒ ͓or Eq. ͑4͔͒, the total microscopic entropy ͐d  19 ͒ as the entropy associated with the fluctuations because, in the collisionless system, turbulent processes produce not S m but S M .From Eq. ͑4͒, we find that the turbulent entropy production rate is given by where v Ek ϵi(c/B)bϫk k is the turbulent EϫB velocity, q Ќ ϵ 1 2 n 0 ͚ k Re(T k v Ek * ) is the turbulent perpendicular heat flux, E ʈ k ϵϪik ʈ k is the parallel electric field, and Q ϵen 0 ͚ k Re(E ʈk u k *) represents the turbulent ion heating.

͑37͒
where where the orthogonality conditions for the Hermite polynomials (2) Ϫ1/2 ͐ Ϫϱ ϱ dx e Ϫx 2 /2 H n (x)H m (x)ϭn!␦ nm (n,m ϭ0,1,2,...) are used.It is found from Eq. ͑36͒ or Eq.͑38͒ that there is no perpendicular heat transport in the direction of the temperature gradient in the collisionless and steady turbulent state where saturation of the entropy ͐dv ʈ ͉ f k ͉ 2 /2F M and the potential amplitude ͉ k ͉ occurs.͓If the collision term and the parallel-nonlinearity term are retained in Eq. ͑4͒, there appear corresponding terms which balance the heat transport term in the right-hand side of Eqs.͑36͒ and ͑38͒, although they are typically small. 2͔ Then, in order to treat the anomalous transport in the collisionless turbulent plasma, we assume a quasisteady state in which the amplitudes of the fluid variables nk with low n (n k ,u k ,T k , q k ,...) and the potential k reach the steady state while the high-n moments included in ͚ nу4 (n!/2) ͉ nk ͉ 2 grow indefinitely in time.The quasisteady state should be regarded as idealization of the real steady state in which those high-n moments eventually saturate as well due to collisional dissipation even if the collision frequency is much smaller than the characteristic frequency of the instabilities causing the turbulence.We should note that the moments with higher n are more effectively dissipated by the Fokker-Planck-type collision operator because they represents finer structure with steep gradients in the velocity space.
Using Eqs.͑5͒-͑8͒, we can derive another balance equation similar to Eq. ͑38͒

͑39͒
which can be regarded as the fluid entropy balance equation.
Comparing Eq. ͑39͒ to Eq. ͑38͒, we easily find

͑40͒
For the quasisteady state in which d( ͚ nр3 (n!/ 2) ͉ nk ͉ 2 )/dtϭ0, we obtain from Eqs. ͑39͒ and ͑40͒, Equation ͑41͒ represents the entropy production rate in the quasisteady state, where the perpendicular heat transport in the presence of the background temperature gradient drives the growth of the fluctuations in the high-n moment variables through the correlation between the parallel heat flux and temperature fluctuations.Thus, it is considered that the entropy production rate and the anomalous perpendicular heat transport are deeply dependent on what closure model is used for the parallel heat flux.If we use the Hammett-Perkins model q k ϭϪn 0 ʈ hp ik ʈ T k with ʈ hp ϭ2(2/ ) 1/2 v t /͉k ʈ ͉, the entropy production rate is given by ͑ for the Hammett-Perkins model͒, ͑42͒ which takes the form of the dissipation for all wave number vectors k.On the other hand, when we use the NCM in Eqs.͑25͒-͑29͒ and Eq.͑30͒, we obtain ͑ for linearly stable modes͒ ͑43͒ and ͑ for linearly unstable modes͒, ͑44͒ where Eqs.͑5͒, ͑6͒, and ͑8͒ are also used.Then, from Eqs. ͑43͒ and ͑44͒, the entropy production rate for the quasisteady state is given by where ͚ k(stable) ( ͚ k(unstable) ) represents the summation over k corresponding to linearly stable ͑linearly unstable͒ modes.
Here, we find that and that C uk Ͼ0 for linearly unstable modes.Compared to Eq. ͑42͒ for the Hammett-Perkins model, Eq. ͑45͒ shows that, in our model, the entropy production rate consists of the dissipation terms for the linearly stable modes ͓the first group of terms in the right-hand side of Eq. ͑45͔͒ similar to Eq. ͑42͒ and the nondissipative terms ͑the second group of terms͒: The latter terms result from the second term in the right-hand side of Eq. ͑25͒ and are associated with the transfer of ͉u k 2 ͉ from the unstable modes to stable modes in the k-space.Due to this nondissipative nature for unstable modes, turbulence simulations using the NCM may give different results on the entropy production rate and the anomalous heat transport from those obtained by using the conventional closure models.In the next section, an example of application of the NCM is shown.

III. THREE-MODE ITG SYSTEM
Here, we apply the collisionless fluid closure given in the previous section to the three-mode ITG problem. 2,14

A. Three-mode ITG equations and exact nonlinear solution
We consider a rectangular domain of L x ϫL y in the x-y plane with a uniform external magnetic-field BϭB 0 (z ˆϩy ˆ) (͉͉Ӷ1), where y ˆand z ˆdenote the unit vectors in the y and z directions, respectively.The system is assumed to be homogeneous in the z direction (‫ץ/ץ‬zϭ0).The background density and temperature gradients are assumed to exist in the x direction, and their gradient scale lengths are given by L n ϭϪ(d ln n 0 /dx) Ϫ1 (Ͼ0) and L T ϭϪ(d ln T i /dx) Ϫ1 (Ͼ0), respectively.We employ the periodic boundary conditions in both x and y directions.Then, in Eq. ͑3͒, we can write k ϭ2͓(m/L x )x ˆϩ(n/L y )y ˆ͔, f k ϭ f m,n , and k ϭ m,n (m,n ϭ0,Ϯ1,Ϯ2,...).In the three-mode ITG system, we only keep (m,n)ϭ(Ϯ1,Ϯ1) and (Ϯ2,0) modes with the symmetry conditions of and Re(f 2,0 )ϭ0.These symmetry conditions imply that the system has boundaries at xϭ(2lϩ1)L x /4 (lϭ0,Ϯ1, Ϯ2,...), where f and vanish.Here, f 20 describes a quasilinear flattening of the temperature profile around xϭ0, which turns off the linear instability drive for the ͑1,1͒ mode.
Other modes with higher wave numbers, which may be nonlinearly destabilized by f 20 steepening the temperature gradient near the boundaries, are not included.Since effects of finite gyroradii are neglected, we have 20 ϭ0 and ignore the nonlinear generation of zonal flows, which can play an important role in ITG turbulence. 5,13,22,23ereafter, f 1,1 , Im(f 2,0 ) and 1,1 are, respectively, denoted by f 1 , h, and 1 for simplicity, where f 1 and 1 are complex-valued while h is real-valued.Then, from Eqs. ͑4͒ and ͑8͒, the three-mode ITG equations are written as 2,14 where L x ϭL y ϭ1/k and T i ϭT e (T e : the electron tempera-ture͒ are assumed and G(v ʈ ) is defined by Here, we have used dimensionless normalized variables x ϭxЈ/ i , yϭyЈ/ i , vϭvЈ/v t , tϭtЈv t /L n , f ϭ f ЈL n v t / i n 0 , and ϭeЈL n /T i i , where prime represents a dimensional quantity, v t ϭͱT i /m i is the ion thermal velocity, i ϭv t /⍀ i is the ion thermal gyroradius, and ⍀ i ϭeB/m i c is the ion gyrofrequency.Two important parameters ⌰ and i in Eq. ͑49͒ are given by ⌰ϭL n / i and i ϭL n /L T , respectively.From Eqs. ͑46͒-͑48͒, the entropy balance equation similar to Eq. ͑36͒ is derived as where q x ϵRe( 1 2 T 1 ik 1 *) with T 1 ϵ͐dv ʈ f 1 (v ʈ 2 Ϫ1) is the heat flux in the x direction, and F M ϵe Ϫv ʈ 2 /2 /(2) 1/2 .
An exact solution of the three-mode ITG equations found by WSS is written as where f Lr (v ʈ ) and f Li (v ʈ ) represent the real and imaginary parts of the eigenfunction Then, from Eqs. ͑51͒ and ͑53͒, we obtain Therefore, this exact nonlinear solution is given by the sum of the linear normal-mode solution with the amplitude ␣ Ͼ (t) and the complex conjugate solution with the amplitude ␣ Ͻ (t).The ordinary differential equations ͑52͒ are rewritten as which show that the normal-mode and complex-conjugate solutions interact with each other through the nonlinearcoupling terms proportional to c. Since the exact nonlinear solution in Eq. ͑54͒ is represented as a linear combination of the normal-mode and complex-conjugate solutions, it is accurately described by the the NCM presented in Sec.II C.
Figure 1 shows ␣ Ͼ (t) and ␣ Ͻ (t) obtained from the exact solution for kϭ0.1, ⌰ϭ1, i ϭ10, ␣ Ͼ (0)ϭ␣ Ͻ (0) ϭ0.005 ͓a(0)ϭ0.01,b(0)ϭ0͔and c(0)ϭ0.For this case, the linear eigenfrequency is given by Lk ϭϪ0.1149 ϩ0.0831i and the period of ͓␣ Ͼ (t),␣ Ͻ (t),c(t)͔ is T ϭ186.78.We find that peaks of the complex-conjugatemode amplitude ␣ Ͻ appear after those of the normal-mode amplitude ␣ Ͼ .Peaks of the potential amplitude ͉(t)͉ ϭa(t)ϭ␣ Ͼ (t)ϩ␣ Ͻ (t) occur between the peaks of ␣ Ͻ and ␣ Ͼ at tϭ(nϩ1/2)T (nϭ0,1,2,...) where ␣ Ͼ becomes equal to ␣ Ͻ .We should note that, in the nonlinear regimes around the potential's peaks, both the normal mode and its complex conjugate have large amplitudes comparable to each other so that the nonlinear behavior cannot be well described by the conventional linear closure models which do not take account of the complex-conjugate mode.The nonlinear closure model proposed by Mattor and Parker succeeded in giving the periodic potential amplitude oscillation in the three-mode ITG problem. 15Thus, the normal mode and its complex conjugate seem to be included in their model.However, since their nonlinear closure model as well as the local PVT model by Mattor 17 assume that the kinetic distribution function f k is written by a single mode structure corresponding to a dominant complex-valued frequency k , they do not accurately describe the nonlinear states which consist of distinct mode structures with comparable amplitudes as found here.In the next subsection, the closed set of fluid equations for the three-mode ITG system, which are obtained by applying the collisionless fluid closure, are examined in detail.

B. Fluid equations for the three-mode ITG system
Taking the velocity moments of Eqs.͑46͒ and ͑47͒, we obtain where ͓n 1 (t), u 1 (t), T 1 (t),q 1 (t We also obtain from Eq. ͑48͒ Here, we should note that h (t)ϭn h (t)ϵ͐ Ϫϱ ϱ dv ʈ h(v ʈ ,t) ϭ0.The entropy balance equation similar to Eq. ͑39͒ is derived from Eqs. ͑56͒-͑61͒ as

͑62͒
where q x ϵRe( 1 2 T 1 ik 1 *) is the heat flux in the x direction.
Comparing Eq. ͑62͒ to Eq. ͑50͒ and using the Hermitepolynomial expansion similar to Eq. ͑37͒, we obtain where ¯represents the high-order-moment terms corresponding to ͚ nу4 (n!/2) ͉ nk ͉ 2 in Eq. ͑40͒.For the case of the exact solution shown in Sec.III A, the time average of the left-hand side of Eq. ͑62͒ and that of the right-hand side of Eq. ͑63͒ vanish and accordingly the time average of the heat flux q x vanishes.Thus, it is not practical to consider the three-mode ITG system for estimating the anomalous heat diffusivity in the steady turbulent state, although it is a useful nonlinear system for examining the validity of kineticfluid closure models.There exists another simple set of fluid equations, which is used to evaluate the anomalous heat diffusivity due to the toroidal ITG mode. 24ow, let us apply the NCM presented in Eq. ͑25͒ and put q 1 ϭC T1 T 1 ϩC u1 u 1 ϩC n1 n 1 , ͑64͒ where C T1 , C u1 , and C n1 are real coefficients given from Eqs. ͑28͒ and ͑29͒ as Im͑ 1 1 *͒ Equations ͑56͒-͑61͒ with the NCM given by Eqs.͑64͒-͑67͒ form a closed system of fluid equations for the threemode ITG system.Figures 2͑a͒ and 2͑b͒ show ͉ 1 (t)͉ and Re͓ 1 (t)/͉ 1 (t)͉͔, respectively, which are obtained by numerically solving these fluid equations.A nondissipative time integration scheme 14 is employed for these numerical calculations.Here, kϭ0.1, ⌰ϭ1, and i ϭ10 are used.For this case, we have C L1 ϭϪ0.381Ϫ1.186i,C T1 ϭϪ2.878ϫ10Ϫ2 and C u1 ϭ1.477 with C n1 ϭ0.The initial conditions are n 1 (0)ϭ 1 (0)ϭ0.01,u 1 (0)ϭϪ2.149ϫn 1 (0), T 1 (0) ϭϪ0.220ϫn 1 (0), and u h (0)ϭt h (0)ϭ0, which are consistent with the initial conditions ͓a(0),b(0),c(0)͔ ϭ(0.01,0,0) for the WSS solution in Eq. ͑51͒ ͑the same initial conditions as used in Fig. 1͒.The fluid equations in Eqs.͑56͒-͑61͒ can reproduce the WSS exact nonlinear solution since the NCM given by Eqs.͑64͒-͑67͒ is completely satisfied by the WSS solution.It is noted that the WSS solution can be reproduced by the closure relation given by Eqs.͑64͒-͑66͒ even with an arbitrary nonzero real value given to C n1 as shown in Appendix B. Also, it can be shown in Appendix B that, even if Eqs. ͑65͒-͑67͒ are not satisfied, the closure relation in Eq. ͑64͒ with real-valued coefficients C T1 , C u1 , and C n1 can give the WSS-type double-periodic nonlinear solution, which consists of the normal mode and its complex conjugate, since the time reversibility is still retained.For example, the case of the most simple closure, in which q 1 ϭ0 (C T1 ϭC u1 ϭC n1 ϭ0), and that of the Hammett-Perkins closure q 1 ϭϪ2(2/) 1/2 ik⌰T 1 are plotted for comparison to the WSS solution in Fig. 2. The time evolution of the potential for the q 1 ϭ0 case shows a qualitatively better agreement with the exact solution than for the Hammett-Perkins model.In the Hammett-Perkins model, the potential ͉ 1 (t)͉ is saturated at a certain amplitude, which is in contrast to the periodic amplitude oscillation shown by the exact solution.Also, in the Hammett-Perkins model, the time average of the heat flux q x in the three-mode ITG system takes on an artificial nonzero value since i q ¯x ϭϪRe( 1 2 T 1 ik⌰q 1 *)ϭ(2/) 1/2 (k⌰) 2 ͉T 1 ͉ 2 Ͼ0.This q ¯x follows from the entropy balance equation ͑62͒ with q 1 ϭϪ2(2/) 1/2 ik⌰T 1 where ¯represents the time average.
Next, we compare solutions of the fluid equations with kinetic solutions for the case in which initial conditions are inconsistent with the WSS solution.The simplest example of such initial conditions is given by the Maxwellian form of the initial distribution function f 1 (v ʈ ,tϭ0)ϰe Ϫv ʈ 2 /2 / (2) 1/2 .Figures 3͑a͒ and 3͑b͒ show ͉ 1 (t)͉ and Re͓ 1 (t)/͉ 1 (t)͉͔, respectively, which are obtained by nu- The initial conditions are given by n 1 (0)ϭ 1 (0)ϭ0.01,u 1 (0)ϭϪ2.149ϫn 1 (0), T 1 (0)ϭϪ0.220ϫn 1 (0), and u h (0)ϭt h (0)ϭ0, which are consistent with the initial conditions used for the WSS exact kinetic solution in Fig. 1.Solid curves in ͑a͒ and ͑b͒ represent ͉ 1 (t)͉ and Re͓ 1 (t)/͉ 1 (t)͉͔, respectively, obtained by using the NCM in Eqs.͑64͒-͑67͒ and completely agree with those of the WSS exact kinetic solution.Results obtained by using the Hammett-Perkins model and the q 1 ϭ0 model are also shown by dashed and dotted curves, respectively.
14. Results obtained by numerically solving Eqs.͑56͒-͑61͒ and ͑64͒-͑67͒ for the initial conditions n 1 (0)ϭ 1 (0) ϭ0.01 and u 1 (0)ϭT 1 (0)ϭu h (0)ϭt h (0)ϭ0 are also plotted in Figs.3͑a͒ and 3͑b͒ for comparison to the kinetic solution.Kinetic and fluid results for this case are in a good agreement except that a small potential phase deviation appears after the potential amplitude recurs to the low level of the initial condition.

IV. CONCLUSIONS
In the present paper, we have presented a new nondissipative closure model ͑NCM͒.We have shown that both linearly unstable modes and their complex conjugate are exact solutions of the linearized version of the collisionless kinetic equation with time reversibility.In order to take account of the nondissipative nature of the kinetic equation, the NCM is derived such that the relation between the fluctuations in the parallel heat flux, the temperature, and the parallel flow is valid both for linearly unstable normal modes and their complex conjugate as well as for any linear combination of them.The fluctuations obtained by the complex conjugate of the linearly stable modes are no longer linear solutions.For the linearly stable modes, we use the dissipative closure relation, which reduces to previous closure models by Hammett and Perkins, 10 and by Chang and Callen 11 in limiting cases.
Kinetic and fluid entropy balances for the collisionless system are investigated to compare effects of the NCM and the Hammett-Perkins model on the entropy production, the turbulence saturation, and the anomalous perpendicular heat transport.The quasisteady state for the collisionless turbulent plasma is defined as the state in which the amplitudes of fluid variables corresponding to low-order velocity moments of the distribution function saturate but the high-order moments FIG. 3. Comparison between kinetic and fluid solutions.Solid curves in ͑a͒ and ͑b͒ represent ͉ 1 (t)͉ and Re͓ 1 (t)/͉ 1 (t)͉͔, respectively, obtained by numerically solving the kinetic system of Eqs.͑46͒-͑48͒ for kϭ0.1, ⌰ ϭ1, and i ϭ10 with the initial conditions f 1 (v ʈ ,0)ϭ0.01ϫeϪv ʈ 2 /2 /(2) 1/2 and h(v ʈ ,0)ϭ0.Dashed curves represent results obtained by numerically solving the fluid system of Eqs.͑56͒-͑61͒ with the NCM in Eqs.͑64͒-͑67͒ for the initial conditions n 1 (0)ϭ 1 (0)ϭ0.01and u 1 (0)ϭT 1 (0)ϭu h (0) ϭt h (0)ϭ0.FIG. 4. Period T of the potential amplitude ͉ 1 (t)͉ as a function of the initial amplitude ͉ 1 (0)͉.Here, kϭ0.1, ⌰ϭ1, i ϭ10, f 1 (v ʈ ,0) ϭn 1 (0)e Ϫv ʈ 2 /2 /(2) 1/2 ͓n 1 (0)ϭ 1 (0)͔, and h(v ʈ ,0)ϭ0.are used.Solid and open circles represent results from solving the kinetic system of Eqs.͑46͒-͑48͒ and those from solving the fluid system of Eqs.͑56͒-͑61͒ with the NCM in Eqs.͑64͒-͑67͒, respectively.and accordingly the total entropy grow indefinitely in time.The entropy production rate given by the product of the perpendicular hear flux and the temperature gradient balance with the correlation between the parallel heat flux and temperature fluctuations in the quasisteady state.In the NCM, this correlation function takes the nondissipative form in the unstable wave number region while, in the Hammett-Perkins model, it is written as a dissipation proportional to the square of the temperature fluctuation amplitude for all wave number regions.Then, it is expected that, when used in turbulence simulations, the NCM predicts different results on the entropy production rate and the anomalous heat transport from those given by the Hammett-Perkins and other dissipative closure models.
The NCM is applied to the three-mode ITG problem as an example.The WSS exact nonlinear kinetic solution can be reproduced from the fluid system of equations with the NCM.Other numerical kinetic solutions can also be accurately described by the fluid system.Especially the oscillatory behaviors of the potential amplitude and the dependence of the oscillation period on the initial amplitude, which cannot be treated by the conventional dissipative closure models, are well reproduced by the NCM.While the three-mode ITG problem is an interesting and successful nonlinear test case, one might wonder whether the NCM also works well for more realistic many-mode turbulent cases which significantly differ from the three-mode case in some aspects.The three-mode system involves interaction between only a single coherent eddy of the ͑1,1͒-mode and a background profile flattening of the ͑2,0͒-mode, and shows fairly regular periodic behavior as given by the exact solution in Eq. ͑51͒.However, the behavior of nonintegrable chaotic systems such as many-mode turbulent systems is very different from such a regular one, and the validity of the use of a nondissipative or a dissipative closure depending on the linear growth rate criterion should be checked for such fully turbulent cases.͑Note that even the fluid system using linear closure models can be highly nonlinear and turbulent because of the EϫB convection terms although it does not necessarily imply the correct description of the original turbulent kinetic system by the fluid model.͒Another point to clarify is collisional effects on the closure.Collisions are necessary for true steady-state turbulence of kinetic systems to be realized, and even a small amount of collisions may have significant effects.In order to investigate the subjects discussed above, applications of the NCM to systems with a higher number of degrees of freedom is work currently in progress.
Here, the complex eigenfrequency L ϭ Lr ϩi␥ is given from Eqs. ͑B2͒ as a solution of the cubic algebraic equation with ␥Ͼ0.Substituting Eqs.͑B1͒ into Eqs.͑56͒-͑60͒, we find that a(t), b(t), and c(t) should satisfy da/dtϭ␥b, db/dtϭ␥aϪ2k2 ac, ͑B3͒ dc/dtϭ4k 2 ab, which are the same as Eqs.͑52͒ for the WSS kinetic solution.
From Eqs. ͑53͒ and ͑B3͒, the same equations as Eqs.͑55͒ for ␣ Ͼ (t), ␣ Ͻ (t), and c(t) are also derived.Now, it turns out that a WSS-type double-periodic solution given by Eqs.
͑B1͒-͑B3͒ appears from this system of fluid equations.Furthermore, when imposing the conditions in Eqs.͑65͒ and ͑66͒ on the real coefficients C T1 , C u1 , and C n1 ͓Eq.͑67͒ is not necessary͔, the complex eigenfrequency L ϭ Lr ϩi␥ also coincides with that given by the kinetic dispersion relation.Then, the system of fluid equations includes the WSS solution as one of their solutions.In other words, the exact kinetic solution of the three-mode ITG problem can be reproduced by the fluid equations with a proper collisionless closure model.
3xS m is conserved without collisions, although the macroscopic entropy S M or ⌬Sϵ 1 ϪS M in our notation͒ the entropy of the system to measure deviation of the fluctuating distribution function f from the ensemble-averaged distribution function ͗ f ͘.Here, we regard ⌬SϭS M ϪS m (ϭϪS ¯ϭF ¯inKrommes and Hu Eq. ͑15͒, respectively, and Lr is the real part of the eigenfrequency Lk where ␥ϵIm( Lk ).These equations are also rewritten in the Hamiltonian form as shown in Appendix A. Functional forms of the solutions a(t), b(t), and c(t) are given in terms of the Jacobi elliptic functions and they are found in Ref.14.We should note that a(t), b(t), and c(t) are periodic func-tions with the period T given by Eq. ͑32͒ in Ref. 14 so that the WSS solution in Eq. ͑51͒ is double-periodic in t with two periods T and 2/ Lr .Let us define ␣ Ͼ (t) and ␣ Ͻ (t) by Ϫ͐dv ʈ kG(v ʈ )/( L1 Ϫk⌰v ʈ )ϭ0 and Im( L1 )Ͼ0 is assumed.