Simulation of ion temperature gradient driven modes with 6D kinetic Vlasov code

With the increase in computational capabilities over the last years it becomes possible to simulate more and more complex and accurate physical models. Gyrokinetic theory has been introduced in the 1960s and 1970s in the need of describing a plasma with more accurate models than fluid equations, but eliminating the complexity of the fast gyration about the magnetic field lines. Although results from current gyrokinetic computer simulations are in fair agreement with experimental results in core physics, crucial assumptions made in the derivation make it unreliable in regimes of higher fluctuations and stronger gradient, such as the tokamak edge. With our novel optimized and scalable semi-Lagrangian solver we are able to simulate ion-temperature gradient modes with the 6D kinetic model including the turbulent saturation. After thoroughly testing our simulation code against analytical computations and gyrokinetic simulations (with the gyrokinetic code GYRO), it has been possible to show first plasma properties that go beyond standard gyrokinetic simulations. This includes the explicit description of the complete perpendicular energy fluxes and the excitation of high frequency waves (around the Larmor frequency) in the nonlinear saturation phase.


I. Introduction
Kinetic models are capable of describing physical phenomena in a tokamak plasma on all scales from the size of the device (∼1m) down to microscales of the order of the electron Larmor radius (∼ 10 −4 m).However, until recent years, the computational capabilities did not allow the computation of the full 6D-kinetic equation for a magnetized plasma with the time resolution required to resolve the fast orbital motion of the articles around the magnetic field.To circumvent this problem, many gyrokinetic models have been developed [10], which reduce the dimensionality and the required time resolution by eliminating frequencies of the order of the Larmor frequency   .Gyrokinetic simulations codes such as GENE [6], GS2 [9], GYRO [2] or ORB5 [7] are in good agreement with experiments in the core of fusion devices where only small perturbation amplitudes and gradients (in density and temperature) are present.Nonetheless, in regimes of large gradients and high fluctuation amplitudes, such as the plasma edge of a tokamak, the gyrokinetic approximations are debatable and at least those models based on a   approximation break down completely.Moreover, the limited frequency range in gyrokinetic simulations precludes the accurate representation of physical processes associated with high-frequency modes.To illustrate this point, Craddock et al. demonstrated the suppression of turbulence by ion Bernstein waves in their study [3].With the rise of modern high-performance computing capabilities, it becomes more and more viable to simulate the full 6D kinetic equation.For this purpose, the massively parallel semi-Lagrangian code BSL6D for the Vlasov equation in 6D phase space has been developed [8].The code, BSL6D, is capable of simulating plasma turbulence across a broad range * Electronic address: Mario.Raeth@ipp.mpg.de of frequencies, including those beyond the Larmor frequency, and with arbitrarily large fluctuation amplitudes.Past efforts to venture into regimes beyond gyrokinetic theory include the drift-cyclotron model by Waltz and Deng [18] and a non-gyrokinetic magnetized plasma turbulence code developed by Sturdevant et al. [11,16,17].However, they either still use a reduced kinetic model, or are limited to a single toroidal model, effectively resulting in a 5D simulation.This work shows unprecedented 6D nonlinear kinetic simulations of the slab ITG instability, examining the energy fluxes inherent in the model.In Section II, we provide a concise introduction to the physical model implemented in the code for slab Ion Temperature Gradient (ITG) simulations.Subsequently, we conduct a verification study in the linear regime, encompassing a comparison with the analytically derived dispersion relation (section III 2) and the computation of quasi-linear energy fluxes (section IV 1).Finally, in section V, we present results from simulations capturing the nonlinear saturation phase of the ITG instability, comparing them to results obtained from the gyrokinetic code CGYRO.Intriguingly, the 6D kinetic simulations reveal that the saturation phase leads to the excitation of high-frequency ion Bernstein waves.We conducted our nonlinear studies employing both a Boussinesq approximation for the gradients (corresponding to the gyrokinetic description of the gradients) and a fully nonlinear treatment.This allows us to investigate the saturation phase under different conditions and gain a deeper understanding of the underlying physics.The discovery of high-frequency ion Bernstein waves during the saturation phase challenges conventional understanding and opens new avenues for research into the nonlinear dynamics of plasma turbulence.

II. 6D semi-Lagrangian kinetic turbulence code for magnetized plasmas
The Vlasov equation describes the motion of a plasma in presence of electromagnetic fields.We consider ions in a constant and homogeneous magnetic field B = ẑ represented by a 6D distribution function  , with an electric field E originating from the interaction with adiabatic electrons.The target of our 6D simulations is the ion kinetic equation in dimensionless variables (  (Ion Larmor radius) =  th (Ion thermal velocity) =  (Ion background density) =  (Ion background temperature) = 1) For the simulation of ion of temperature gradient (ITG) driven modes, the gradient is introduced in Boussinesq approximation.For this purpose, we assume that a background distribution function  0 exists, which fulfills This condition is met, when the background distribution is parameterized by an arbitrary function is the location of the gyrocenter of a given particle, where r is the configurations space coordinate and ρ = ẑ × v is the Larmor radius vector).The background distribution is chosen to resemble a Maxwellian with varying temperature  (R), and thus reads When splitting the distribution function in the Vlasov equation into a perturbation   =  −  0 and the background  0 , one obtains where the gradient on the right-hand side is given by The limit that the gradient length is infinite  (R) →  = 1 is taken, such that the background distribution is independent of R. Thus,  0 is replaced with the homogeneous Maxwellian For simplicity, we assume that the electrons are adiabatic   =  −    ≈   −    and the plasma is quasineutral   =  (with the electron density   ), the complete system of equation is given by For all presented results, the electron temperature has been chosen to be the same as the ion temperature   = 1.The system is integrated in time using a semi-Lagrangian method [15], which combines the advantages of grid-based (Eulerian) and particle-based (Lagrangian) methods.In this method, the distribution function is represented as a discrete function on a grid.However, instead of explicitly computing partial derivatives (as in Eulerian methods), the distribution function is advected along the characteristics of the Vlasov equation.This is achieved by solving the characteristic equations for particle trajectories and interpolating the distribution function at their new positions.This approach circumvents the strict CFL (Courant-Friedrichs-Lewy) condition imposed on Eulerian methods, allowing for the selection of a larger time step.Additionally, the noise levels are significantly lower compared to particle methods, making the semi-Lagrangian method a more robust and efficient choice for simulating the dynamics of plasmas.Implementation details can be found in the work of Kormann et al. [8].

III. Dispersion relation 1. Analytical description
The response of the non-adiabatic distribution, given by ℎ =  −   +    , to a present electrostatic potential  in the linearized system is given by the Gordeev integral [5].This computation yields an expression for the velocity distribution (v = (  ,   ,   )  ) of a specific Fourier mode of the distribution function ℎ  , influenced by an electrostatic wave  =     (k•r−  ) with frequency  and wave vector k (where where   () represents the Bessel function of the first kind, and v * is defined in equation (5).The density response To simplify the integration, the source term k • v *   , defined in equation ( 5) can be written in terms of a derivative with respect to the helping variable  (assuming ∇  ∝ ŷ) which is set  = 1 after the derivative is calculated.The parallel velocity integral is resolved using the plasma dispersion function After solving the perpendicular integral, we introduce The density response is introduced in the quasi-neutrality equation which gives the electrostatic potential form equation (7).
The resulting dispersion relation, given by solution (k) to the equation describes all modes for the 6D kinetic system.In the gyrokinetic limit  ≪ 1 and   ≪ 1 all terms of the sum vanish except the  = 0 contribution, simplifying the dispersion relation to

Verification of dispersion relation
For the numerical verification of the code against the dispersion relation, a linear simulation has been conducted.The removal of the non-linear term from the Vlasov equation ( 1) results in indefinite linear growth, simplifying the determination of growth rates.In space, a box with dimensions   = 128 × 128 × 8 has been selected.The box size ensures the smallest perpendicular wave number is  ,0 =  ,0 = 0.3 (  =   = 20 3 ).The parallel length is determined by the wave number of the fastest-growing mode in the system, dependent on the temperature gradient.A temperature gradient of ∇  = 0.05 has been employed, leading to an ideal parallel wave number of   = 1 240 and, consequently,   = 480.The velocity space is symmetric in all directions, with  max = 4 and   = 32 × 32 × 32.For interpolation, a 7th-order Lagrange interpolation has been applied to the velocity directions, and an 8th order for the spatial directions.
The simulation commences with a small white noise density perturbation, and growth rates are determined by fitting the time evolution of various Fourier modes of the density.Figure 1 presents a comparison between the numerical results and the analytical dispersion relation.The code demonstrates excellent agreement for both the growth rate and frequency.However, discrepancies are more pronounced for very small and larger wave numbers.For small   , the growth rates are 0.0000 0.0005 0.0010 0.0015 0.0020 0.0025  13) in solid lines minimal, or even negative, making the numerical computation challenging.As for larger perpendicular wave numbers, the numerical damping, attributed to the inherent diffusion of the Lagrange interpolation, is comparable to the ITG growth rates.Consequently, the code tends to underestimate the growth rate.This discrepancy can be mitigated by either enhancing the resolution or adopting a higher interpolation order [14].

IV. Contributions to the energy flux in the 6D system 1. Analytical description
The total energy density  of the system results from the sum of the kinetic energy of ions and the energy stored in the electrons The time derivative of the total energy is computed using the Vlasov equation ( 1) Applying integration by parts to the second integral yields which farther can be modified to be expressed as a divergence of the entire integral By introducing and utilizing the field equation ( 7), the expression can be succinctly summarized as Apart from Q = 1 2 ∫ v 2  d 3 , an additional energy flux S =  appears which, can be identified as a Poynting like energy flux.The energy flux in the gyrokinetic system is completely described by the E × B heat flux.Thus, in the gyrokinetic limit, the energy flux in the 6D system has to be equal the E × B heat flux.The various contributions to the energy flux 2  d 3  are obtained by computing the time derivative with the Vlasov equation (1).The Lorentz force term of the Vlasov equation is altered by partial integration After introducing the expression for kinetic energy density 2 vv  d 3  and computing the cross product with the magnetic field   Q × ẑ, the Grassmann identify can be applied to obtain an expression for the perpendicular energy flux For low frequency waves  ≪ 1, the time derivative can be neglected   Q ≈ 0 and without boundary contributions, the flux surface ( − −plane, ⟨•⟩ := 1 The overall change in energy density can be written as In conclusion, the mean perpendicular energy flux consists of three contributions In the following, we analytically compute the contributions to the energy flux for a given solution of the dispersion relation. 1.The Poynting flux S =  is given by the particle flux multiplied by the electrostatic potential.The derivation of the particle flux, denoted by  = ∫ v  d 3 , follows the same approach as for the density response in equation ( 12) using the ansatz for the distribution function from equation (8).For the perpendicular flux, computations are carried out independently in the two directions ( 1 ,  2 ) =  ⊥ (− sin , cos ), where  is the azimuth angle in velocity space, ensuring that the angle of ẑ × k corresponds to  =  2 .For a more general calculation, the integral is evaluated for an arbitrary complex Fourier mode   = ∫  ⊥    d 3 .We recall the definition of the non-adiabatic perturbation of the distribution function ℎ =  −   +    , which allows us to compute the integral   from equation ( 8) and subsequently assemble the particle flux .The double sum simplifies when the angle integral is evaluated with ((  −  + ))d = 2 , + which leads to the expression The velocity integrals are solved numerically for a given set of parameters and the respective solution (, k) of the dispersion relation (13).When computing the flux in a spatial direction, the definition of the velocity angle  has to be kept in mind.
The angle is defined such that k • ẑ corresponds to  =  2 .Thus, the fluxes point in the direction perpendicular (for Γ 1 ) and parallel (for Γ 2 ) of the wave vector.In the case of k =  ŷ, Γ 1 and Γ 2 corresponds to the  and  direction respectively.For a general wave vector, the fluxes need to be combined to compute the correct fluxes in  and  direction.The resulting velocity integrals where Ξ =  −   ∇    , are then computed numerically.2 ℎ  d 3  is computed from the distribution function (8) resulting in

The
3. The stress induced energy flux is computed from the stress tensor response   = ∫ vvℎ  d 3 .When the velocity components ( 1 ,  2 ,   ) = (− ⊥ sin ,  ⊥ cos ,   ) are introduced, the expression for the energy flux component reads The energy flux is assembled using the response derived in equation (25).The velocity integrals have been computed numerically for the given parameters.

Verification of energy fluxes
In the numerical testing, the contribution to the energy flux derived in section IV 1 are computed from the numerically determined particle flux  = ∫ v  d 3 , the energy density 2  d 3 , and the components of the stress tensor  = ∫ vv  d 3 .Subsequently, these numerical results are compared to the analytical descriptions for the individual contributions to the total energy flux as derived in equations (23).The analytical calculations involve the introduction of the particle flux (28), energy density (30), and the stress tensor response, as described in equation ( 23).In the simulations detailed in section III 2, where the growth rate and frequency of various modes have been determined, a linearized model is used.In the linearized limit, the velocity advection reduces to ∇ • ∇   ≈ −∇ • v   to remove nonlinear effects.However, while this simplified model accurately reproduces the particle flux, higher moments of the distribution function are not faithfully represented.For these simulations, the full advection of the background distribution by the electric field must be considered.To prevent saturation, a very small initial perturbation amplitude (here,  ∼ 10 V. Non-linear ITG simulations

Comparison with gyrokinetic simulation
Having demonstrated the effective performance of the simulations in the linear regime, we have delved into the investigation of the nonlinear saturation of the ITG instability.For this purpose, a simulation similar to section IV 2 has been conducted.However, the grid size has been adjusted from  = 128 × 128 × 8 to  = 64 × 64 × 64, specifically targeting the necessity for a higher resolution in -direction due to the emergence of high-frequency waves during the nonlinear phase.This adjustment ensures that the various  − −slices remained connected.Figure 3 displays snapshots of the density perturbation in the  −  plane during three distinct phases of the simulation.In the first plot, a distinct mode structure is evident, characteristic of a slab ITG instability where the system is dominated by modes with k ⊥ ∇ (in this case, ∇/|∇ | = x).The growth rates decrease rapidly as the wave number parallel to the temperature gradient increases.The second snapshot depicts the density during the transition from the linear to the nonlinear phase.The third image represents the perturbation in a fully saturated turbulent state.The nonlinear saturation of the simulation is compared to the gyrokinetic code CGYRO [2].To achieve this, a simulation with the same spatial parameters (length and resolution) is conducted.In gyrokinetic codes, the velocity space is not fully resolved as in 6D kinetic simulations but is parameterized by two coordinates: the kinetic energy,  = [0,  max ], and the cosine of the pitch angle,  ∈ [−1, 1].The pitch angle characterizes the ratio between parallel and perpendicular velocity.The chosen parameters for the simulation are  max = 8,   = 12, and   = 32.This simulation is initiated with a random noise density perturbation having the same amplitude as the kinetic simulation.For this comparison, emphasis will be placed on the heat flux induced during both the linear and nonlinear phases.Section IV 1 has established that the sole heat flux present in gyrokinetic theory is the E × B heat flux.Therefore, this serves as the metric for comparison.In Figure 4 (top), the heat flux parallel to the gradient is depicted for the two simulations.Following a linear phase (up to  = 3000), during which the flux experiences exponential growth, the energy flux reaches saturation.The saturation levels are similar for both simulations; however, in both cases, the levels are susceptible to significant fluctuations.This is partially attributed to the small domain size.The smallest wave number in the system is  0 = 0.3, whereas the dominant growing mode has a wave number   = 0.9.Consequently, only three wavelengths of the dominant mode fit into the system.Increasing the domain size would lead to a stronger averaging and thus, a more stable saturation level.In addition to expanding the domain, a longer runtime would enhance the comparison by providing more data and enabling a statistical analysis.However, due to the substantial computational cost associated with the 6D kinetic simulation, this option has been foregone, and the obtained results are considered satisfactory.
For further verification, we conducted a comparison of the ratio between the heat flux and the  2 -norm of the electrostatic potential, as depicted in Figure 4 (bottom).For the dominant mode, with k = ( 9 10 , 0, 1 240 ), a ratio of approximately 0.9 is anticipated, as indicated in Figure 2 (2nd row, left), which is in close agreement with the simulated results.Following the nonlinear phase, the ratio decreases to approximately 0.25, a consistent result for both models.In summary, a very good agreement is observed between the two simulations.
A notable distinction between the two simulations is evident.When examining the system in the turbulent state, one observes the presence of high-frequency oscillations with a frequency closely aligned with the Larmor frequency.The existence of such high-frequency oscillations in the saturation of nonlinear ion temperature gradient simulations has been previously acknowledged in the literature [11].However, a comprehensive description, let alone an investigation into their excitation, is still lacking.Figure 5 displays a spectrogram of the electrostatic potential.The ITG intensity is visible at the bottom of the figures with frequencies close to zero.Modes with frequencies close to the harmonics of the Larmor frequency are excited during the nonlinear saturation phase (at  ∼ 3500).The modes are clustered in frequency bands close to the harmonics with a mean frequency slightly larger than the harmonic.The distribution of the high frequency waves suggest the presence of ion Bernstein waves (IBWs) [2].More precisely neutralized IBWs are electrostatic waves in the ions with the presence of adiabatic electrons.Various mechanisms have been proposed to explain the instability of Ion Bernstein Waves (IBWs) [12,13,19].Our analysis indicates that the excitation is triggered by local negative velocity space gradients of the distribution function, induced by the temperature gradient source term [14].These findings underscore the necessity of a nonlinear treatment of the temperature gradient in the full-f 6D model.

Nonlinear treatment of gradients
This section transitions from a treatment of the gradients in a local limit, using the Boussinesq approximation to a nonlinear approach for handling gradients in our simulation code.The chosen initial condition in the simulation is designed to allow the distribution function to exhibit a density and temperature profile The background profiles (r − ρ) and  (r − ρ) are defined in gyrocenter coordinates R := r − ρ (where ρ := ẑ × v is the Larmor radius vector) to establish a background that does not oscillate with the Larmor frequency.To simplify the treatment of boundary conditions, the profiles are periodically set up using a sine-profile in the -direction To prevent the background density gradient from generating an electric field, all modes with wave numbers parallel to of the perturbation in the linear phase is displayed in figure 6.The perturbation exhibits two peaks situated at the maxima of the normalized gradient  (  )  (  ) .To determine the wavenumber, a sine wave is fitted to the envelope of the unstable within the full width half maximum (FWHM) of the perturbation (in x-direction), resulting in a wavenumber of approximately   ≈ 0.44.To facilitate a comparison between the nonlinear gradient and the analytical calculation, the mean of the gradient is computed across the mode profile.Upon determining the wavenumber and effective gradient, the solution of the dispersion relation is obtained from the analytical expression in equation (13).The anticipated complex frequency for the fastest-growing mode is  = 0.01850 + 0.00420.Figure 7 displays the root-mean-square (RMS) of the electrostatic potential on a logarithmic scale (shown in blue), alongside the expected linear growth rate represented by a dashed gray line.The results indicate that the analytically predicted growth rate is slightly higher than the expected value.The growth rate and frequency for the simulation are derived from the Fouriertransformed electrostatic potential (k, ), resulting in The growth rate and frequency exhibit an error of approximately 2% when compared to the analytically calculated values.
After the linear phase, the simulation gradually reaches saturation around  ≈ 6000 and settles into a fully saturated turbulent state (compare figure 6

VI. Summary
Six-dimensional kinetic simulations, which enable the representation of modes with frequencies around the Larmor frequency (such as ion Bernstein waves), serve as a valuable tool for investigating physics beyond conventional gyrokinetic models.The development of our semi-Lagrangian solver for the Vlasov system marks the initial stride in constructing a comprehensive kinetic simulation code for plasma simulations at all tokamak relevant parameters.
The verification tests conducted demonstrated excellent agreement in the growth rate and frequency for the gyrokinetic modes induced by the Ion Temperature Gradient (ITG) instability, with the values computed from the analytical dispersion relation (13).This consistency extends across a broad range of wave numbers.Additionally, we established that the code accurately reproduces the quasi-linear fluxes throughout the domain.In the course of this, a novel formulation of the energy in the 6D kinetic system has been developed, representing a sum of its distinct components (namely, Poynting flux, E × B energy flux, and stress-induced energy flux).The explicit expression for the energy flux establishes a connection to gyrokinetic theory, as the E × B heat flux can be recognized as one of its components.Additionally, the individual contributions can be analytically computed, facilitating a comparison with the simulation results.Beyond the linear verification of our model, we have conducted simulations extending well into the nonlinear phase.Comparisons with the nonlinear gyrokinetic code CGYRO indicate that the code accurately depicts the nonlinear saturation of a slab ITG instability, yielding consistent saturation amplitudes.Examination of the saturation phase and the ensuing turbulent state has led to the identification of the excitation of high-frequency modes.The gyrokinetic model, by definition, does not account for these excited ion Bernstein waves.The potential impact of these modes on the saturation process and energy transport remains unexplored.Some studies from the 1990s suggested that these waves could suppress turbulence levels in the plasma edge [3,4].Additionally, alongside simulations employing a Boussinesq limit of the gradients, we have demonstrated the reproducibility of ITG simulations with a nonlinear treatment of the gradients.The possession a tool that enables the simulation of such modes presents a valuable opportunity to uncover new physics in regimes beyond commonly utilized models.
and  is the azimuth angle in velocity space with ẑ × k corresponding to  = 2 )[1]

4 -FIG. 1 :
FIG.1: Comparison of the frequency (left) and the growth rate (right) for various modes with   = 1 240 determined from linear simulations (•) and computed by numerical root finding from the analytical dispersion relation(13) in solid lines

1 .
Poynting flux S =  2. E × B -energy flux Q E×B = −(∇ × ẑ) 3. Stress tensor induced energy flux Q  = ( • ∇) × ẑ This depiction of the energy flux offers two significant advantages.Firstly, it provides a direct mean of calculating the individual contributions to the energy flux from the code.In the gyrokinetic limit, the Poynting flux S and the stress-induced heat flux Q  nullify each other, leaving the E × B-heat flux as the sole remaining form of energy transport.The equivalency between these two energy fluxes serves as a useful test for verifying the accuracy of the energy flux in the gyrokinetic regime.Secondly, it facilitates the direct computation of contributions to the energy flux, rendering it less susceptible to gyro-oscillations and, consequently, easier to calculate.
E × B heat flux Q  × ⊥ = −∇ × ẑ can be computed analogously to the Poynting flux.In Fourier space, the flux can be written as Q E×B ⊥ = k × ẑ  , where the energy density response   = ∫  2

FIG. 2 :
FIG. 2: Comparison (left column) of energy fluxes (Poynting flux S, E ×B heat flux Q E×B , and stress induced heat flux Q Π ) opposite to the temperature gradient for various wave numbers between (1) results from numerical simulations (•) and (2) analytical computations (-).Respective deviation from the analytical results are shown in the second column

− 8 )
must be chosen.All other parameters remain consistent with the simulation of the dispersion relation in Section III 2 ( = 128 × 128 × 8 × 32 × 32 × 32,   =   = 20 3 ,   = 480, Δ = 0.03), and the distribution function is initialized with a white noise density perturbation.The results for the contributions of the energy flux are shown in figure 2 in comparison to the analytically obtained results together with the absolute error in the left column.The E × B heat flux is an order of magnitude larger than the other two fluxes.For large wave numbers, the E × B heat flux is proportional to the -component of the wave vector   .The values of the Poynting flux and the stress induced heat flux are nearly identical, as anticipated, given that the gradient has been selected close to the gyrokinetic limit.In gyrokinetic theory, the entire energy flux is determined by the E × B heat flux.The code results demonstrate strong agreement with analytical predictions.Absolute errors are less than 10% for both the pointing flux and the stress-induced heat flux while the relative error of the E × B heat flux is approximately 1%.

FIG. 3 :
FIG. 3: Snapshots of the particle density in the --plane (with fixed ) for various points in time

4 FIG. 7 :
FIG. 7:Root mean square of electrostatic potential compared to expected growth rate