First principles gyrokinetic analysis of electromagnetic plasma instabilities

TL;DR

This paper presents a gyrokinetic theoretical framework and numerical analysis of electromagnetic plasma instabilities in tokamaks, including energy conservation and destabilizing mechanisms, using the Orb5 code.

Contribution

It provides a first-principles derivation of the gyrokinetic model, energy conservation diagnostics, and analysis of instability mechanisms in electromagnetic tokamak plasmas.

Findings

Energy conservation verified via Noether theorem in Orb5

Analysis of destabilizing mechanisms for ITG and electrostatic instabilities

Transition to Kinetic Ballooning Modes at high electromagnetic beta

Abstract

A two-fold analysis of electromagnetic core tokamak instabilities in the framework of the gyrokinetic theory is presented. First principle theoretical foundations of the gyrokinetic theory are used to explain and justify the numerical results obtained with the global electromagnetic particle-in-cell code Orb5 whose model is derived from the Lagrangian formalism. The energy conservation law corresponding to the Orb5 model is derived from the Noether theorem and implemented in the code as a diagnostics for energy balance and conservation verification. An additional Noether theorem based diagnostics is implemented in order to analyse destabilising mechanisms for the electrostatic and the electromagnetic Ion Temperature Gradient (ITG) instabilities in the core region of the tokamak. The transition towards the Kinetic Ballooning Modes (KBM) at high electromagnetic $\beta$ is also investigated.

Tables

Table 1. Table 1: List of benchmark parameters

$r_0/a$	$0.5$
$a/L_{\mathrm{ref}}$	$0.36$
$R_0/L_{\mathrm{ref}}$	$1.0$
$T_i\left(r_0\right)/T_{\mathrm{ref}}=T_e\left(r_0\right)/T_{\mathrm{ref}}$	$1.0$
${\kappa }_{T_i}={\kappa }_{T_e}$	$6.96$
$w_{T_i}=w_{T_e}$	$0.3$
$n_i\left(r_0\right)/T_{\mathrm{ref}}=n_e\left(r_0\right)/T_{\mathrm{ref}}$	$1.0$
${\kappa }_{n_i}={\kappa }_{n_e}$	$2.23$
$w_{n_i}=w_{n_e}$	$0.3$
$m_i/m_{\mathrm{ref}}$	$1.0$
$m_e/m_{\mathrm{ref}}$	$5.44617\cdot {10}^{-4}$

Table 2. Table 2: Nominal reference and derived reference values based on the low elongation magnetic surfaces case (CBC) [ 21 ] , Fig.5, discharge # ​ 81499 # 81499 \#81499 , at time t = 4000 𝑡 4000 t=4000 ms and ρ = 0.5 𝜌 0.5 \rho=0.5 , which after the rescaling of magnetic surfaces shape towards concentric surfaces corresponds to r / a = 0.5 𝑟 𝑎 0.5 r/a=0.5 .

$m_{\mathrm{ref}}(=m_D)/m_p$	$2.0$
$n_{\mathrm{ref}}(=n_e)/{10}^{19}{m}^{-3}$	$4.66$
$T_{\mathrm{ref}}(=T_e)/\mathrm{KeV}$	$2.14$
$B_{\mathrm{ref}}(=B_t\left(R_{\mathrm{mag}}\right))/T$	$2.0$
$L_{\mathrm{ref}}(=R_0=R_{\mathrm{mag}})/m$	$1.67$
${\beta }_{\mathrm{ref}}$	$0.0045$
${\rho }^{\ast }={\rho }_s/a$	$0.00555\sim 1/180.2$

Table 3. Table 3: Data corresponding to analysis of linear electromagnetic simulations provided on Fig. 2.

$\beta$	${\omega }_{\mathrm{Orb5}}$	${\gamma }_{\mathrm{Orb5}}$	$\Delta E_F$
in $\%$	in $c_s/R_0$	in $c_s/R_0$
${10}^{-5}$	$0.667$	$0.533$	$0.999$
${10}^{-4}$	$0.661$	$0.527$	$0.999$
$0.025$	$0.665$	$0.540$	$0.794$
$0.05$	$0.613$	$0.544$	$0.313$
$0.055$	$0.613$	$0.544$	$0.182$
$0.06$	$0.597$	$0.511$	$0.038$
$0.065$	$0.630$	$0.486$	$-0.117$
$0.068$	$0.620$	$0.555$	$-0.217$
$0.07$	$0.677$	$0.516$	$-0.284$
$0.075$	$0.677$	$0.516$	$-0.4649$
$0.1$	$0.651$	$0.497$	$-1.576$
$0.15$	$0.710$	$0.477$	$-4.967$
$0.25$	$0.773$	$0.464$	$-17.870$
$0.45$	$0.741$	$0.391$	$-86.200$
$0.5$	$0.850$	$0.345$	$-117.734$
$0.6$	$0.910$	$0.270$	$-205.114$
$0.65$	$1.260$	$0.225$	$-256.218$
$0.68$	$3.255$	$0.079$	$-283.435$
$0.689$	$3.189$	$0.064$	$-291.325$
$0.69$	$3.456$	$0.0067$	$-293.262$
$0.7$	$3.205$	$0.152$	$-292.720$
$0.75$	$2.870$	$0.196$	$-197.360$
$0.8$	$2.701$	$0.438$	$-202.060$
$0.9$	$2.516$	$0.809$	$-208.500$
$1.000$	$2.365$	$1.083$	$-206.237$
$1.125$	$2.171$	$1.337$	$-204.897$
$1.25$	$2.191$	$1.532$	$-205.668$

Equations

$${\mathcal{A}}^{\text{mbox}\text{sc}Orb5}$$

$$L_s(\mathbf{Z},\dot{\mathbf{Z}})=\frac{q_s}{c}{\mathbf{A}}_s^{\ast }\cdot \dot{\mathbf{X}}+\frac{m_{s}c}{q_s}\mu \dot{\theta }-\left(H_{0,s}+{ϵ}_{\delta }{H}_{1,s}\right),$$

$${\mathbf{A}}_s^{\ast }=\mathbf{A}+\frac{c}{q_s}{p}_z\hat{\mathbf{b}}.$$

$${\mathbf{B}}_s^{\ast }=\nabla \times {\mathbf{A}}_s^{\ast }\text{mbox}and{B}_{\Vert ,s}^{\ast }=\hat{\mathbf{b}}\cdot {\mathbf{B}}_s^{\ast }.$$

$$f_{0,s}({\mathbf{z}}_0)=f_s(\mathbf{Z}({\mathbf{z}}_0,t);t).$$

$$H_{0,s}=\frac{p_z^2}{2m_s}+\mu B.$$

$$H_{1,s}=q_s⟨{\varphi }_1(\mathbf{X}+{\mathbf{\rho }}_{0,s})-A_{1\Vert }(\mathbf{X}+{\mathbf{\rho }}_{0,s})\frac{p_z}{m_{s}c}⟩,$$

$$⟨\varphi (\mathbf{X}+{\mathbf{\rho }}_{0,s})⟩=\frac{1}{2\pi }{\int }_0^{2\pi }\mathrm{d}\theta \varphi (\mathbf{X}+{\mathbf{\rho }}_{0,s}).$$

$$H_{1,e}=e\left({\varphi }_1(\mathbf{X})-A_{1\Vert }(\mathbf{X})\frac{p_z}{m_{e}c}\right).$$

$$H_{2,s}=-\frac{m_s{c}^2}{2B^2}{\left|{\nabla }_{\perp }{\varphi }_1\right|}^2+\frac{q_s^2}{2m_s{c}^2}\left(A_{1\Vert }^2+\left(\frac{m_s{c}^2}{q_s^2}\right)\frac{\mu }{B}{A}_{1\Vert }{\nabla }_{\perp }^2{A}_{1\Vert }\right).$$

$$H_{2,e}$$

$$\delta {\mathcal{A}}^{\text{mbox}\text{sc}Orb5}=\delta {\int }_{t_0}^{t_1}\mathcal{L}\mathrm{d}t=0.$$

$$\int {\mathrm{d}}^n\mathbf{r}\frac{\delta \mathcal{L}}{\delta \eta }\cdot \chi (\mathbf{r})$$

$$\sum _{\mathrm{s}\ne e}{q}_s\int \mathrm{d}{\mathbf{z}}_0{f}_{0,s}({\mathbf{z}}_0)⟨{\hat{\varphi }}_1⟩+e\int \mathrm{d}{\mathbf{z}}_0{f}_{0,e}({\mathbf{z}}_0){\hat{\varphi }}_1$$

$$={ϵ}_{\delta }\sum _{\mathrm{s}\ne e}\int \mathrm{d}{\mathbf{z}}_0{f}_{C,s}(\mathbf{Z}({\mathbf{z}}_0,t))\frac{m_s{c}^2}{B^2}{\nabla }_{\perp }{\varphi }_1\cdot {\nabla }_{\perp }{\hat{\varphi }}_1.$$

$$\sum _{\mathrm{s}\ne e}{q}_s\int \mathrm{d}\mathbf{z}{f}_s(\mathbf{z})⟨{\hat{\varphi }}_1⟩+e\int \mathrm{d}\mathbf{z}{f}_e(\mathbf{z}){\hat{\varphi }}_1$$

$$={ϵ}_{\delta }\sum _{\mathrm{s}\ne e}\int \mathrm{d}\mathbf{z}{f}_{C,s}(\mathbf{z})\frac{m_s{c}^2}{B^2}{\nabla }_{\perp }{\varphi }_1\cdot {\nabla }_{\perp }{\hat{\varphi }}_1,$$

$${ϵ}_{\delta }\int \frac{{\mathrm{d}}^3\mathbf{x}}{4\pi }{\nabla }_{\perp }{A}_{1\Vert }\cdot {\nabla }_{\perp }{\hat{A}}_{1\Vert }=\sum _{\mathrm{s}\ne e}\int \mathrm{d}\mathbf{z}{f}_s(\mathbf{z})\frac{q_s{p}_z}{m_{s}c}⟨{\hat{A}}_{1\Vert }⟩+\int \mathrm{d}\mathbf{z}{f}_e(\mathbf{z})\frac{e{p}_z}{m_{e}c}{\hat{A}}_{1\Vert }$$

$$+{ϵ}_{\delta }\sum _{\mathrm{s}\ne e}\int \mathrm{d}\mathbf{z}{f}_{C,s}(\mathbf{z})\left(\frac{q_s^2}{m_s{c}^2}{A}_{1\Vert }{\hat{A}}_{1\Vert }+\frac{\mu }{B}\left[A_{1\Vert }{\nabla }_{\perp }^2{\hat{A}}_{1\Vert }+{\hat{A}}_{1\Vert }{\nabla }_{\perp }^2{A}_{1\Vert }\right]\right)$$

$$+{ϵ}_{\delta }\int \mathrm{d}\mathbf{z}{f}_{C,e}(\mathbf{z})\frac{e^2}{m_e{c}^2}{A}_{1\Vert }{\hat{A}}_{1\Vert }.$$

$$\int \mathrm{d}{\mathbf{z}}_0\frac{\delta \mathcal{L}}{\delta \mathbf{Z}}\cdot \hat{\mathbf{Z}}=\sum _{\mathrm{s}}\int \mathrm{d}{\mathbf{z}}_0{f}_{0,s}({\mathbf{z}}_0)\left(\frac{\partial L_s}{\partial \mathbf{Z}}\cdot \hat{\mathbf{Z}}+\frac{\partial L_s}{\partial \dot{\mathbf{Z}}}\cdot \dot{\hat{\mathbf{Z}}}\right)$$

$$=\sum _{\mathrm{s}}\int \mathrm{d}{\mathbf{z}}_0{f}_{0,s}({\mathbf{z}}_0)\left(\frac{\partial L_s}{\partial \mathbf{Z}}\cdot \hat{\mathbf{Z}}-\frac{\mathrm{d}}{\mathrm{dt}}\frac{\partial L_s}{\partial \dot{\mathbf{Z}}}\cdot \hat{\mathbf{Z}}+\frac{\mathrm{d}}{\mathrm{dt}}\left[\frac{\partial L_s}{\partial \dot{\mathbf{Z}}}\cdot \hat{\mathbf{Z}}\right]\right).$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial L_s}{\partial \dot{\mathbf{Z}}}=\frac{\partial L_s}{\partial \mathbf{Z}}.$$

$$0=\frac{\mathrm{d}{f}_s}{\mathrm{dt}}=\frac{\partial f_s}{\partial t}+\dot{\mathbf{X}}\cdot \nabla f_s+\dot{p_z}\frac{\partial f_s}{\partial p_z},$$

$$\dot{\mathbf{X}}$$

$${\dot{p}}_z$$

$$\dot{\mathbf{X}}$$

$${\dot{p}}_z$$

$$\dot{\mathbf{X}}$$

$${\dot{p}}_z$$

$$\dot{\mathbf{X}}{|}_0$$

$${\dot{p}}_z{|}_0$$

$$\frac{\mathrm{d}\mathcal{L}}{\mathrm{d}t}$$

$$\int \mathrm{d}{\mathbf{z}}_0\frac{\delta \mathcal{L}}{\delta \mathbf{Z}}\cdot \hat{\mathbf{Z}}$$

$$0=\int \mathrm{d}\mathbf{x}\frac{\delta \mathcal{L}}{\delta {\varphi }_1}\cdot \frac{\partial {\varphi }_1}{\partial t}=\int \mathrm{d}\mathbf{x}\frac{\delta {\mathcal{L}}_s}{\delta A_{1\Vert }}\cdot \frac{\partial A_{1\Vert }}{\partial t}.$$

Full text

N Tronko1,2, A Bottino1, C Chandre3, E Sonnendrücker1,2,

E Lanti4, N Ohana4, S Brunner4, and L Villard4

1 Max-Planck Institute for Plasma Physics, 85748, Garching, Germany

2 TU Munich, Mathematics Center, 85747, Garching, Germany

3 Aix Marseille Univ, CNRS, Centrale Marseille, I2M, Marseille, France

4 Swiss Plasma Center, Ecole Polytechnique Fédérale de Lausanne, CH-1015, Lausanne, Switzerland

1 Introduction

Strongly magnetised fusion plasmas represent a paradigmatic example of out-of-equilibrium systems, in which turbulence is ubiquitous. This omnipresence originates from the concept of magnetic fusion itself: Bringing the mix of hydrogen isotopes into the confinement mode implies by construction the existence of strong spatial gradients. Typically, a difference of three orders of magnitude for the temperature is present, with $\sim 10^{8}K$ at the centre of the device, where the plasma is hot and relatively dense and with $\sim 10^{5}K$ close to the edge, where the plasma is more rarefied and colder. The strong gradients or space inhomogeneities of the temperature, velocity and current are intrinsic to fusion plasmas, and provide sources of free energy and therefore represent sources for instabilities. The instabilities manifest themselves via the exponentially growing perturbations of the electromagnetic fields. When the critical amount of instabilities has been developed, the system moves to a turbulent state with strongly unpredictable behaviour in space and time. In turn, the transport associated with turbulence is extremely dangerous for the plasma confinement. The understanding of its origins is a subject of numerous investigations [1].

Extensive studies in the framework of different formalisms, from the single fluid MHD model and the multi-fluid approaches to kinetic models, aiming to identify the instability mechanisms, have been carried out over the past decades, both analytically and numerically. A detailed overview of these studies is presented in [2], where the transition between the electrostatic and the electromagnetic regimes is also discussed. In particular, in low $\beta=p/[B^{2}/(8\pi)]$ plasmas, where $p$ is the kinetic pressure and $B^{2}/(8\pi)$ the magnetic pressure, the ion temperature gradient instability, known as the ITG mode, has almost an electrostatic polarization, and therefore carries almost only electrostatic energy for the coupled drift waves and the ion acoustic waves that gives rise to a collective instability.

In higher pressure plasmas with $\beta>m_{e}/m_{i}$, the electromagnetic energy is injected, and the ITG mode couples with the shear Alfvén wave becoming a dispersive oscillation with an electromagnetic polarization. As the plasma pressure increases the inductive electric field from the fluctuating magnetic field $\delta B_{\perp}$ begins to cancel part of the electrostatic component of the parallel electric field. This cancellation reduces the energy transfer rate $\left\langle j_{\|}E_{\|}\right\rangle$ and reduces the growth rate of instability.

There are numerous studies of these weakly electromagnetic ITG modes including [3], [4], [5] and [6] which detail how the ITG modes change with increasing $\beta$. At certain intermediate values of $\beta$, both the ITG and the kinetic ballooning modes are present with different frequencies and growth rates.

In this paper, we investigate the transition between the low $\beta$ and finite $\beta$ regimes in the framework of the gyrokinetic theory, both analytically and numerically.

Strongly magnetized plasma exhibits multi-scaled dynamics: The fast rotation around magnetic field lines, called gyromotion, is at least three orders of magnitude faster than the slow drifts across the magnetic field lines. The gyrokinetic dynamical reduction [7],[8] aims to simplify the dynamical description in which fast gyration is systematically and reversibly eliminated, resulting in considerable simplifications and a gain of computational time.

Nowadays, the gyrokinetic (GK) codes play a significant role in the understanding of the development and the saturation of turbulence and the prediction of the subsequent transport properties [9]. Electrostatic gyrokinetic simulations have been the topic of numerous studies during the last decades [10], [11], [12], so that the properties of the electrostatic instabilities in the framework of the gyrokinetic theory are rather well known. However, global electromagnetic simulations are more recent [13], [14], [15] and some elements related to the GK electromagnetic instability mechanisms still need to be clarified.

In order to provide a better understanding of the global electromagnetic GK simulations, we present an analysis of the instability mechanisms by performing global linear electromagnetic simulations with the particle-in-cell code Orb5. The numerical set-up is similar to the one used for benchmarking the global electromagnetic codes [15]. The Orb5 code is based on the GK Lagrangian model [16], [14] and allows one to use diagnostics tools issued from first principles, exactly corresponding to the theoretical model.

This article is organised as follows: In Sec. 2, the field Lagrangian model including all approximations implemented in the code Orb5 is presented. Section 3 provides the derivation of the Orb5 energy invariant through the Noether theorem. In Sec. 4, the diagnostics issued from Noether’s method are derived, and their implementation in the code Orb5 is discussed. Finally, in Sec. 5, the analysis of the electromagnetic instability mechanisms is presented.

2 Variational formulation of the Orb5  code model

A detailed derivation of the Orb5 model in the framework of the Eulerian variational principle [17] is given in [16], where CGS units are used. In this work, we use the Lagrangian variational formulation accordingly to [18], since the Lagrangian formulation of the gyrokinetic field theory lends itself to a discretisation by finite element methods. We notice that the choice of formalism, Eulerian or Lagrangian, does not affect the final expressions for the gyrokinetic Maxwell-Vlasov equations.

The Lagrangian action functional for the code Orb5 depends on the fluctuating electromagnetic potentials $(\phi_{1},A_{1\|})$ and on the coordinates of the Lagrangian particle trajectories $\mathbf{Z}(\mathbf{z}_{0},t)=\left[\mathbf{X}(\mathbf{z}_{0},t),p_{z}(\mathbf{z}_{0},t),\mu(\mathbf{z}_{0},t),\theta(\mathbf{z}_{0},t)\right]$, where ${\bf X}$ are the positions of the gyrocentres, $p_{z}$ their momenta parallel to the magnetic field lines, $\mu$ their magnetic moments and $\theta$ their gyroangles. These coordinates are labelled by their initial conditions $\mathbf{z}_{0}$ such that ${\bf Z}({\bf z}_{0},0)={\bf z}_{0}$. The expression of the action is given by [16]:

[TABLE]

where the particle Lagrangian is given by:

[TABLE]

the generalized vector potential is defined as:

[TABLE]

The volume element in the reduced phase space is ${\rm d}\mathbf{z}_{0}=mB_{\|,s}^{*}({\bf z}_{0}){\rm d}\mathbf{X}_{0}{\rm d}p_{z,0}{\rm d}\mu_{0}{\rm d}\theta_{0}$ (see Ref. [19] for more detail) where

[TABLE]

The dynamical distribution function is denoted by $f_{s}$ for each species, and the background distribution function by $f_{C,s}$. The property of the Vlasov distribution function being conserved along the particles trajectories translates into

[TABLE]

We assume that the background distribution $f_{C,s}$ is frozen and corresponds to a good approximation of $f_{s}$ at all times. The first and the second terms of the action (1) correspond to the gyrocentre reduction, and the last term is a contribution from the perturbed magnetic field. We notice that only the perpendicular part of the perturbed magnetic field $\mathbf{B}_{\perp}=\widehat{\mathbf{b}}\times\bm{\nabla}A_{1\|}$ is implemented.

The background Hamiltonian contains information on the kinetic energy of a charged particle moving in a magnetic field of amplitude $B$:

[TABLE]

The first order correction of the Hamiltonian model for the ions is given by

[TABLE]

where the gyroaveraging operator $\langle\cdot\rangle$ is the average over the fast gyroangle $\theta$ contained in the fast rotating Larmor vector ${\bm{\rho}}_{0,s}$ measuring the distance between the initial particle position and the guiding-centre position. More precisely, in an orthonormal basis $(\hat{\bf b}_{1},\hat{\bf b}_{2},\hat{\bf b})$ at ${\bf X}$, the Larmor vector is given by ${\bm{\rho}}_{0,s}=mc/q_{s}\sqrt{2\mu/m_{s}B}[\hat{\bf b}_{1}({\bf X})\cos\theta-\hat{\bf b}_{2}({\bf X})\sin\theta]$, and the expression of the gyroaveraging is given by

[TABLE]

For the electrons, the first order correction to the Hamiltonian only contains the first order Finite Larmor Radius (FLR) correction. It corresponds to the drift-kinetic dynamics, and it is given by

[TABLE]

In what follows the bracketed quantities are evaluated at the position $\mathbf{X}+\bm{\rho}_{0,s}$, all the other quantities are evaluated at the gyrocentre position $\mathbf{X}$.

There exist several nonlinear models in the Orb5 code. The most complete model considers the nonlinear Hamiltonian model $H_{2}$ for the ions, including all order FLR corrections in its electrostatic part and up to second order FLR terms in its electromagnetic part. In this work, in addition to the drift-kinetic model for the electrons, we consider the long-wavelength approximation of the nonlinear model for the ions (see Eq. (68) in Ref. [16]), i.e. :

[TABLE]

The second order Hamiltonian model for the electrons only contains the drift-kinetic correction:

[TABLE]

Before presenting the equations of motion implemented in Orb5, we discuss all necessary approximations included in the gyrokinetic action given by Eq. (1). The first term of the action involves the full distribution function $f_{s}$, while the second term involving the nonlinear Hamiltonian $H_{2}$ contains a “canonical” distribution function $f_{C,s}$, which is by definition invariant under the unperturbed Hamiltonian dynamics, i.e., it satisfies the condition $\{f_{C,s},H_{0,s}\}_{\mathrm{gc}}=0$, where the guiding-center Poisson bracket is defined accordingly to Eq. (7) in Ref. [16]. This approximation brings several simplifications to the model. First, it results in the linearisation of the gyrokinetic Poisson and Ampère equations. Second, it simplifies the gyrokinetic Vlasov equation by excluding some nonlinear terms from the gyrocentre characteristics associated with Hamiltonian $H_{2}$.

2.1 Gyrokinetic field equation

The equations of motion are derived from

[TABLE]

We use functional derivatives for evaluating the r.h.s. of this expression explicitly. As a reminder, for a functional $\mathcal{L}[\eta]=\int{\rm d}^{n}{\bf r}\ {\mathsf{L}}\left(\eta,\bm{\nabla}\eta\right)$ depending functionally on a scalar field $\eta$ and its gradient $\bm{\nabla}\eta$, the functional derivative is defined as $\delta{\mathcal{L}}/\delta\eta$ acting on the test function $\chi$ as:

[TABLE]

The corresponding quasineutrality equation in the weak form with the test function $\widehat{\phi}_{1}$ is obtained by calculating the functional derivative $\int{\rm d}^{3}\mathbf{x}\ \left(\delta\mathcal{L}/\delta{\phi_{1}}\right)\cdot\widehat{\phi}_{1}$:

[TABLE]

We perform the change of variables $\mathbf{z}=\mathbf{Z}(\mathbf{z}_{0},t)$, and the quasineutrality equation becomes

[TABLE]

where $\mathrm{d}{\mathbf{z}}=mB_{\parallel,s}{\rm d}\mathbf{X}{\rm d}p_{z}{\rm d}\mu{\rm d}\theta$. The Ampère equation obtained from the same variational principle is derived from the computation of $\left(\delta\mathcal{L}/\delta{A_{1\|}}\right)\cdot\widehat{A}_{1\|}$ and the change of variables $\mathbf{z}=\mathbf{Z}(\mathbf{z}_{0},t)$:

[TABLE]

2.2 Nonlinear gyrokinetic Vlasov equation

We now proceed with derivation of the particles dynamics implemented in the Orb5 code. The equations of motion for the particles are obtained by setting to zero the functional derivatives with respect to the phase space coordinates $\mathbf{Z}=\left(\mathbf{X},p_{z},\mu,\theta\right)$:

[TABLE]

Since in the action functional (1) the nonlinear part of reduced particle dynamics, i.e., Hamiltonian $H_{2}$ is only coupled to the non-dynamical part of the distribution function, i.e., $f_{C}$, $H_{2}$ does not contribute to the particle characteristics used for reconstructing the gyrokinetic Vlasov equation. The last term in Eq. (16) vanishes when integrating the Lagrangian in time to get the action integral. However, this term is used later in order to determine the conserved energy of the system by Noether’s theorem. As a consequence, the functional derivatives vanish for all test functions $\widehat{{\mathbf{Z}}}$ if and only if the Euler-Lagrange equation for the particles is satisfied:

[TABLE]

The gyrokinetic Vlasov equation is reconstructed from the linearised gyrocentre characteristics according to the approximations performed on the action functional given by Eq. (1).

[TABLE]

where the linearised gyrocentre characteristics depend on the linearised Hamiltonian model:

[TABLE]

where $H_{s}=H_{0,s}+\epsilon_{\delta}H_{1,s}$. For the ordering considered above, the equations for the ion characteristics become:

[TABLE]

while the equations for the electron characteristics contain the first FLR corrections only:

[TABLE]

We notice that the equations for the unperturbed characteristics for both species coincide:

[TABLE]

3 Noether theorem for the Orb5 code model

In order to derive the expression for the energy, we calculate the time derivative of the Lagrangian density $\mathcal{L}=\mathcal{L}(\mathbf{X},p_{z},\theta,\mu,\phi,A_{\|})$:

[TABLE]

Using the Euler-Lagrange equations given by Eq. (16), we get:

[TABLE]

For the field equations, we choose the test function $\widehat{\phi}_{1}=\partial{\phi}_{1}/\partial t$ and $\widehat{A}_{1\|}=\partial A_{1\|}/\partial t$ and we use the corresponding Euler-Lagrange equations in the weak form:

[TABLE]

Finally, using the fact that the total time derivative of the Vlasov density vanishes, we get the expression for the energy invariant:

[TABLE]

where

[TABLE]

Using the same change of variables ${\mathbf{z}}=\mathbf{Z}(\mathbf{z}_{0};t)$, as in the Poisson and Ampère equations, we get:

[TABLE]

The procedure allowing one to get the power balance diagnostics is the following one: First, we directly substitute the expression for Hamiltonians $H_{0,s}$, $H_{1,s}$ given by Eqs. (6-9) and $H_{2}$ given by Eqs. (10-11). Then we define the unperturbed kinetic energy of the particles:

[TABLE]

and the remaining terms are referred to as field energy:

[TABLE]

Next, the nonlinear term containing $H_{2}$ in the expression for the energy is rewritten using the corresponding quasineutrality and Ampère equations in the weak form. This is achieved by choosing a particular test function $\widehat{\phi}_{1}=\phi_{1}$ and by substituting it in Eq. (14). Similarly, the test function $\widehat{A}_{1\|}=A_{1\|}$ is substituted to the corresponding Ampère equation given by Eq. (15). The quasineutrality equation (14) with $\widehat{\phi}_{1}=\phi_{1}$ is written as:

[TABLE]

The Ampère equation (15) with $\widehat{A}_{1\|}=A_{1\|}$ is written as:

[TABLE]

Now using Eqs. (29-30) we substitute the expressions for the electrostatic and electromagnetic contributions into Eq. (3) and we get Eq. (31):

[TABLE]

For clarity, we define a function for each component of $\mathcal{E}_{\mathrm{EM}}$:

[TABLE]

We remark that these expressions are general for all electromagnetic models and is independent from the choice of the nonlinear model, i.e., the second order Hamiltonian $H_{2}$.

4 Energetically consistent diagnostics for the Orb5 code

The derivation of dynamical invariants via the Noether’s method is naturally included in the Lagrangian framework. It gives an opportunity to construct code diagnostics, allowing one to control the quality of the simulations and to get information about the mechanisms triggering the instabilities.

In particle-in-cell codes, the dynamics of particles and fields is computed in two different ways: Particles are advanced along their characteristics without the use of any grid, while fields are evaluated on a grid. Within one calculation cycle, both sides are communicating: Particle are pushed along their renewed characteristics by using the values of the electromagnetic fields evaluated on the grid. Then the new values of the particle positions are deposed on the grid in order to provide the new values for the charge and current density entering into the electromagnetic field equations (14)-(15).

Considering the energy exchange between particles and fields, i.e., independently calculating the time derivatives of $\mathcal{E}_{\mathrm{kin}}$ and $\mathcal{E}_{F}$, allows one to control the consistency of the algorithm and the quality of the simulation by verifying the energy conservation. Moreover, further application of Noether theorem makes it possible to analytically calculate a simplified expression for the time derivative of the field energy $\mathcal{E}_{F}$, given by the non-perturbed characteristics of the particles only. Such a simplification gives a possibility to access the underlying instability mechanisms through the particle characteristics.

In this section, we provide the detailed derivation of the Orb5 diagnostics developed from the field-particles energy balance equation:

[TABLE]

where the time derivative of the l.h.s. can be evaluated through the particle characteristics and the r.h.s. from the field contributions evaluated on the grid.

Two diagnostics issued from Eq. (36) are defined: First, the power balance diagnostics is defined as the energy balance equation (36), normalised by the field energy $\mathcal{E}_{F}$. Second, the $\Delta\mathcal{E}_{F}$ diagnostics is defined as the energy balance equation (36) normalised by an electrostatic component of the field energy, i.e., $\mathcal{E}_{\mathrm{es}}$.

Using the definition of the kinetic part of the energy $\mathcal{E}_{\mathrm{kin}}$, given by Eq. (35), we explicitly calculate the contributions to Eq. (36):

[TABLE]

where

[TABLE]

The direct implementation of the linear gyrocentre characteristics for $\mathbf{X}$ and $p_{z}$ given by Eqs. (2.2-2.2) leads to the cancellation of all nonlinear terms related to the perturbed electromagnetic fields, i.e., the final expression contains only the contributions corresponding to the unperturbed Hamiltonian dynamics given by Eq. (2.2):

[TABLE]

where

[TABLE]

The geometric contribution $\bm{\nabla}\times\widehat{\mathbf{b}}$ to ${\mathbf{B}}_{s}^{*}$ given by Eqs. (4, 3) is expressed by using the projection on the parallel and perpendicular directions, following the calculations given in Appendix B of [16]:

[TABLE]

where the scalar $\tau$ represents the magnetic twist and the vector $\mathbf{G}$ is referred to as the magnetic curvature. Since $\mathbf{B}\times\left(\bm{\nabla}\times{\mathbf{B}}\right)=-\bm{\nabla}p$ in single fluid MHD equilibrium, we rewrite the curvature vector $\mathbf{G}$ in order to evidence the pressure-like contributions into the characteristics:

[TABLE]

We also decompose the geometric magnetic field in the parallel and perpendicular components in the following way:

[TABLE]

The unperturbed characteristics in the power balance equation are given by:

[TABLE]

where we have used the divergence free property of magnetic field: $\widehat{\mathbf{b}}\cdot\bm{\nabla}B=-B\ \bm{\nabla}\cdot\widehat{\mathbf{b}}$.

4.1 The power balance diagnostics

In order to understand and analyse possible sources of plasma deconfinement, one aims to investigate mechanisms, triggering the growth of microinstabilities and turbulent transport. The mechanisms contributing to the development of microinstabilities are directly related to the exponential growth of electromagnetic field fluctuations. Considering that electromagnetic instabilities have an exponential growth: $\mathcal{E}_{F}=\bar{\mathcal{E}}_{F}e^{\gamma t}$, we derive the expression for the power balance diagnostics:

[TABLE]

Therefore, for practical reasons, in numerical simulations, it is useful to consider the power balance equation in the following form (i.e., normalized by the field energy $\mathcal{E}_{F}$):

[TABLE]

The power balance diagnostics is suitable for quality verification in linear and nonlinear simulations. In addition to that in the case of the linear simulations, the power balance equation not only gives an indication about the quality of the simulation but also can be used to measure the instantaneous growth rate of instability [16].

4.2 $\Delta{\mathcal{E}}_{F}$ diagnostics

In the case of electrostatic simulations, the power balance diagnostics is sufficient for investigating the stabilizing and destabilizing mechanisms. The situation is slightly different in the case of electromagnetic simulations when $\mathcal{E}_{\mathrm{es}}=\mathcal{E}_{\mathrm{em}}$ therefore, $\mathcal{E}_{F}=0$ and the diagnostics defined by Eq. (42) is not defined. In order to investigate the transition between electrostatic and electromagnetic instabilities, we introduce the following diagnostics:

[TABLE]

where functions ${\mathcal{E}}_{\mathrm{es}}$ and ${\mathcal{E}}_{\mathrm{em}}$ are defined accordingly to Eqs. (33)-(34).

This diagnostics allows one to investigate the properties of the electromagnetic simulations from different viewpoints: First, the sign of the function $\Delta{\mathcal{E}_{F}}$ determines if the instability is electrostatic (positive) or electromagnetic (negative). Moreover, $\Delta{\mathcal{E}}_{F}$ allows one to access all stabilizing/destabilizing mechanisms through Eq. (38) even in the situation with $\mathcal{E}_{F}=0$. In addition, investigating functional properties of $\Delta{\mathcal{E}}_{F}$ as a function of magnetic $\beta$, i.e., the zeros and extremum points (minimum for example) gives the possibility to analyse bifurcations and exchanges of stabilizing/destabilizing mechanisms. We provide this analysis for linear electromagnetic simulations in the next section.

5 Analysis of electromagnetic simulations

5.1 Numerical setup

The parameters used in simulations are derived from the Cyclone Base Case (CBC), which is a well established set of parameters for the flux tube (simulations for one magnetic field line) and global studies (simulations covering the full radial range in the small section of the device). First, it has been used [20] for the benchmark of different flux tubes codes and recently for the benchmark of global electromagnetic codes [15]. The original discharge (H- mode shot $\#81499$ taken at $t=4000$ ms and minor radius $r=0.5a$, where $a$ is the the minor radius) of the DIII-D device which serves as a basis for the CBC has naturally more complex shaped flux surfaces. In our case, the equilibrium magnetic configuration is circular and concentric with the inverse aspect ratio $a/R_{0}=0.36$ and the safety factor profile:

[TABLE]

Here $a$ is the minor tokamak radius and $R_{0}$ is a major one, $r$ is the local radius of a flux surface. The temperature and density profiles and their normalized logarithmic gradients are given by:

[TABLE]

which gives us a peaked gradient profile of density and temperature $A=(n,T)$ centred at $r=r_{0}$ with maximal amplitude $\kappa_{A}$ and characteristic width $w_{A}$. The macroscopic reference length $L_{\mathrm{ref}}$ is fixed to the major radius $R_{0}$ in what follows.

The values of parameters used for the benchmark are summarized in Table 1. We remark that the profiles for ions and electrons are chosen to be identical.

The nominal reference values issued from the original experimental work [21] are given in Table  2.

In order to reduce the resolution requirement and the computational effort, the ion-electron mass ratio is set to the proton-electron mass ratio, i.e. the electrons are considered being twice heavier than in reality. Concerning the spatial resolution, the associated finite-size parameter $\rho^{*}=\rho_{i}/a$, defined as the ratio between the ion gyroradius $\rho_{i}$ and the minor radius $a$, is set to $1/180$. In fact, considering a hydrogen aspect ratio for ions would need smaller spatial scales by a factor 2.

The nominal $\beta$ value at the reference position $r/a=0.5$ is close to $0.5\%$ assuming we have taken in consideration the following definition of $\beta_{\mathrm{ref}}=8\pi n_{\mathrm{ref}}T_{\mathrm{ref}}/(B^{2}_{\mathrm{ref}})$ and the normalisation of $n_{\rm{ref}}$, corresponding to the Orb5 code. Here $n_{\mathrm{ref}}$ and $T_{\mathrm{ref}}$ are the density and temperature taken at the reference position $r/a$.

5.2 Electrostatic and electromagnetic instability analysis

Introducing the electromagnetic effects into the gyrokinetic simulations adds significant complexity compared to the electrostatic simulations, and requires a more detailed analysis for the implementation of diagnostics. It has not been remarked in electrostatic simulations that the normalization of the power balance diagnostics (41) on the field energy $\mathcal{E}_{F}$ together with using this diagnostics for the growth rate calculation in linear simulations may introduce some inconsistencies. It starts to be evident for the electromagnetic simulations when the amount of magnetic energy is equal to the amount of electrostatic energy, i.e., $\mathcal{E}_{F}=0$. It is evident that in this case, Eq. (41) cannot be used for measuring the growth rate. However, thanks to a small modification, this diagnostic can be adapted for investigating the behaviour of the instability triggering mechanisms at this transitional point [see Eq. (43)].

The contributions to the growth rate $\gamma$ arising from the different terms in the unperturbed guiding-center characteristics $\dot{\mathbf{X}}|_{0}$ and $\dot{p}_{z}|_{0}$ can be separated in the power balance equation and give a clear vision of which type of instability is present in the system: This diagnostics is suitable for both linear and nonlinear electromangetic simulations:

[TABLE]

where

[TABLE]

Consistently with the choice of magnetic equilibrium effects, we are neglecting the $\bm{\nabla}P$ contribution to the curvature or $\bm{\nabla}B$ drift. Therefore, $\bm{v}_{{\bm{\nabla}P},s}=0$ and $\bm{v}_{\mathrm{curv},s}=\mu B\bm{\nabla}\cdot\widehat{\mathbf{b}}$. This gives us the possibility to follow easily the dynamics of the remaining contributions to ${\rm d}\mathcal{E}_{F}/{\rm dt}$. Accordingly to the sign, the contributions $\bm{v}_{\|,s}$, $\bm{v}_{\rm{curv},s}$ and $\bm{v}_{{\bm{\nabla}B},s}$ to the time derivative of the field energy are considered stabilizing when it is negative or destabilizing when it is positive.

Following the electromagnetic $\beta$ scan of the linear electromagnetic simulations, summarized in the Table 3, we focus on the sign of $\bm{v}_{\|,s}$,$\bm{v}_{\rm{curv},s}$ and $\bm{v}_{\bm{\nabla}B,s}$ contributions. In Fig. 1, different cases of the instability triggering mechanisms are presented. We notice that for the case with $\mathcal{E}_{F}=0$ or close to zero, additional standard linear fit diagnostics for the growth rate is used to avoid numerical errors due to the division by a small number in the denominator of Eq. (41). At the same time, for each value of $\beta$, we monitor the value of $\Delta E_{F}$, defined by Eq. (43), as the value of the electromagnetic field energy normalized by the electrostatic energy (see Fig. 2). The change of sign for $\Delta E_{F}$ corresponds to the transition from the electrostatic to the electromagnetic regime.

In Fig. 1, the ${\rm d}\mathcal{E}_{F}/{\rm d}t$ diagnostics is presented for different values of $\beta$. Figure 1a) represents the electrostatic ITG with $\beta=0.00001\%$, $\Delta E_{F}=0.99$ which is mainly destabilized by curvature contribution with $\bm{v}_{\bm{\nabla}B,s}>0$. Figure 1b): $\beta$ is increased up to $0.0025\%$, $\Delta E_{F}=0.79$, which activates an additional destabilizing mechanism with $\left\langle A_{1\|}\right\rangle\bm{v}_{\rm{curv},s}>0$. Figure 1c): $\beta=0.06\%$ corresponds to $\Delta E_{F}=0.03$, i.e., the amount of electrostatic energy is almost equal to the amount of electromagnetic energy in the system: This corresponds to the exchange of the stabilizing and destabilizing mechanisms, i.e., the mode is now destabilized by kinetic effects with $\bm{v}_{\|,s}>0$. However, the growth rate and the frequency of the mode exhibit no bifurcation. Figure 1d): $\beta=0.5\%$ corresponds to the electromagnetic ITG mode with $\Delta E_{F}=-86.20$, destabilized by kinetic effects $\bm{v}_{\|,s}>0$. Figure 1e): $\beta=0.69\%$ with $\Delta E_{F}=-296.76$, a minimal value for the normalized field energy corresponds to the ITG to KBM bifurcation in frequencies and the growth rates. Figure 1f): $\beta=1.125\%$ $\Delta E_{F}=-204.89$ corresponds to the high frequency KBM mode with $\bm{v}_{\|,s}>0$.

In Fig. 2, we perform the $\beta$-scan with both energy conservation based diagnostics, the power balance diagnostics given by Eq. (41) and the $\Delta{\mathcal{E}}_{F}$ diagnostics. In Fig. 2a), the growth rate of the instabilities as a function of $\beta$ is calculated according to the power balance diagnostics given by Eq. (41). The minimum corresponds to the bifurcation from ITG to KBM. On Fig. 2b), the corresponding frequencies are presented: the transition between the slow (ITG) and the fast mode (KBM) happens at the value $\beta=0.65\%$, which corresponds to the minimum of the power balance diagnostics. Figures 2c) and 2d) represent the $\Delta{\mathcal{E}}_{F}$ diagnostics. The functional dependency on $\beta$ is investigated: On Fig. 2c), the zero of the $\Delta{\mathcal{E}}_{F}(\beta)$ function corresponds to the transition from the electrostatic ITG to the electromagnetic ITG, which corresponds to the exchange of stabilizing/destabilizing mechanisms, from the curvature drift $\bm{v}_{\rm curv,s}$ to the kinetic mechanisms driven by $\bm{v}_{\|,s}$. The minimum of the $\Delta{\mathcal{E}}_{F}(\beta)$ function corresponds to the bifurcation between the ITG and KBM modes, i.e., the minimum of the power balance diagnostics from Fig. 2a). Figure 2d) represents the part of Fig. 2c), with low $\beta$, featuring zero of $\Delta\mathcal{E}_{F}(\beta)$ function.

6 Conclusions

In this work, the energy exchange channels leading to destabilisation of the electromagnetic instabilities in global gyrokinetic simulations with Orb5 code have been identified. First, it has been observed through the simulations that with increasing $\beta$ the contribution of the magnetic curvature in the ITG mode destabilisation decreases together with the mode growth rate. It confirms the results of previous studies, indicating the stabilizing role of the magnetic $\beta$ on the ITG instability. Second, the implementation of the diagnostics, issued from Noether’s method, allows one to investigate the transition from the electrostatic towards the electromagnetic ITG regime. Following the contributions to the energy time variation, the transition from the destabilizing role of curvature drift towards its stabilizing role is identified. The implementation of the energy invariant based diagnostic $\Delta{\mathcal{E}}_{F}$ allowed to systematically analyse essential features of the electromagnetic instabilities, for instance, by looking at the zeros of $\Delta{\mathcal{E}}_{F}$ corresponding to the ITG instability mechanism exchange and the minimum of that function corresponding to the bifurcation between the ITG and the KBM mode. We observe that the transition of the instabilities generating mechanisms in the framework of the gyrokinetic theory follows the results obtained in previous studies realised in the framework of fluid and kinetic approaches [2]. This transition happens at $\beta=0.06\%$, which follows previous theoretical predictions for $\beta>m_{e}/m_{i}\sim 0.055\%$. Finally, it has been identified that the mode bifurcation towards the kinetic ballooning mode dominated regime happens at $\beta=0.65\%$, i.e., at a $\beta$ value $10$ times higher than the $\beta$ value at which mixing between the ITG (slow) and KBM (fast) modes occurs together with the destabilizing mechanisms exchange.

Acknowledgements

The first author thanks P.J. Morisson and J.D. Meiss for useful advices and discussions. This material is based upon work supported by the National Science Foundation under Grant No. DMS-1440140 while NT and CC were in residence at the Mathematical Sciences Research Institute in Berkeley, California, during the Fall 2018 semester. This work has been carried out within the framework of the French Federation for Magnetic Fusion Studies (FR-FCM) and of the Eurofusion consortium, and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement No.  633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.