Nondissipative kinetic simulation and analytical solution of three-mode equations of the ion temperature gradient instability

A nondissipative drift kinetic simulation scheme, which rigorously satisfies the time-reversibility, is applied to the three-mode coupling problem of the ion temperature gradient ~ITG! instability. It is found from the simulation that the three-mode ITG system repeats growth and decay with a period which shows a logarithmic divergence for infinitesimal initial perturbations. Accordingly, time average of the mode amplitude vanishes, as the initial amplitude approaches zero. An exact solution is analytically given for a class of initial conditions. An excellent agreement is confirmed between the analytical solution and numerical results. The results obtained here provide a useful reference for basic benchmarking of theories and simulations of the ITG modes. © 2000 American Institute of Physics.@S1070-664X~00!02003-6#


I. INTRODUCTION
Understanding the anomalous heat transport mechanism in high-temperature plasmas has been a central subject in magnetic confinement fusion studies.The ion temperature gradient ͑ITG͒ instability 1,2 is widely recognized as one of the candidates for the cause of the anomalous ion thermal transport in the core of tokamaks.Many first-principle simulations, as well as theoretical predictions, have been carried out on ITG turbulence transport.Among them, the gyrokinetic [3][4][5][6] and gyrofluid [7][8][9][10] simulations have largely contributed to development of transport modeling in the last decade.However, they have a discrepancy in the ion thermal diffusivity by a factor of two or more. 5,11Since the reason for the difference between the two methods has not yet been completely understood, simple nonlinear problems, for which reliable solutions with sufficient accuracy are established, are preferable for benchmark studies of the numerical schemes.
The three-mode coupling of the ITG mode, 3 as well as the drift waves, [12][13][14] which is represented by a reduced set of the drift kinetic equation in a two-dimensional slab geometry, has been used to explain qualitatively a nonlinear saturation mechanism in the gyrokinetic simulation. 3The threemode system has also been examined in comparison of the gyrokinetic and gyrofluid simulations 11 and in a test of a nonlinear kinetic fluid closure method. 15In the pioneering work by Lee et al., 12 they derived a nonlinear dispersion relation of the three-mode drift waves.However, so far, no exact solution of the three-mode coupling equations has been obtained analytically, except for a steady state solution with the Maxwellian velocity distribution. 143][14][15] The Vlasov simulation results, however, seem to suffer from numerical dissipation, because a dissipative integrator such as the predictor-corrector is employed for the collisionless drift kinetic equation.Therefore, it would be meaningful to develop a nondissipative Vlasov simulation method, and to find nu-merical and analytical solutions of the three-mode equations.Results of the present study give a useful reference for benchmarking of various simulations and theories which are employed to investigate more complicated systems.
The remainder of the present paper is organized as follows.In Sec.II, we will analytically derive an exact solution for a certain class of initial conditions from the three-mode ITG equations.Section III gives the numerical simulation results of the three-mode coupling of the ITG modes, where a newly developed drift kinetic simulation method is also briefly explained.Comparisons of the numerical and theoretical results are also presented in Sec.III.Finally, we summarize the results in Sec.IV.

A. Model configuration
We start from the electrostatic drift kinetic equation in a slab geometry, where the long wavelength limit (k Ќ i Ӷ1) is assumed, and v EϫB ϭEϫB/B 2 .v is the parallel velocity.Here, we consider the same model as that used in comparison of the gyrokinetic and the gyrofluid simulations 11 except for truncation of the distribution function and a /2 phase difference in x ͑direction of density and temperature gradients͒.The same model was also used by Mattor and Parker for examining a nonlinear kinetic fluid closure. 15A rectangular domain of L x ϫL y is considered in the x-y plane with a uniform external magnetic field perpendicular to the x-axis.Neglecting the parallel nonlinearity and expanding density and temperature scale length, L n and L T , of an assumed Maxwellian background, F M (v), the drift kinetic equation for ions, Eq. ͑1͒, leads to where f ˜denotes a perturbed distribution function normalized by f ˜ϭ f ˜ЈL n v ti / i n 0 .Prime means a dimensional quantity.v ti , i , and n 0 are the ion thermal velocity, the ion gyroradius, and the background plasma density.⌰ is defined as ⌰ϭL n / i , where an inclination of the magnetic field Ӷ1 is assumed.Other quantities are normalized as x ϭxЈ/ i , yϭyЈ/ i , vϭvЈ/v ti , tϭtЈv ti /L n , i ϭL n /L T , and ϭeЈL n /T i i with the elementary charge e and the background ion temperature T i (ϭm i v ti 2 ; m i means the ion mass͒.We have taken T i ϭT e throughout this paper.We also assume the adiabatic electron response and the quasineutrality.Thus, We employ the periodic boundary conditions in both x and y directions.Then, f ˜and can be written as where k x ϭ2m/L x and k y ϭ2n/L y for mϭ0,Ϯ1,Ϯ2, . . .and nϭ0,Ϯ1,Ϯ2, . . . .We also take L x ϭL y .For studying the three-mode coupling, we only keep (m,n)ϭ(Ϯ1,Ϯ1) and (Ϯ2,0) modes with the following symmetry conditions of f ˜1,1 ϭ f ˜Ϫ1,1 ϭ f ˜1,Ϫ1 * ϭ f ˜Ϫ1,Ϫ1 * and f ˜2,0 ϭ f ˜Ϫ2,0 * .Also, Re(f ˜2,0 )ϭ0.

B. Three-mode equations and linear solution
We derive exact nonlinear solutions of the three-mode problem of the ITG modes.Hereafter, in Sec.II, f ˜1,1 , Im(f ˜2,0 ), and 1,1 are, respectively, denoted by f , h, and for simplicity.From Eqs. ͑2͒ and ͑3͒ the three-mode ITG system is described by 3

͑8͒
where f and are complex-valued while h is real-valued.G(v) is defined as

͑9͒
We easily find The linearized version of the three-mode equations has a linear solution of the form Here, the linear eigenfunctions are given by h L ͑ v ͒ϵ0, and L ϵ1͑normalization͒, ͑12͒ and the complex eigenfrequency ϭ r ϩi␥ is determined by the dispersion relation where ␥Ͼ0 is assumed.

C. Nonlinear solution
Now, let us consider a certain class of exact solutions of the nonlinear three-mode ITG equations, which are written in terms of the real and imaginary parts of the eigenfunction f L (v) and the real eigenfrequency r as where a(t), b(t), and c(t) are real-valued functions of the time t.The linear solution given by Eqs.͑11͒-͑13͒ corresponds to the case in which a(t)ϭb(t)ϰexp(␥t) and c(t) ϭ0.Substituting these into Eqs.͑6͒-͑8͒ and using Eq.͑11͒, we obtain the ordinary differential equations for These equations have two types of stationary solutions, which are written as respectively, where a s and c s are arbitrary constants.The former solution in Eq. ͑16͒ is the c-axis, and the latter one is parallel to the a-axis.From Eq. ͑15͒, we easily find

͑18͒
Also, combining these two equations, we have

͑19͒
which corresponds to Eq. ͑10͒.Thus, as shown in Fig. 1, the orbit of the solution in the (a,b,c)-phase space is given by the intersection between the two surfaces, which are written as where C 1 and C 2 are constants.The orbit is also on the spheroid surface, We also note that the stationary solutions (a s ,0,␥/2k 2 ) with varying a s form the central axis of the elliptic column given by Eq. ͑21͒, which is also shown in Fig. 2. Using Eqs.͑20͒-͑22͒, the solution of Eq. ͑15͒ for the initial condition (a,b,c) tϭ0 ϭ(a 0 ,b 0 ,c 0 ) can be obtained by and

͑25͒
where the Jacobi elliptic functions are defined by dnuϭ(1 Here, the parameters ␣, ␤, and 2 are given by ␣ϭ Ά respectively, where The solutions given above are periodic functions of the time t.For rϾ0, a(t) takes any value in the range ␤(r/␣) 1/2 р͉a͉р␤, and its sign is a constant determined by the initial value a 0 .On the other hand, for rϽ0, the value FIG. 2. A separatrix surface of rϭ0 for the same parameters as Fig. 1 with stationary solutions of (0,0,c s ) and (a s ,0,␥/2k 2 ), where two solid curves represent typical solutions for rϾ0 and rϽ0.
range of a(t) is given by Ϫ␤рaр␤ and the sign of a changes with time.Thus, the shape of the closed orbit in the (a,b,c)-phase space is like a butterfly for rϽ0 ͑see Figs. 2  and 9͒.The period T of the solution ͓a(t),b(t),c(t)͔ is written as where K() is the complete elliptic integral of the first kind.We find that, as (a 0 ,b 0 ,c 0 )→(0,0,0), r→0 and Tϳ ͭ ␥ Ϫ1 log͑1/͉r͉͒ for r→ϩ0 2␥ Ϫ1 log͑1/͉r͉͒ for r→Ϫ0 .͑33͒ Thus, the period T shows the logarithmic divergence for r →0.The solution ͓a(t),b(t),c(t)͔ stays most of the period in the neighborhood of the stationary point ͑0,0,0͒.For the case of rϭ0, the solution is no longer periodic ͑or Tϭϱ) and its orbit emerges from ͑approaches to͒ the stationary point ͑0,0,0͒ for t→Ϫϱ (ϩϱ).The set of (a 0 ,b 0 ,c 0 ) satisfying rϭa 0 2 Ϫ (2k 2 /␥) a 0 2 c 0 ϩ (2k 4 /␥ 2 ) a 0 4 Ϫb 0 2 ϭ0 forms a separatrix surface, which separates the phase space into the two regions (rϾ0 and rϽ0; see Fig. 2͒ which are filled with the two different types of orbits given by Eqs.͑23͒-͑25͒.

A. Simulation scheme
A basic scheme of our drift kinetic simulation is briefly presented here.It is easily found that Eq. ͑1͒ is reversible in time.Avoiding the numerical noise inherent in the particle simulation, we employ an Eulerian scheme which keeps the time-reversibility.A discrete spectral representation in the phase space or discretization on numerical grids makes Eq. ͑1͒ a set of ordinary differential equations.In a vector form, it can be written as In order to keep the time-reversibility of Eq. ͑1͒, a numerical time-integration scheme of Eq. ͑34͒ should also be reversible in time.It is well known that the symplectic scheme often used for integration of a Hamiltonian system, preserving a symplectic 2-form exactly, is nondissipative, namely, time-reversible. 16,17One of the simplest examples is the leapfrog integrator, which is a standard scheme in particle simulations, 18 where particle motions are given by the Newton-Lorentz equation ͑not by drift motion of the guiding center͒.To trace the EϫB drift particle motion with the time-reversibility, however, one needs to use an implicit symplectic method.This is because the Hamiltonian, H(x,y,z,v)ϭmv 2 /2ϩe(x,y,z), for the drift motion is nonseparable for perpendicular coordinates x and y, which are a conjugate pair of the coordinates ͑where dx/dtϭ Ϫ‫ץ/ץ‬y and dy/dtϭ‫ץ/ץ‬x), while it is separable for the another conjugate pair, that is, the parallel coordinate z and the parallel velocity v.One of the simplest implicit schemes is the implicit midpoint rule, where n means a time step.This scheme is apparently reversible in time.When Eq. ͑35͒ is applied to the Hamilton's equation in canonical coordinates, Uϭ(q,p), it leads to a canonical transform from (q n ,p n ) to (q nϩ1 ,p nϩ1 ). 19e have employed Eq. ͑35͒ to integrate the drift kinetic equation, where f ¯ϭ( f nϩ1 ϩ f n )/2 and ͕,͖ means the Poisson brackets.
H ¯depends on f ¯through the electrostatic potential .Although Eq. ͑36͒ is not a symplectic transform of f generated by a particle Hamiltonian, it preserves the time-reversibility, namely, is nondissipative.It is also noteworthy that Eq. ͑36͒, which can be solved by iteration, is regarded as a discretized form of ‫ץ‬ f ¯/‫ץ‬tϭϪ͕ f ¯,H ¯͖ with second-order accuracy in time ͑namely, a time-centered finite difference͒.Construction of a fourth-order scheme is straightforward by successive operations of Eq. ͑35͒. 20he nondissipative simulation scheme given above is applied to the three-mode coupling system of the ITG modes in the two-dimensional shearless slab geometry.Starting with an initial condition of for L x ϭL y ϭ2/k, we numerically follow a time evolution of f ˜in Eq. ͑2͒.
The physical parameters used here are kϭk x ϭk y ϭ0.1 and i ϭ10.We have carried out several runs for different ⌰ scanning from 0.25 to 3. Amplitude of the initial perturbation, , is also changed from 10 Ϫ5 to 1.The box size of the simulation is L x ϭL y ϭ20 with 32ϫ32 grid points.Spatial derivatives in Eq. ͑2͒ are calculated in the Fourier space.The velocity space of Ϫ5рvр5 is discretized by 129 grid points.A time step is taken to be ⌬tϭ0.25.
We have made convergence checks for the time step, the resolution in the velocity space, the maximum velocity, and the accuracy of integration scheme, all of which give the same results as shown below.

B. Linear results
In Fig. 3, the linear frequency r and growth rate ␥ resulting from simulations for different ⌰ are shown by circular and triangular marks in comparison with the linear analysis.Simulation results shown here are obtained by the second-order implicit scheme in Eq. ͑36͒.Solving the linear dispersion relation given in the last section ͓see Eq. ͑13͔͒ numerically, we have calculated the theoretical values, which agree with the simulation results.By the Nyquist criterion, we have also confirmed that only one eigenmode is unstable for a given value of ⌰.
In Fig. 4 are plotted profiles of f ˜1,1 for ⌰ϭ1 and ϭ10 Ϫ5 in the velocity space during the linear growth phase at tϭ100.Here, a phase of f ˜1,1 is shifted in y so that 1,1 is real.Amplitude of f ˜1,1 ͑and also f ˜2,0 in Fig. 6͒ is normalized by 1,1 .Circular and triangular marks in Fig. 4 show the simulation data on each grid point.Solid and dashed lines represent real and imaginary parts of a linear eigenfunction of the ͑1,1͒ mode, f Lr (v) and f Li (v), for the eigenfrequency r ϩi␥, which are defined in Eqs.͑11͒ and ͑13͒ in Sec.II.In the linear growth phase, f ˜1,1 is well fitted by the eigenfunction, while f ˜2,0 is negligible.

C. Nonlinear results
A time evolution of mode amplitudes for ⌰ϭ1 and ϭ10 Ϫ5 is shown in Fig. 5, where the 0th-moment of f ˜1,1 and 2nd-moment of f ˜2,0 are represented by solid and dashed lines, respectively.After the initial linear growth phase, the amplitude of the ͑1,1͒ mode peaked at tϭ196 due to appearance of the nonlinear ͑2,0͒ mode.Then, the ͑1,1͒ mode exponentially decays with the same rate as the linear growth phase, and reaches its minimum at tϭ393.The minimum amplitude is nearly equal to the initial perturbation level.After that, the growth and decay are repeated with a regular oscillation period of Tϭ394.
Figure 6 shows profiles of f ˜1,1 and f ˜2,0 in the velocity space in the nonlinear phase.When the ͑1,1͒ mode amplitude is peaked at tϭ196, Im(f ˜1,1 ) disappears ͑upper panel͒.On  FIG. 6.The same as Fig. 4, but for tϭ196 ͑upper͒ and tϭ290 ͑lower͒.
the other hand, one can see a finite amplitude of Im(f ˜2,0 ), of which profile is scaled as & f Li (v).During the decay phase, f ˜1,1 is fitted by the complex conjugate of f L (v).Neglecting small fluctuations of order , therefore, the linear and nonlinear evolutions of f ˜1,1 and f ˜2,0 are described in terms of f Lr and f Li as is considered in the last section.
The peak amplitude of 1,1 provides a good benchmark for the simulation scheme, since the same test has been done by the gyrokinetic and gyrofluid codes 11 as well as theoretical predictions. 3,11A systematic scan for ⌰ from 0.25 to 3 has been made with ϭ10 Ϫ5 .The peak levels are summarized in Fig. 7 with a theoretical prediction of ͉ 1,1 ͉ peak ϭ␥/&k 2 given as follows: Substitute C 1 ϭC 2 ϭ0 to Eqs. ͑20͒ and ͑21͒ for infinitesimal initial perturbation.Since da/dtϭ0 at the peak of the mode amplitude, bϭ0 ͓see Eq. ͑15͔͒.Then, one finds aϭ␥/&k 2 , bϭ0, and cϭ&a ϭ␥/k 2 , which also gives the scaling of Im(f ˜2,0 )ϭ& f Li shown in Fig. 6.The simulation and the theory give the same peak level and scaling.
The scaling of ͉ 1,1 ͉ peak ϭ␥/&k 2 agrees with previous theoretical works 3,11 except for a factor of &.The peak amplitude for ⌰ϭ1 is almost the same as the drift kinetic simulation shown in Ref. 15.In the early work on the gyrokinetic and gyrofluid comparison, 11 the peak amplitude given by the gyrokinetic code is about 15% lower than the present results but with the similar dependence on ⌰.The gyrofluid code shows an agreement at small ⌰, but overestimates at large ⌰.Since f ˜or density and temperature perturbations were not truncated in Ref. 11, more detailed benchmarking would be desired.
As found in Fig. 5, after the peaking, the mode amplitude decreases to the initial perturbation level.Thus, it is considered that the oscillation period T depends on .3][14] In order to examine dependence of the oscillation period on the initial perturbation amplitude, we have performed six runs for ⌰ϭ1, changing as 10 Ϫ5 , 10 Ϫ4 , . . ., 1.The observed periods T in the simulations are plotted in Fig. 8 versus the minimum amplitude ͉ 1,1 ͉ min , which shows a logarithmic dependence of T on ͉ 1,1 ͉ min .The analytical solution for rϾ0 given by Eq. ͑32͒ predicts the same period.The ͉ 1,1 ͉ min dependence of T is easily explained by noting that a smaller initial perturbation stays for a longer time at the linear growing and decaying phases.Thus, a time average of ͉ 1,1 ͉ approaches zero as →0, which is consistent with Eq. ͑33͒.
Finally, Fig. 9 shows the two kinds of orbits in Eqs.͑23͒-͑25͒ for rϾ0 with (a 0 ,b 0 ,c 0 )ϭ(0.25,0,0)͑upper͒ and rϽ0 with (a 0 ,b 0 ,c 0 )ϭ(0,0.25,0)͑lower͒ by solid lines.We also plot circular marks representing the Vlasov simulation results at every 2L n /v ti which are, respectively, begun with initial conditions of As shown in Fig. 9, a perfect agreement is found between the analytical solutions and the Vlasov simulation results.It is clearly confirmed that our numerical scheme strictly preserves the time-reversibility, the periodicity, and the initial value dependence of the exact solution.Here, we have assumed the initial conditions in Eqs.͑38͒ and ͑39͒ which satisfy Eq. ͑14͒, while a different type of the initial condition in Eq. ͑37͒ has been employed for the simulations shown in Figs.3-8.Nevertheless, all of the simulation results are well described by the analytical solution.It is because, if the initial perturbation is sufficiently small, only one linearly unstable eigenmode in the form of Eq. ͑14͒ becomes dominant at the early time stage.Thus, the amplitude oscillation found in Fig. 5, where aϾ0, is explained by the periodic motion along the orbit for rϾ0 in the (a,b,c)-phase space.

IV. SUMMARY
We have developed a nondissipative Eulerian kinetic simulation method which rigorously preserves the time-FIG.7. Peak amplitudes of 1,1 for different ⌰ ͑marks͒ with a theoretical prediction of ␥/&k 2 ͑solid͒ where ϭ10 Ϫ5 with kϭ0.1 and i ϭ10.
reversibility inherent in the collisionless drift kinetic equation.The simulation method is applied to the three-mode ITG problem.An analytical solution of the three-mode ITG equations is also derived here, and successfully describes the nonlinear evolution of the three-mode system.We have found an excellent agreement between the simulation results and the analytical solution.
The main results obtained here are summarized as follows: Amplitudes of (Ϯ1,Ϯ1) and (Ϯ2,0) modes repeat the nonlinear oscillation, and evolve in time along a closed orbit which encircles the stationary solution (a s ,0,␥/2k 2 ) in the phase space.This fact reflects the time-reversibility of Eq. ͑2͒.Thus, the three-mode system never reaches a steady state, unless the initial condition coincides with the stationary solutions in Eq. ͑16͒.Since the solution stays most of the period in the neighborhood of the stationary point ͑0,0,0͒ with the exponential time dependence, the oscillation period has a logarithmic dependence on the initial mode amplitude.Thus, the time average of the mode amplitude vanishes, as the infinitesimal initial perturbation approaches zero.
The success of the present simulation, the validity of which is confirmed by the theoretical analysis, is achieved by a proper choice of the numerical scheme with no dissipation.If the scheme is dissipative, it affects estimates of the oscillation period and the averaged amplitude.Application of the present numerical method to a multimode coupling system and the gyrokinetic equations remains for future works.

FIG. 1 .
FIG. 1.An orbit of the exact solution for ⌰ϭ1, kϭ0.1, and i ϭ10 in the (a,b,c)-phase space shown by the intersection of parabolic and elliptic column surfaces which is also on a spheroid.