2DESR: a two-dimensional Fourier-space gyrokinetic eigenvalue code for the ion-temperature-gradient modes in tokamaks

A two-dimensional (2D) gyrokinetic eigenvalue solver, 2DESR, has been developed to solve the 2D gyrokinetic eigenvalue problem in the poloidal Fourier space for the ion-temperature-gradient (ITG) modes in tokamaks. With full kinetic effects of ions retained, the 2D gyrokinetic eigenvalue equations in the poloidal Fourier space have been derived and numerically solved in the 2DESR code. In the linear ITG Cyclone test with adiabatic electrons, the 2DESR code benchmarks well against the gyrokinetic initial-value codes GENE and NLT. It is found that two branches of ITG modes coexist in the system.

💡 Research Summary

The paper introduces a novel two‑dimensional gyrokinetic eigenvalue solver, named 2DESR, designed to treat ion‑temperature‑gradient (ITG) modes in tokamak plasmas directly in the poloidal Fourier space. Unlike traditional ballooning‑space eigenvalue codes that rely on translational invariance to reduce the problem to one dimension, 2DESR retains the full (r, θ) dependence by expanding perturbations in the coupled Fourier basis (z = n q(r) − m, m). The authors retain the complete kinetic description of ions while adopting an adiabatic electron model, thereby solving the Vlasov–Poisson system without the simplifying assumptions (e.g., well‑circulating particles) often employed in existing 2D eigenvalue codes.

The theoretical development starts from the linear gyrokinetic Vlasov equation and the quasineutrality condition, applies a Laplace transform in time and Fourier transforms in toroidal angle ζ and the modified poloidal angle χ, and derives the coupled equations (8) and (9) for the non‑adiabatic ion distribution g_i and the electrostatic potential φ̂ in Fourier space. The gyro‑average operator is expressed through an inverse Fourier integral, and the coupling among neighboring poloidal harmonics is truncated to |p| ≤ 3 after convergence tests.

Numerically, the domain (z, m, v_∥, μ) is discretized with uniform grids in z and v_∥, Gauss‑Laguerre quadrature (N_μ = 16) for μ, and a finite‑difference stencil for the derivative operators. The Vlasov equation becomes a linear system A(ω) g = B(ω) φ̂, where A is an N³ × N³ sparse matrix (N³ = N_z × N_m × N_v). Using the sparse direct solver PARDISO, the authors compute C(ω) = A⁻¹B and eliminate g, yielding the nonlinear eigenvalue problem M(ω) φ̂ = 0. Newton’s method, consistent with previous eigenvalue solvers (ESR), iterates on the complex frequency ω until convergence. Parallelization is performed over the μ dimension with MPI; for a typical run (N_z = 81, N_m = 22, N_v = 32, N_μ = 16) the matrix size is 57 024 and the total wall‑clock time on 64 Intel Xeon Gold cores is about 140 s.

Benchmarking uses the standard Cyclone base case (circular flux surfaces, q(r) profile given in Eq. 13, identical ion/electron temperature and density gradients). The toroidal mode number n is scanned from 5 to 45. 2DESR identifies two distinct unstable ITG branches: “mode 1” peaks near n ≈ 20, while “mode 2” peaks near n ≈ 30. Growth‑rate and real‑frequency curves from 2DESR agree with the initial‑value gyrokinetic codes GENE and NLT, but each code captures only one of the branches (GENE matches mode 1, NLT matches the most unstable branch, which switches from mode 1 to mode 2 around n ≈ 35). The proximity of the two growth‑rate curves near n ≈ 35 explains why initial‑value simulations struggle to resolve the dominant mode without very long integration times (≈ 2 500 R₀/C for n = 35 versus ≈ 60 R₀/C for n = 20).

Eigenmode structures are examined for n = 5 and n = 20. The electrostatic potential exhibits the characteristic up‑down asymmetry of tilted ballooning modes. For n = 20 the mode is radially localized between r/a ≈ 0.3 and 0.6, with maximum amplitude at r/a ≈ 0.45. For the lower‑n case (n = 5) the mode spreads over a broader radial interval and shows noticeable differences among adjacent poloidal harmonics. These structures are in quantitative agreement with those obtained from the initial‑value code NLT, confirming the physical fidelity of the eigenvalue approach.

A discussion on coordinate choice emphasizes why the (z, m) Fourier representation is advantageous for capturing the radial localization of each poloidal harmonic on its rational surface, while an (E, μ) representation (used in earlier ESR work) is less suitable for trapped ions because the trapped‑particle distribution vanishes outside a limited bounce angle. By retaining (v_∥, μ) as velocity coordinates, the authors can treat the full poloidal domain (θ ∈