Kinetic simulation of steady states of ion temperature gradient driven turbulence with weak collisionality

Statistically steady states of the ion temperature gradient driven turbulence with weak collisionality, where the collision frequency is much lower than characteristic ones of the turbulence, are investigated by means of a Eulerian kinetic simulation with high resolution. In the saturated state of the entropy variable, the ion heat transport balances with the collisional dissipation that is indispensable to realizing a steady-turbulence state of perturbed distribution function d f . The kinetic simulation deﬁnitely conﬁrms the conventional hypothesis that, in a low-collisionality limit, the low-order velocity-space moments of d f as well as the ion heat transport ﬂux agree with those in the quasisteady state of the collisionless turbulence with the constant entropy production. A spectral analysis of d f in the velocity-space clariﬁes the transfer and dissipation processes of the entropy variable associated with ﬂuctuations, where the phase mixing, the E 3 B nonlinearity, and the ﬁnite collisionality are taken into account. A power-law scaling predicted by the theoretical analysis is also veriﬁed by the simulations in a subrange of the power spectrum which is free from the entropy production and the collisional dissipation. © 2004 American Institute of Physics. @ DOI:


I. INTRODUCTION
Turbulent transport in high-temperature plasmas has long been a key issue in the magnetic confinement fusion research, 1 since it is considered as a main cause of the anomalous transport of particles and energy. Understanding the microscopic turbulence is important as the first step to prediction and control of the anomalous transport. Extensive simulation studies 2 on drift wave turbulence, such as the ion temperature gradient ͑ITG͒ mode, have revealed several important aspects of the turbulent transport in magnetically confined plasmas, for example, the transport suppression by the self-generated zonal flow. 3 Nevertheless, saturation mechanism of the collisionless turbulence has been an open question. Since the collisionless gyrokinetic equation has time-reversal symmetry, one needs to consider a coarsegrained form of the one-body velocity distribution function f with small-scale fluctuations in order to define an irreversible transport process in collisionless turbulence.
It has been pointed out that, when a steady transport flux is observed in the collisionless turbulence, a quasisteady state should be realized, 4 -6 where high-order velocity-space moments of the perturbed distribution function ␦ f continue to grow but the low-order ones are constant in average. Here, ␦ f ϵ f ϪF M is deviation from the equilibrium given by the Maxwellian velocity distribution F M . Existence of the quasisteady state in the collisionless ITG turbulence has been confirmed by means of a Eulerian ͑so-called Vlasov͒ numerical simulation of the gyrokinetic equation. 7 The phase mixing generates fine-scale fluctuations of ␦ f , and leads to continuous growth of the high-order moments, as well as an entropy variable associated with fluctuations ͑that is defined by a square integral of ␦ f ͒ of which the growth rate balances with the transport flux multiplied by a normalized ion temperature gradient ͑see Sec. III A for more detail͒.
If the collisionless assumption has a practical meaning for considering the steady anomalous transport in an actual plasma with weak but finite collisionality, the quasisteady state should be an idealization of a real statistically steady state in a weak-collisionality limit. Namely, for sufficiently low collision frequencies, statistical behaviors of the loworder moments such as the transport flux should agree with those in the collisionless turbulence. Our concern here is, thus, to find the collisionality dependence of the kinetic ITG turbulent transport, and to examine the above conjecture on the relation between the steady and quasisteady states. Introduction of the finite collisionality allows the system to approach the real steady state. Even if the collision frequency is much lower than characteristic ones of the ITG modes, it definitely affects evolution of the system through dissipation of the fine-scale fluctuations of ␦ f in the velocity space.
In simulations of the ITG turbulence shown below, we employ a two-dimensional slab model without complication of interaction between the zonal flow and the turbulence, as has been done in our previous work. 7 The Eulerian kinetic simulation with the simplified model setting enables us not only to investigate fundamental processes in plasma turbulence, such as the EϫB advection, the phase mixing, and the collisional dissipation, but also to give a useful reference for construction of kinetic-fluid closure models. 6,8 This is because collisionless fluid simulations of the steady turbulence transport are based on the conjecture on existence of the quasisteady state of turbulence. A detailed comparison between the collisionless kinetic and fluid simulations of the slab ITG turbulence has recently been carried out, 9 where the transport coefficient given by the fluid simulation with the nondissipative closure model is in good agreement with the kinetic results.
In the latter part of this paper, a velocity-space power spectrum of ␦ f represented by a quadratic form of the Hermite-polynomial expansion coefficients is investigated in analogy with the passive scalar convection in the homogeneous isotropic turbulence of a neutral fluid. 10 The spectral analysis elucidates the entropy transfer process from macroto microscales in the phase space through the phase mixing and the EϫB nonlinearity. The entropy variable is damped by collisions in a microvelocity scale. Similarly to the viscous-convective subrange in a power spectrum of the passive scalar, we identify a subrange in the power spectrum of ␦ f which is free from the entropy production and the collisional dissipation. A scaling law for the spectrum of the entropy variable in this subrange will also be derived in this paper. This paper is organized as follows. Our simulation model is described in Sec. II. Simulation results are presented in Sec. III, where the statistically steady state of the weakly collisional turbulence is discussed in Sec. III A. The collision frequency dependence of the transport flux and a convergence check to mode truncation in the wave number space are given in Secs. III B and III C, respectively. The spectral analysis of the distribution function is described in Sec. IV. The obtained results are summarized in Sec. V.

II. SIMULATION MODEL
Our simulation model considered here is the same as in the previous work 7 on the collisionless ITG turbulence except for an ion-ion collision term. We consider a periodic two-dimensional slab configuration with translational symmetry in the z direction, where the uniform magnetic field is set in the y -z plane such that BϭB(ẑ ϩŷ ) with Ӷ1. Equations numerically solved here are derived from a v Ќ integral of the gyrokinetic equations 11 by neglecting the parallel nonlinear term and by assuming where F M denotes the Maxwellian velocity distribution. We also assume constant density and temperature gradients of the background ions in the x direction with much larger scale lengths ͓L n ϵϪd(ln n)/dx and L T ϵϪd(ln T i )/dx] than the fluctuation wave lengths. This assumption enables us to simulate the local turbulent transport process without the quasilinear flattening of the background temperature profile, and is considered to be relevant to simulating realistic experimental conditions where the background inhomogeneities tend to persist as quasistatic profiles because of the continuous supply of particles and energy. 12 Therefore, we arrive at the following equations represented in the wave number space kϭ(k x ,k y ) as

͑2͒
where the electric potential k is related to ⌿ k by ⌿ k ϭe Ϫk 2 /2 k , with k 2 ϭk x 2 ϩk y 2 . The background electron temperature T e ϭT i and the adiabatic electron response are also assumed, such that ñ e,k ϭ k for k y 0. Here, for comparison between the collisionless and weakly collisional turbulence, we consider a limiting case with no zonal flow component of k y ϭ0 by fixing f k y ϭ0 ϭ k y ϭ0 ϭ0, as has been done in our previous simulation 7 with the aim of simulating a large transport level observed in a toroidal geometry. In the slab configuration, otherwise, the turbulence is too severely suppressed by the zonal flow to cause a mean transport. 7 The assumption of no zonal flow also enables us to examine a finite-collisionality effect on turbulence without complication of the collisional damping of the zonal flow and its interplay with turbulence. 13 In addition, we neglect k x ϭ0 modes of f k and k , since they are included in the background part with constant density and temperature gradients in the x direction. 12 Equations ͑1͒ and ͑2͒ are normalized as follows: and T i are the ion thermal velocity, the ion thermal gyroradius, the ion cyclotron frequency, the background plasma density, the elementary charge, and the background ion temperature (T i ϭm i v ti 2 ; m i means the ion mass͒, respectively. Prime means a dimensional quantity. ⌰ is defined as ⌰ ϭL n / i . i and ⌫ 0 (k 2 ) are given by i ϭL n /L T and ⌫ 0 (k 2 )ϭexp(Ϫk 2 )I 0 (k 2 ), respectively. I 0 (z) is the zeroth modified Bessel function of z.
The parallel advection term on the left-hand side of Eq. ͑1͒ contributes to generation of fine-scale fluctuations of f k in the velocity space, that is, the phase mixing. The instability drive is contained in the first term on the right-hand side of Eq. ͑1͒. The second term on the right-hand side denotes the ion-ion collision term for which we employ the Lenard-Bernstein model collision operator with the collision frequency normalized by v ti /L n . The collision operator in Eq. ͑3͒ makes f k approach F M preserving the mass. Although the momentum and the energy are not conserved, it doesn't cause a significant influence on the results in the present study. Velocity-space derivatives in the collision term are calculated in the velocity wave number space l into which f k (v ʈ ) is Fourier transformed from the v ʈ space discretized by a uniform grid in a range of Ϫv max рv ʈ рv max , with v max ϭ10v ti . Then, they are transformed back to the v ʈ space. f k is fixed to zero at v ʈ ϭϮv max , since the fluctuation amplitude at the velocity-space boundary is negligibly small. For the collisionless case, we set v max ϭ5v ti . We have employed sufficient resolution of the velocity space in accordance with the magnitude of ; finer grid spacing for v ʈ is necessary for a smaller value of . Numerical time integration is carried out by the fourth-order Runge-Kutta-Gill method, with careful convergence checks to the time step ⌬t so as to keep enough accuracy, while a nondissipative timeintegration scheme is employed for the collisionless simulation. 14,15 The minimum and maximum values of the wave number are set to k min ϭ0.1 and k max ϭ3.2, respectively, for both of the k x -and k y directions with the 3/2 rule for dealiasing in the spectral method. Results of a convergence check to k max are given in Sec. III C.

A. Steady state of weakly collisional turbulence
In the system described above, we note a balance equation of entropy variable defined by a functional, ␦S which is derived from Eq. ͑1͒ by multiplying f k */F M ͑where the asterisk denotes complex conjugate͒ and taking the velocity-space integral and summation over k. Here, Q i , W, and D are defined as the perpendicular ion heat flux  16 In the collisionless system with ϭ0, as we have shown in Ref. 7 for the no zonal flow case, the quasisteady state characterized by monotonical increase of ␦S is realized in turbulence while keeping W and Q i constant in average, that is where¯indicates time averaging on a certain period longer than a characteristic time of the turbulence. On the other hand, even if is much smaller than inverse of the characteristic time of instabilities, introduction of the collision term may lead to statistically steady turbulence, where not only low-order moments but also the distribution function itself are statistically steady. Therefore, in the case with finite collisionality, it is expected that d(␦S)/dtϷ0 and

͑6͒
In order to examine effects of the finite collisionality, we have performed several simulations for different 's. Throughout the simulation runs shown below, we set i ϭ10 and ⌰ϭ2.5. For these parameters, in the collisionless case, the angular frequency r and the linear growth rate i of the most unstable mode with k x ϭ0.1 and k y ϭ0.3 are r ϭϪ0.957 and i ϭ7.73ϫ10 Ϫ2 ͑normalized by v ti /L n ), respectively, of which changes due to finite values of are negligible in the present parameter range. ␦S given by the weakly collisional simulation of the ITG turbulence for ϭ1.25ϫ10 Ϫ4 is plotted as a function of time in Fig. 1, where the collisionless simulation result is also shown as a reference. Time evolutions of ␦S for different values of ͓which is changed from (1/512)ϫ10 Ϫ3 to 8ϫ10 Ϫ3 ] are similar to that for ϭ1.25ϫ10 Ϫ4 in Fig. 1. In the quasisteady state of the collisionless turbulence, one finds the monotonical increase of ␦S, while the potential energy W and the heat flux Q i are saturated. 7 In the finite collisionality case, however, the growth of ␦S ceases in the turbulence as well as W and Q i . It means that not only the low-order but also the highorder moments of ␦ f are statistically steady in the weakly collisional case, since ␦S is also represented by the sum of squares of the velocity space moments, ␦Sϭ͚ n ␦S n ͓see Eq. ͑8͒ for definition of ␦S n ]. 6 Time histories of each term in Eq. ͑4͒, d(␦S)/dt, dW/dt, i Q i , and D for ϭ1.25ϫ10 Ϫ4 , are plotted in Fig. 2, where data are running averaged for a time period of ϭ50. One can see that the collisional dissipation D balances with the mean transport, that is, i Q i ϷϪD , while d(␦S)/dtϷdW/dtϷ0. This also means that the weakly collisional turbulence is in the real statistically steady state. The high accuracy in calculation of the entropy balance is achieved by the sufficient velocity-space resolution of the Eulerian kinetic simulation.

B. Collision frequency dependence of transport
Collision frequency dependence of the ion heat transport coefficient, i ϵQ i / i , is summarized in Fig. 3, where the time average is taken from tϭ1000 to 3000. According to the value of , the time step ⌬t is changed from 1/80 to 1/320 so that the numerical error in Eq. ͑4͒ should be much smaller than i Q i and D . The error bars are estimated from the standard deviation of running-averaged i for a time period ϭ10. In a range of 1.25ϫ10 Ϫ4 ϽϽ8ϫ10 Ϫ3 , i has a logarithmic dependence on . The dependence of i becomes quite weak for lower collision frequencies (р1.25 ϫ10 Ϫ4 ), where i approaches a level of the collisionless one shown by a horizontal dashed line in Fig. 3, that is, i Ϸ0.36 i 2 v ti /L n . From the results shown in Figs. 1-3, it is summarized that, if is small enough, then the collision term does not influence the low-order moments of ␦ f as well as the transport coefficient i , while it is indispensable to realizing the statistically steady turbulence through damping of the high-order ones generated by the phase mixing. These facts agree with a concept that the quasisteady state is regarded as an idealization of the real steady state in the weakcollisionality limit. 6

C. Convergence check to mode truncation
We have also carried out simulation runs with different values of the maximum wave number k max in mode truncation such as k max ϭ1.6, 3.2, 4.8, 6.4, and 12.8, for ϭ1.25 ϫ10 Ϫ4 . The results are summarized in Fig. 4 in terms of the ion heat transport coefficient i . The numerical simulation is carried out up to tϭ3000, and the time-averaging period is the same as those in Fig. 3. However, the simulation with k max ϭ12.8 could only be run up to tϭ1200, because it requires a large amount of computational cost. The time aver-aging is, thus, taken from tϭ1000 to 1200, although it does not affect the results. One can see good convergence of i for k max у3.2, while i for k max ϭ1.6 is about 30% larger than the others. One of the reasons is that, in the case of k max ϭ1.6, the potential amplitude at kϭk max is not sufficiently damped by the finite Larmor radius ͑FLR͒ effect that is represented as exp(Ϫk 2 /2) in the definitions of k and ⌿ k . Convergence of the simulation results to k max is also discussed in the next section.

A. Theoretical framework
The parallel advection term in Eq. ͑1͒ generates finescale fluctuations of the perturbed distribution function ␦ f in the velocity space ͑phase mixing͒. The small-scale components of ␦ f are effectively damped by the collision term with the second derivative in v ʈ . Here, we investigate a power spectrum of the velocity distribution function ␦ f . For the spectral analysis of ␦ f , it is meaningful to pay attention to the transfer of the entropy variable ␦S from macro to micro velocity scales. Using the basic equations, Eqs. ͑1͒ and ͑2͒, we obtain d dt ͫ ␦S n ϩ␦ n, 1  n!͉ f k,n ͉ 2 , ͑8͒ J nϪ1/2 ϵ ͚ k ⌰k y n! Im͑ f k,nϪ1 f k,n * ͒, In the steady state, the left-hand side of Eq. ͑7͒ vanishes. We see that J nϪ1/2 (J nϩ1/2 ) represents the entropy transfer from the (nϪ1)th (nth) to the nth ͓(nϩ1)th͔ Hermitepolynomial portion. The third and fourth terms on the righthand side of Eq. ͑7͒ represent the entropy production due to the downward turbulent heat flux in the temperature gradient and the collisional dissipation, respectively. It is important to note that, in the range nу3, the entropy production rigorously disappears, which is the reason why the Hermitepolynomial expansion is employed here. A clear cutoff of the entropy production like this never occurs if we use the Fourier expansion in terms of exp(ilv ʈ ) (ϪϱϽlϽϱ) as basis functions.
Note that where the function e Ϫl 2 /2 l n of l has the maximum absolute value at lϭϮͱn and is expanded around it as e Ϫl 2 /2 l n Ӎe Ϫn/2 n n/2 ͓1Ϫ͑ lϯͱn ͒ 2 ͔. ͑13͒ Thus, the main contribution to the nth component of the Hermite expansion is from the Fourier components with ͉l/ͱnϯ1͉Ͻ1/ͱn, so that we may use the relation nӍl 2 for nӷ1. The inverse of Eq. ͑12͒ is given by The phase mixing process described in Eq. ͑1͒ causes the factor exp͓il(t)v ʈ ͔ in the distribution function, where l(t) satisfies dl(t)/dtϭϪ⌰k y . Then, the sign of l(t) is opposite to that of k y for large t so that we write l(t) ӍϪ(k y /͉k y ͉)ͱn(t) with the order n(t) of the Hermitepolynomial expansion as a function of t.
For nу3 in the steady state, we find from Eq. ͑7͒ that where n is treated like a continuous variable and the finite difference is approximated by the derivative. The ratio of J n to ␦S n is written as where (nϩ1/2)!/n!ϭ⌫(nϩ3/2)/⌫(nϩ1)Ӎͱn for large n and the averaging operator ͗•͘ n ϵ( ͚ k ͉ f k,n ͉ 2 •)/(͚ k ͉ f k,n ͉ 2 ) are used. In Eq. ͑16͒, we put f k,nϪ1/2 f k,nϩ1/2 * Ӎi(k y /͉k y ͉)͉ f k,n ͉ 2 by assuming the n dependence of the phase of f k,n to be described mainly by Eq. ͑14͒ with l ӍϪ(k y /͉k y ͉)ͱn.
Now, let us examine the role of the EϫB convection term. In analogy with the study by Batchelor on the spectrum of the passive scalar for wavelengths smaller than the Kolmogorov scale in the large Prandtl number case, 10 we postulate that, for large n, f k,n varies so rapidly that EϫB flow acting on f k,n is regarded as a steady one which is statistically independent of f k,n . Then, the strain of the steady flow is considered to cause the exponential growth of the wave number of the convected variable. Thus, writing k y (t)ϰe ␥t and using ͱnϭϪ(k y /͉k y ͉)lϭ⌰͐dt͉k y ͉, we have ⌰͉k y ͉ Ӎ␥ͱn, which is substituted into Eq. ͑16͒ to yield J n /␦S n ϭ2␥n. ͑17͒ Substituting Eq. ͑17͒ into ͑15͒ and integrating it with respect to n, we obtain where coincides with the entropy production rate as shown by the constraint derived from Eq. ͑7͒ ϭ2 ͵ 0 ϱ n␦S n dnӍ2 ͚ n n␦S n ϭ i Q i . ͑19͒ Thus, in the range where neither entropy production nor collisional dissipation occurs (1ӶnӶ␥/), we expect the power law of ␦S n ϰ1/n with J n ϭϭconst which is analogous to the passive scalar spectrum and its power transfer in the viscous-convective subrange.
In the above analytical treatment, ͉͗k y ͉͘ n ϰͱn increases infinitely with n. However, in the numerical simulation performed in the present study, there exists the upper limit of ͉k͉. Even if the potential amplitude is sufficiently damped at the maximum wave number in the simulation, still f k,n for large ͉k͉ and large n is continuously produced by the combination of the EϫB convection and the phase mixing process.
Therefore, saturation of ͉͗k y ͉͘ n with increasing n is anticipated due to the upper limit of ͉k͉. In this case, taking ⌰͉͗k y ͉͘ n ϭ␥ M as independent of n, we obtain from Eqs. ͑15͒ and ͑16͒ where is the same as given by Eq. ͑19͒.

B. Application to the simulation results
Let us compare the theoretical analysis given above with the simulation results. For different collision frequencies , we have plotted ␦S n and J nϪ1/2 / in Figs. 5 and 6, respectively. The maximum wave number is k max ϭ3.2. One can see that ␦S n at nϭ1, 2, and 3 has almost the same values for different , and is followed by a power-law profile but with -dependent exponents ͓␦S n ϰn Ϫ␣ where ␣()Ͼ0]. Then, ␦S n exponentially decays in the dissipation range, where the collision term is dominant, of which a typical value of n is scaled as n d ϰ Ϫ2/3 in consistent with Eq. ͑21͒. Thus, k max affects the spectra for high n's.
Profiles of ␦S n obtained by simulations depend on as well as k max and seem to be described by a mixture of Eqs. ͑18͒ and ͑21͒. Certainly, for sufficiently small (р1.25 ϫ10 Ϫ4 ), we have confirmed 0.5Ͻ␣Ͻ1. If the collisional dissipation on the left-hand side of Eq. ͑15͒ is neglected, one finds that J nϪ1/2 Ӎconst. In correspondence, for р1.25 ϫ10 Ϫ4 , there exists a flat profile of J nϪ1/2 where neither the entropy production nor the collisional dissipation is seen. For larger , the subrange with constant J nϪ1/2 ends at relatively small n, which means that the entropy production region on the low-n side and the dissipation range are not separated well. Thus, lower-order moments and transport are influenced by . Therefore, it results in the spectrum of ␦S n with ␣Ͼ1 for the relatively large values of (Ͼ1.25ϫ10 Ϫ4 ). It is also noteworthy that the logarithmic dependence of i on is observed in Fig. 4 for the same parameter range of as that discussed here.
The entropy variable ␦S k,n before taking summation over k in Eq. ͑7͒ is transferred in the (k,n) space by the phase mixing as well as the EϫB advection. The latter effect can be found in profiles of ␦S k y ,n ϵ ͚ k x ␦S k,n , of which cross-sectional plots are shown in Fig. 7 for the case with ϭ1.25ϫ10 Ϫ4 and k max ϭ12.8. A steeply peaked profile of ␦S k y ,n at around k y ϭ0 in nϽ100 ͑as shown by the upper three lines in Fig. 7͒ disappears with increasing n. This is caused by the growth of ͉͗k y ͉͘ n due to the strain of the E ϫB flow. Then, the whole profile gradually broadens, until the saturation of ͉͗k y ͉͘ n occurs due to the finite k max .
The growth and saturation of the spectrum-averaged wave number ͉͗k y ͉͘ n defined in Eq. ͑16͒ is more clearly recognized in Fig. 8. For k max ϭ12.8, ͉͗k y ͉͘ n grows nearly in proportion to ͱn from nϭ3 to ϳ100, which agrees with the estimate of ⌰͉k y ͉Ӎ␥ͱn used for derivation of Eq. ͑17͒. As n increases, then, the growth of the wave number slows down due to the upper limit k max . A similar evolution of ͉͗k y ͉͘ n is also found for smaller k max , although the slow-  Even if k max is small, thus, the values of ␦S n at nϭ1, 2, and 3 are unchanged, except for the case of k max ϭ1.6, where ͉͗k y ͉͘ n increases slower than ͱn. The above result is consistent with the convergence check of the transport coefficient i to k max given in Sec. III C. Therefore, we can conclude that the obtained i as well as the entropy transfer process on the macro velocity scale but in a small ͉k͉ region is insensitive to the maximum wave number of k max у3.2.
It is meaningful to estimate ␥ and ␥ M from the simulation result. Let us suppose ͉͗k y ͉͘ n ϳ0.2ͱn according to Fig.   8, and substitute it into ⌰͉͗k y ͉͘ n Ӎ␥ͱn. Thus, we find ␥ ϳ0.5 for ⌰ϭ2.5. In the dissipation range of nտ10 3 , however, ͉͗k y ͉͘ n ϳ6 for k max ϭ12.8 as shown in Fig. 8, and thus, ␥ M ϭ⌰͉͗k y ͉͘ n ϳ15. The spectrum ␦S n obtained by the simulation with k max ϭ12.8 is then compared with those given by the theoretical analysis with ϭ36 in Fig. 9, where the solid, dashed, and dotted lines indicate the simulation result, Eq. ͑18͒ with ␥ϭ0.5, and Eq. ͑21͒ with ␥ M ϭ15, respectively. One can see that the simulation result is nearly proportional to that of Eq. ͑18͒ in a range of 3рnՇ10 3 . In this range n Շ10 3 , the constant factor /2␥ in Eq. ͑18͒ gives smaller values of ␦S n than that obtained by the simulation. The reason for this is understood as follows. We should recall that the constant factor is derived from the constraint 2͐n␦S n dnϭ, where is evaluated from the simulation result. Then, in the range nՇ10 3 , the spectrum ␦S n in Eq.
͑18͒ needs to be smaller than that in the simulation in order to satisfy the constraint because, for nտ10 3 , the latter spectrum is smaller than the former due to the effect of finite k max . However, in the dissipation range of nտ10 3 , the spectrum found in the simulation is well fitted by Eq. ͑21͒. The above results show that the transfer process of the entropy variable observed by the Eulerian kinetic simulations of the slab ITG turbulence can be well described by combining the analytical expressions in Eqs. ͑18͒ and ͑21͒ in Sec. IV A. It is expected that, if we can employ a sufficiently high value of k max , the simulation will reproduce the spectrum ␦S n in Eq.
͑18͒ for the whole range nу3.

V. CONCLUDING REMARKS
We have carried out Eulerian kinetic simulations of the slab ITG turbulence with weak collisionality, where we have employed the gyrokinetic equation ͑integrated for v Ќ ) with the Lenard-Bernstein model collision operator and the quasineutrality condition. Introduction of finite collisionality enables us to find the real statistically steady state of turbulence, where not only the turbulence energy and the transport flux but also the entropy variable ␦Sϵ͐␦ f 2 /2F M dv are constant in average. Then, the ion heat transport flux Q i multiplied by i is balanced with the collisional dissipation. It is in contrast to the quasisteady state of the collisionless turbulence, where d(␦S)/dt balances with i Q i but with constant energy. 7 A parameter survey for the collision frequency shows the logarithmic dependence of i on relatively large values of . For sufficiently low collision frequency, however, the transport coefficient i approaches a value in the collisionless case, which means that the low-order velocityspace moments of the distribution function in the quasisteady state of the collisionless turbulence agree with those in the real steady state of weakly collisional one. It is also confirmed that i has a good convergence to the maximum wave number k max у3.2 employed in the simulations.
We have also done the spectral analysis of the distribution function in the velocity space by means of the Hermitepolynomial expansion with the Maxwellian weight function.
The entropy variable ␦S n defined by a power spectrum of the distribution function has almost the same values at nϭ1, 2, and 3 even for different values of , where n denotes the order of the Hermite-polynomial expansion. This is consistent with the dependence of i for the low-collisionality case described above as well as the conjecture by Krommes and Hu such as ''flux determines dissipation.'' 4 The entropy variable produced on the low-n side is transferred in the (k,n) space by the phase mixing and the EϫB nonlinearity. Then, it is dissipated by the collision term on a high-n side. Our theoretical analysis of the spectrum describes well how ␦S n (nӷ3) depends on n, , and k max . In analogy with the passive scalar convection for wavelength smaller than the Kolmogorov scale in the large Prandtl number case, the scaling law, ␦S n ϰ1/n, is found in the subrange of the n space where neither entropy production nor collisional dissipation occurs.
The two-dimensional slab ITG turbulence is considered in the present study, while the entropy balance equation similar to Eq. ͑4͒ can also be derived from the toroidal gyrokinetic equation. 17 The relationship between the steady and quasisteady states shown in Sec. III is expected to be valid in the toroidal configuration as well, although further investigations are required. In the toroidal gyrokinetic case, the phasemixing process becomes more complicated because the toroidal magnetic drift should be taken into account. Extension of the present study to the toroidal configuration is currently in progress and will be reported elsewhere.