Slablike ion temperature gradient driven mode in reversed shear tokamaks

The ion temperature gradient driven (ITG) mode in reversed shear tokamaks is analysed using a gyrokinetic toroidal particle code. It is found that the ITG mode in the reversed shear configuration shows a coupled mode structure between the slab and toroidal ITG modes. Especially in the qmin-region, a slablike feature due to the reversed shear slab ITG mode (Idomura Y, Tokuda S and Wakatani M 1999 Phys. Plasmas 6 4658) becomes strong. This coupled eigenmode structure is changed from a slab mode to a toroidal mode depending on ηi = Ln/Lti and Lti. Results show that in reversed shear tokamaks the ITG mode is determined from a competition between the slab and toroidal ITG modes.


Introduction
It is an important issue to analyse the ion temperature gradient driven (ITG) mode in reversed shear tokamaks for understanding transport properties of the internal transport barrier (ITB). The ITB region is characterized by a steep pressure profile and weak magnetic shear around the q min -surface. In this weak or zero magnetic shear region, the ITG mode has a global radial eigenmode structure, and the distance between the neighbouring mode rational surfaces becomes large. Thus, the scale length ordering between the radial eigenmode structure and equilibrium profiles, and the assumption of a strong toroidal mode coupling, which are used in a kinetic ballooning theory, break down. It is, therefore, difficult to apply kinetic ballooning codes to this region, and a gyrokinetic global analysis is required. To this end, we have developed a gyrokinetic toroidal particle code.
The toroidal ITG mode has been considered as a main candidate to account for the anomalous ion thermal transport in high temperature tokamak plasmas. However, in [1], it was reported that the growth rate of the toroidal ITG mode decreases as the sign of the magnetic shear is changed from positive to negative, based on a kinetic ballooning theory. On the other hand, we have shown 101.2 analytic and numerical solutions of the negative sheared slab ITG mode or the reversed shear slab ITG (RS-ITG) mode, based on a nonlocal gyrokinetic analysis in a sheared slab geometry [2]. The RS-ITG mode has a global radial eigenmode structure in a weak magnetic shear region around the q min -surface. Since a stabilizing effect of magnetic shear disappears, the RS-ITG mode becomes strongly unstable in the q min -region. This kind of global slab modes cannot be treated using a local ballooning theory. In [3], linear global analyses using a gyrokinetic global spectral code have shown that a slablike feature appears in a negative or zero magnetic shear region. The above results suggest that both the toroidal and slab ITG modes coexist in reversed shear tokamaks.
In the present work, gyrokinetic global analyses have been performed particularly for the ITG mode which is excited in the q min -region. From numerical results, it is found that the ITG mode in the reversed shear configuration shows a coupled mode structure between the slab and toroidal ITG modes. Moreover, especially in the q min -region, a slablike feature becomes strong. From a detailed study of the eigenfunction, it has been shown that this slablike feature comes from the RS-ITG mode. In local theory, the stability of the slab ITG mode is determined by η i = L n /L ti , while the toroidal ITG mode is characterized by L ti , provided that a relatively flat density profile is assumed, where L n and L ti are the characteristic scale lengths of the density and temperature gradients. In the analysis, it is found that the eigenmode structure is changed from a slab mode to a toroidal mode depending on η i and L ti . These results show that the ITG mode in the reversed shear configuration is determined from a competition between the slab and toroidal ITG modes.
The reminder of the paper is organized as follows. In section 2, model equations and numerical methods used in the gyrokinetic toroidal particle code are described. In section 3, properties of the RS-ITG mode are reviewed based on a nonlocal analysis in a sheared slab geometry. In section 4, numerical results obtained from gyrokinetic global analyses are shown and properties of the ITG mode in reversed shear tokamaks are discussed. Finally, in section 5, obtained results are summarized.

Gyrokinetic toroidal PIC code
In the present work, we have developed a gyrokinetic toroidal particle code for a threedimensional nonlinear simulation (GT3D). In this code, a nonlinear gyrokinetic Vlasov-Maxwell system [4] is solved for a circular concentric tokamak with electrostatic perturbations. In the present work, this code is applied only for linear global calculations. In this section, we will show a physical model and numerical methods used in the code.
To describe a charged particle motion in a strong magnetic field, we introduce gyro-centre coordinates, Z = (t; R, v , M, α), where R is the position of the guiding centre; ; Ω = qB/mc; m and q are the mass and charge of the particle, respectively; c is the velocity of light; b = B/B is the unit vector in the direction of the magnetic field; the gyro-phase angle is given by α ≡ tan −1 (v · e 1 /v · e 2 ) and the unit vectors, e 1 and e 2 , are chosen as e 1 × e 2 = b. In analysing low frequency fluctuations in tokamak plasmas, we can apply the usual gyrokinetic ordering: where ω is the characteristic frequency of the fluctuation, k and k ⊥ are wavenumbers in the parallel and perpendicular directions to b, respectively, T is the temperature, and ρ is the Larmor radius. Under the gyrokinetic ordering, fast non-secular perturbations related to the gyration of a charged particle is removed from a Vlasov-Maxwell system. We then have the gyrokinetic 101.3 Vlasov-Maxwell system in the gyro-averaged coordinates,Z = (t;R,v ,M,ᾱ): where B * = B + (Bv /Ω)∇R × b, B * = b · B * , φ ᾱ is the gyro-averaged electrostatic potential,F s =F 0s + δF s is the velocity distribution function,F 0s and δF s are the equilibrium and perturbed parts of the distribution function, respectively, n s and δn s are the equilibrium and perturbed parts of the particle density in a real space x,D s = |q s |m s B * /c is the Jacobian of the gyro-averaged coordinates, I n is the nth order modified Bessel function, ρ ts is the Larmor radius evaluated with the thermal velocity and s denotes the particle species. Since the phase velocity of the ITG mode, which is basically an ion sound wave, is much smaller than the electron thermal velocity, ω/k ∼ c s v te , the electron response is assumed to be adiabatic, where c s is the ion sound velocity. In the present study, we ignore the non-adiabatic response of trapped electrons, because trapped particle modes tend to be stabilized in the reversed or negative shear configuration [5]. The ion polarization density (the second term in equation (4)) is treated using a long wavelength approximation, k 2 ⊥ ρ 2 ti 1. This approximation is valid for k θ ρ ti = (nq/r)ρ ti ≤ 1, where k θ is a wavenumber in the poloidal direction, n is a toroidal mode number, q is the safety factor and r is the minor radius. For parameters used in section 4, this approximation is valid for a turbulent spectrum with n ≤ 100. From equations (4) and (5), we obtain the gyrokinetic Poisson equation for the ITG turbulence, where φ θ is the flux surface averaged potential and λ Ds is the Debye length. It is noted that the gyrokinetic Vlasov-Maxwell system, equations (1) and (6), has an energy conservation law, The magnetic configuration used in the present version of GT3D is a circular concentric tokamak configuration,

101.4
where B 0 is the magnetic field on the magnetic axis, = r/R 0 , R 0 is the major radius, ϕ is the toroidal angle, θ is the poloidal angle and e ϕ and e θ are the unit vectors in the toroidal and poloidal directions, respectively. In solving the gyrokinetic equation (1), we adopt a conventional nonlinear characteristic δf method [6]. In this method, the perturbed part δF i is solved along the nonlinear characteristics of the gyrokinetic equations (2) and (3). Although a new δf method based on the canonical Maxwellian has been developed [7], we have used a conventional δf method based on the local Maxwellian in the present linear calculations. The gyro-average for the electrostatic potential is approximated by the four-point average in a poloidal plane [8], which is sufficient for perturbations with k 2 ⊥ ρ 2 ti 1. The time integration of the nonlinear characteristics is performed by using a predictor-corrector method. To avoid difficulties arising in the treatment of particles which pass through the magnetic axis, or in the gyro-average across the magnetic axis, particle orbits are calculated in cylindrical coordinates (R, ϕ, Z), where R = R 0 + r cos(θ), and Z = r sin(θ). After the time integration, the particle position is transformed into flux coordinates (r, θ, ζ), where ζ = −ϕ. Use of the flux coordinates provides an accurate expression for k ∼ 0 perturbations as well as the edge boundary condition.
The gyrokinetic Poisson equation, equation (6), is solved using a finite element PIC method [9]. Since the system is axisymmetric, we solve each of the n components independently by expanding perturbations into a toroidal Fourier space. We then extract a ballooning phase factor or a k = 0 component analytically as where the safety factor in the straight-field-line coordinates is defined asq = −q/ √ 1 − 2 , and r s is the position of the reference magnetic surface where the mode is excited. In equation (9), the ballooning phase factor, e inqχ , is chosen to be constant in the radial direction in order to avoid a spatially varying jump condition at the boundary, θ = 0, 2π. By applying this extraction method, high m oscillations are removed from a numerical representation, and a low cost calculation of high n modes in a large device parameter is enabled, where m is a poloidal mode number. By introducing two-dimensional finite elements Λ(r, θ) with a quadratic spline function S, the slowly varying part, φ n (r, θ), is written as where k = i θ + N θ (i r − 1), i r = 1 ∼ N r , i θ = 1 ∼ N θ and N r and N θ are the system sizes in the rand θ-directions, respectively. As a boundary condition, φ n (a, θ) = 0 is imposed at the plasma edge, r = a. In order to treat the m = 0 component, which is generated from nonlinear processes, a natural boundary condition is imposed at the magnetic axis. In the poloidal direction, a jump condition, φ n (r, 0) = φ n (r, 2π)e i2πnq(rs) , is imposed at θ = 0, 2π. In the radial direction, a non-uniform mesh spacing is used to reduce the particle noise around the magnetic axis. This treatment is also useful to improve the numerical resolution in the region where the mode is excited. A matrix form of the gyrokinetic Poisson equation is then written as .
where j denotes the index of marker particles, andW j is the weight of the jth marker particle. The four-point averaging method is used also for the gyro-average in equation (15) so as not to generate a self-force. It is noted that in this algorithm, point particles are used. In a finite element PIC method, sub-grid-scale noises are suppressed by spline functions for the two-dimensional basis Λ k , while in the standard PIC simulation, such a noise is excluded by a form factor of finite size particles. In the code, the time-independent linear operator (14) is decomposed by using a Choleski decomposition at the initial time step. In each time step, φ k is obtained efficiently by calculating a back-substitution. Benchmark calculations have been performed for the following tokamak parameter: hydrogen plasma, R 0 = 1.19 m, a = 0.21 m, B 0 = 1 T and q(r) = 1.25 + 3(r/a) 2 . The density and temperature profiles are given as where T i0 = T e0 = 1 keV, L ti /R 0 = 0.13, L te /R 0 = ∞, L n /R 0 = ∞, ∆r n = ∆r t = 0.31a and r 0 = 0.5a. These benchmark parameters are the same as those used in [9]. The frequency and growth rate spectra of the ITG mode obtained from benchmark calculations are shown in figure 1. The results agree well with figure 8 in [9], and the validity of GT3D is successfully demonstrated.

Slab ITG mode in reversed shear tokamaks
In this section, we review the properties of the reversed (or negative) shear slab ITG (RS-ITG) mode based on analytic solutions obtained from a non-local linear calculation in the q min -region.

101.6
In this calculation, analytic solutions are given particularly for the double rational surface and nonresonant RS-ITG modes [14], because these modes are the most unstable slab modes in the reversed shear configuration. It is noted that according to the numerical analyses based on the gyrokinetic integral eigenvalue code, the growth rate of the single rational surface RS-ITG mode (k = 0 at the q min -surface) is much smaller than the above two cases. This is because a wide marginally stable region with v ti v te < |ω/k | appears in the q min -region, and unstable regions of the slab ITG mode, where v ti < |ω/k | v te is satisfied, are shifted into finite magnetic shear regions on both sides of the q min -surface [2].
We consider here a sheared slab geometry, where the x-direction corresponds to the radial direction, the z-direction is chosen in the direction of the magnetic field at the reference magnetic surface x = 0 or r = r s and the y-direction is chosen to be normal to both the xand z-directions. We assume the periodic boundary condition in the yand z-directions, and the fixed boundary condition with conducting walls in the x-direction. By expanding the q-profile around the position x = 0, we write the q-profile as q(x) = q s + q s x + 1 2 q s x 2 + · · ·, where q s , q s and q s are evaluated at x = 0. For the q min -region with q s = 0, we choose the model magnetic configuration as where L s = (2q 2 s R)/(q s r s ), and x = 0 corresponds to the position of the q min -surface. By linearizing the gyrokinetic Vlasov-Maxwell system, equations (1) and (6), we have the second-order ordinary differential equation, where the diamagnetic drift frequency is ω * e = (k y T e /m e Ω e L n ), Z i = Z(ξ i ) is the Fried-Conte plasma dispersion function, ξ i = ω/( √ 2|k |v ti ), τ ≡ T e /T i and ρ s = √ τ ρ ti . Each variable in equation (19) is normalized as follows:Ω ≡ ω/ω * e ; x ≡ x/ρ s ;k y ≡ k y ρ s andφ ≡ eφ/T e . Using the asymptotic expansion for the plasma dispersion function, i − · · ·, under the fluid limit approximation, ξ i > 1, we have a reduced form of the eigenmode equation [10], For the model magnetic field, equation (18), the eigenmode equation is rewritten as

101.7
By considering c 1 as a perturbation parameter, a perturbation theory [11] can be applied to equation (22). The unperturbed solution of equation (22) is obtained as a solution for the standard Weber equation, where l denotes a radial mode number, H l is the lth order Hermite polynomial and the eigenfunction, equation (23), is normalized as |φ l | 2 dx = 1. Solving the perturbed eigenfunctionφ (1) l and the perturbed energy level (1) l yields the eigenfunction, and the dispersion relation, whereφ (0) * l denotes a complex conjugate ofφ (0) l . Equation (25) involves two characteristic solutions depending on the relative sign betweenk z andk y . Whenk z > 0 (andk y > 0), the mode has two mode-rational surfaces on both sides of the q min -surface. Assuming the eigenfrequency satisfying |Re(Ω)| |Im(Ω)|, the potential in equation (22) is recognized as a parabolic potential well with a fourth-order perturbation, and the eigenfunction, equation (25), becomes basically a bounded solution in thex space. On the other hand, whenk z < 0 (andk y > 0), the configuration has no mode-rational surface. For the eigenfrequency with |Re(Ω)| |Im(Ω)|, a potential function in the eigenmode equation behaves as a parabolic potential hill with a fourth order perturbation. In this case, the eigenfunction, equation (25), corresponds to an oscillatory solution. For an unstable mode, the asymptotic solution is written as [12] lim where C is a constant. This solution has the group velocity of the outgoing wave: This result shows a stabilizing effect of shear convective damping as seen in the conventional slab ITG mode [13]. Therefore, we expect that the growth rate of the double rational surface RS-ITG mode is larger than that of the nonresonant RS-ITG mode. These modes are also confirmed from non-local analyses of the slab ITG mode using a gyrokinetic integral eigenvalue code [2]. Figure 2 shows typical eigenfunctions of the double rational surface and nonresonant RS-ITG modes. Corresponding to the analytic solutions, the double rational surface (nonresonant) mode shows a bounded (oscillatory) solution. In [14], it was shown that both modes give much larger growth rates than the slab ITG mode in a normal shear configuration, because a stabilizing effect of magnetic shear disappears in the q min -region. Moreover, the growth rate of the nonresonant mode is slightly smaller than that of the double rational surface mode, because of the shear convective damping. It is noted that in [15], a gap mode structure at the q min -surface was predicted based on a kinetic ballooning theory. However, the slab ITG mode was not included in that analysis.

Gyrokinetic global analysis in toroidal geometry
In this section, global analyses of the ITG mode using GT3D are shown for tokamak-fusiontest-reactor-(TFTR-) like parameters [16]: deuterium plasma, R 0 = 2.6 m, a = 0.94 m and B 0 = 4.6 T. The density and temperature profiles are given by equations (16) and (17) with T i0 = 26 keV, T e0 = 8 keV, L ti /R 0 = 0.11, L te /R 0 = 0.54, L n /R 0 = 0.54, ∆r n = ∆r t = 0.3a and r 0 = 0.5a. It is noted that a destabilizing effect for the ITG mode is largest at r = r 0 = 0.5a in this configuration. In these plasma parameters, a characteristic device size parameter is given as a/ρ ti ∼ 320, and we will discuss properties of the ITG mode in reactor relevant parameters. In the analysis, the normal and reversed shear configurations shown in figure 3 are assumed, and both results are compared to understand properties of the ITG mode. Figure 4 shows eigenfunctions of the ITG mode with n = 20. In the normal shear configuration, the mode shows a typical toroidal mode structure that balloons on the weak field side. However, the mode structure observed in the reversed shear configuration is quite different from the normal shear case. In the q min -region, the mode shows a slablike mode structure which also has a finite amplitude on the strong field side. This slablike mode has a large radial correlation length. On the outside of the q min -surface, where the magnetic shear is positive, the mode structure shows a toroidal feature. In the reversed shear configuration, these two types of ITG mode are coupled to form a single global eigenfunction. The m-spectra observed at the q = 2 (r = 0.5a) surface or the q min -surface are shown in figure 5, where sample points are 101.9  In the reversed shear configuration, the ITG mode shows a slablike mode structure which also has a finite amplitude on the strong field side. chosen in the straight-field-line coordinates. In the normal shear configuration, the m-spectrum has a sharp peak at m = nq with small side band components which come from a toroidal mode coupling. The m-spectrum in the reversed shear configuration also peaks at m = nq. However, the spectrum contains much larger side band components than the normal shear case. The side band components consist of a strong double rational surface mode with m = nq + 1 and a weak nonresonant mode with m = nq − 1. Since the distance between the neighbouring mode rational surfaces becomes large in a weak magnetic shear region, the toroidal mode coupling is weak in the q min -region. Therefore, it is concluded that these double rational surface and nonresonant modes show a contribution from the RS-ITG modes shown in section 3. It is noted that as for a resonant mode with m = nq, the contribution from the slab ITG mode is very weak, because it has a wide marginally stable region with v ti v te < |ω/k | around the q min -surface. Thus, it is considered that the resonant component mainly shows a contribution from the toroidal ITG mode. Based on local theory, the threshold condition for the slab ITG mode is given by [17] where b = k 2 ⊥ ρ 2 ti . On the other hand, the threshold condition for the toroidal ITG mode is written as [18] where n = L n /R 0 . For the flat density case, the toroidal ITG mode is characterized by L ti and condition (34) becomes where ti = L ti /R 0 . Thus, we expect that the coupled mode structure in the reversed shear configuration is changed from a slab mode to a toroidal mode, depending on η i and L ti . The η i dependence of the growth rate is shown in figure 6. In this parameter study, density and temperature profiles are chosen so that the average temperature and L ti are fixed. In the normal shear configuration, where the toroidal ITG mode is dominant, the growth rate is almost constant for η i = 1.5-5. On the other hand, in the reversed shear case, the growth rate goes down on decreasing η i . This is because the contribution from the RS-ITG mode, which has a critical value around η i ∼ 1.5 [2], becomes weak at a small η i parameter. This is also seen in the eigenfunction in figure 7. In the normal shear configuration, the mode structure is not affected by varying η i . However, in the reversed shear configuration, a slablike mode structure disappears at a small η i parameter. Figure 8 shows the L ti -dependence of the growth rate. In figure 9, the growth rates at each L ti are obtained without changing the average temperature or η i . In both the normal and reversed shear configurations, the growth rates become small as the temperature (and density) gradient is decreased. However, their mode structures are qualitatively different. In the normal shear Figure 6. The η i dependence of the growth rate is plotted for the ITG mode with n = 20. In this study, L n is changed and L ti , L te and the average temperature are fixed. In the reversed shear configuration, the growth rate goes down with decreasing η i , and the slablike mode disappears. show the normal and reversed shear cases, respectively. In the reversed shear case, a slablike mode structure around the q min -region disappears (cf figure 3(b)). configuration, the mode keeps a typical toroidal mode structure for all L ti parameters. On the other hand, in the reversed shear configuration, a slablike mode structure becomes distinct at a large L ti -parameter.
All these results show that compared with the normal shear configuration, where the toroidal ITG mode is dominant, a contribution from the slab ITG mode becomes relatively important in the reversed shear configuration.

Summary
The toroidal ITG mode has been believed to be a main candidate for the anomalous ion thermal transport in high temperature tokamak plasmas. However, in this paper, we have shown that a contribution from the slab ITG mode also becomes important in the q min -region of reversed  In the reversed shear case, a toroidal mode structure becomes weak, and the contribution from the RS-ITG mode becomes dominant (cf figure 3(b)).
From gyrokinetic global analyses using a newly developed toroidal particle code, it is found that the ITG mode in reversed shear tokamaks is characterized by a coupled mode structure between the slab and toroidal ITG modes. Moreover, especially in the q min -region, a slablike feature is clearly seen in the eigenfunction. As for the slab ITG mode, the RS-ITG mode, which was found from nonlocal gyrokinetic analyses in a sheared slab geometry, is also identified in toroidal calculations. As predicted from local theory, the strengths of contributions from the slab and toroidal ITG modes depend on η i and L ti , respectively. When η i is decreased without changing L ti , the slablike feature becomes weak and the growth rate goes down only in the reversed shear configuration, while the toroidal ITG mode in the normal shear configuration is not so affected. On the other hand, when L ti is increased without changing η i , the contribution from 101.13 the slab ITG mode becomes relatively strong compared with that from the toroidal ITG mode.
In [15], a gap mode structure across the q min surface was reported. It seems that there is some inconsistency between [15] and the present paper because, in the present analyses, such a gap mode structure has not been found as a linear eigenfunction. However, in [15], the gap mode structure is observed only in a case where the maximum temperature gradient region is localized inside the q min -surface, and this structure consists of two independent modes, which are excited in normal and reversed shear regions on both sides of the q min -surface. Thus, the situation discussed in [15] is a quasilinear saturation phase of several linearly independent modes with close growth rates, and problems treated in both papers are essentially different. In order to check the existence of the gap mode structure in the quasilinear or fully nonlinear phase of the ITG mode, it is important to study nonlinear properties of above qualitatively different ITG modes in reversed shear tokamaks. In the future work, nonlinear simulations will also be addressed using GT3D.