Field theory and structure-preserving geometric particle-in-cell algorithm for drift wave instability and turbulence

TL;DR

This paper introduces a structure-preserving geometric Particle-In-Cell algorithm for simulating drift wave instability and turbulence in magnetized plasmas, ensuring long-term accuracy and conservation properties.

Contribution

It develops a novel geometric PIC algorithm based on discrete exterior calculus and Hamiltonian splitting that preserves symplectic structure and gauge symmetry.

Findings

The algorithm accurately simulates ion Bernstein and drift waves.

Simulation results show turbulence energy diffusion scales between Bohm and gyro-Bohm.

Density blobs form as prominent structures in fully developed turbulence.

Abstract

A field theory and the associated structure-preserving geometric Particle-In-Cell (PIC) algorithm are developed to study low frequency electrostatic perturbations with fully kinetic ions and adiabatic electrons in magnetized plasmas. The algorithm is constructed by geometrically discretizing the field theory using discrete exterior calculus, high-order Whitney interpolation forms, and non-canonical Hamiltonian splitting method. The discretization preserves the non-canonical symplectic structure of the particle-field system, as well as the electromagnetic gauge symmetry. As a result, the algorithm is charge-conserving and possesses long-term conservation properties. Because drift wave turbulence and anomalous transport intrinsically involve multi time-scales, simulation studies using fully kinetic particle demand algorithms with long-term accuracy and fidelity. The structure-preserving geometric PIC algorithm developed adequately servers this purpose. The algorithm has been implemented in the \textsl{SymPIC} code, tested and benchmarked using the examples of ion Bernstein waves and drift waves. We apply the algorithm to study the Ion Temperature Gradient (ITG) instability and turbulence in a 2D slab geometry. Simulation results show that at the early stage of the turbulence, the energy diffusion is between the Bohm scaling and gyro-Bohm scaling. At later time, the observed diffusion is closer to the gyro-Bohm scaling, and density blobs generated by the rupture of unstable modes are the prominent structures of the fully developed ITG turbulence.

Equations

$$-\frac{q_i}{q_e}{n}_i=n_e=n_{e0}\exp \left(-\frac{q_e\varphi }{T_e}\right)\text{leavevmode},$$

$${\ddot{\mathbf{x}}}_p$$

$$S\left[{\mathbf{x}}_p,\varphi \right]$$

$$\frac{\delta S}{\delta {\mathbf{x}}_p}$$

$$\frac{\delta S}{\delta \varphi }$$

$$S^{\prime }\left[{\mathbf{x}}_p\right]$$

$$\rho \left(\mathbf{x}\right)$$

$$\varphi \left(\rho \right)$$

$$\frac{\delta S^{\prime }}{\delta {\mathbf{x}}_p}$$

$$S_{sd}\left[{\mathbf{x}}_p,{\varphi }_I\right]=\int \mathrm{d}t{L}_{sd}\left[{\mathbf{x}}_p,{\varphi }_I\right]\text{leavevmode},$$

$$L_{sd}\left[{\mathbf{x}}_p,{\varphi }_I\right]$$

$$\frac{\delta S_{sd}}{\delta {\mathbf{x}}_p}$$

$$\frac{\delta S_{sd}}{\delta {\varphi }_I}$$

$${\varphi }_I$$

$${\rho }_I$$

$${\ddot{\mathbf{x}}}_p$$

$${\mathbf{E}}_0$$

$${\mathbf{B}}_0$$

$${\ddot{\mathbf{x}}}_p$$

$${\mathbf{E}}_J=-\sum _I{{\nabla }_{\mathrm{d}}}_{J,I}{\varphi }_I\text{leavevmode}.$$

$$L_{sd}^{\prime }\left[{\mathbf{x}}_p\right]$$

$$V$$

$${\rho }_I$$

$$H_p\left({\mathbf{x}}_p,{\mathbf{v}}_p,W_p,t_p\right)=\frac{1}{2}{m}_i{\mathbf{v}}_p^2-W_p.$$

$$\displaystyle\left\{F,G\right\}_{p}=\nabla_{p}F\left[\begin{array}[]{cccc}0&\frac{1}{m_{i}}I&0&0\\ -\frac{1}{m_{i}}I&\frac{q_{i}}{m_{i}^{2}}\hat{B}_{0}\left(\mathbf{x}_{p},t_{p}\right)&\left(\frac{q_{i}}{m_{i}}\frac{\partial\mathbf{A}_{0}\left(\mathbf{x}_{p},t_{p}\right)}{\partial t_{p}}\right)^{T}&0\\ 0&-\frac{q_{i}}{m_{i}}\frac{\partial\mathbf{A}_{0}\left(\mathbf{x}_{p},t_{p}\right)}{\partial t_{p}}&0&-1\\ 0&0&1&0\end{array}\right]\left(\nabla_{p}G\right)^{T}\leavevmode\nobreak\ ,$$

$${\nabla }_{p}F$$

$${\hat{B}}_0$$

$$H$$

$$\left\{F,G\right\}$$

$$\dot{f}$$

$$H=H_x+H_y+H_z+H_W+H_V\text{leavevmode},$$

$$\dot{f}$$

$$\displaystyle\left\{\begin{array}[]{ccl}\dot{\mathbf{x}}_{p}&=&0\leavevmode\nobreak\ ,\\ \dot{\mathbf{v}}_{p}&=&\partial V/\partial\mathbf{x}_{p}\leavevmode\nobreak\ ,\\ \dot{W}_{p}&=&0\leavevmode\nobreak\ ,\\ \dot{t}_{p}&=&0\leavevmode\nobreak\ ,\end{array}\right.$$

$$\displaystyle\Theta_{V}\left(\Delta t\right):\left\{\begin{array}[]{ccl}\mathbf{x}_{p}&\rightarrow&\mathbf{x}_{p}\leavevmode\nobreak\ ,\\ \mathbf{v}_{p}&\rightarrow&\mathbf{v}_{p}+\Delta t\partial V/\partial\mathbf{x}_{p}\leavevmode\nobreak\ ,\\ W_{p}&\rightarrow&W_{p}\leavevmode\nobreak\ ,\\ t_{p}&\rightarrow&t_{p}\leavevmode\nobreak\ ,\end{array}\right.$$

$$\text{For all }p$$

Full text

Structure-preserving geometric algorithm, particle-in-cell, drift wave instability, ion temperature gradient turbulence

1 Introduction

Drift wave turbulence and associated anomalous transport Krall and Rosenbluth (1965); Coppi (1967); Tang (1978); Horton (1999) are important physical processes in magnetic fusion devices. They have been intensively studied using the Particle-In-Cell (PIC) methods Dawson (1983); Hockney and Eastwood (1988); Birdsall and Langdon (1991), which numerically solve the Vlasov-Maxwell or Vlasov-Poisson equations. The dynamics of charged particles in a magnetic field described by the Vlasov equation contains multiple timescales, e.g., the cyclotron frequencies of electrons and ions, plasma frequency, and the drift wave frequency. When simulating low frequency phenomena directly using the PIC method, the time-step must be chosen small enough to resolve the high frequency dynamics of charged particles. Thus, the total number of time-steps required is large, often exceeding computer resource available. The low frequency drift wave instability is such a case, where the ratio between wave frequency and the electron gyro-frequency is in the order of $10^{-5}$. To overcome this difficulty, simplified models which eliminate some of the high-frequency processes while properly describing the slow ion dynamics are developed. A commonly adopted such kinetic model is based on adiabatic electron assumption and quasi-neutrality condition. In this model, for plasmas with one ion species, the ion density $n_{i}$, electrons density $n_{e}$, and the electrostatic potential $\phi$ are linked as

[TABLE]

where $q_{i}$ and $q_{e}$ are the charges of ions and electrons, and $T_{e}$ is the electron temperature. The ion dynamics is governed by Newton’s equation with the Lorentz force,

[TABLE]

where $\mathbf{B}_{0}$ is the background magnetic field, $\mathbf{E}=-\nabla\phi+\mathbf{E}_{0}$ is the electric field, and $\mathbf{E}_{0}$ is the background electric field.

To further decrease computational complexity, gyrokinetic particle simulation methods have been developed and applied to study low frequency instabilities and turbulent transport Lee (1983); Dimits and Lee (1993); Parker and Lee (1993); Parker et al. (1993); Hu and Krommes (1994); Dimits et al. (1996); Lin et al. (1998); Parker et al. (1999). In spite of the success of gyrokinetic simulations, it was pointed out recently that the basic ordering of the gyrokinetic theory Frieman et al. (1966); Catto (1978); Frieman and Chen (1982); Dubin et al. (1983); Hahm (1988); Brizard (1989); Qin et al. (1998, 2000); Qin and Tang (2004); Qin (2005); Qin et al. (2007); Burby (2015); Burby et al. (2015) is not always valid in certain parameter regimes for modern magnetic fusion devices, especially for the H-mode pedestal physics Wan and Parker (2012); Wan et al. (2013) and when density perturbations are large Deng and Waltz (2015). Moreover, due to the requirement of accuracy and numerical stability, the time-step $\Delta t$ in gyrokinetic simulations are often restricted to the same order of ion gyro-period Chen et al. (2008); Chowdhury et al. (2016) already. In these situations, the gyrokinetic method has no significant computational advantage over fully kinetic methods. Recently, a fully kinetic ion scheme was developed Miecnikowski et al. (2018); Sturdevant et al. (2017); Hu et al. (2018). However, even though adiabatic electron model removes the fast electron dynamics from the system, drift wave instabilities and turbulence still evolve in a slow timescale, about one thousandth of ion gyro-period. Simulating these low frequency physics using fully kinetic ions requires a large number of time-steps, and the long-term conservative properties of the numerical schemes become crucial. Conventional PIC methods are based on direct discretization of differential equations, for which numerical errors in general accumulate coherently during the iterations, and long-term simulation results are not reliable.

In the present study, we use a very different approach to construct an explicit high-order structure-preserving geometric PIC algorithm for simulating low frequency drift wave instabilities and turbulence in magnetized plasmas. First, a field theory for low frequency electrostatic dynamics is established with fully kinetic ions and adiabatic electrons. Then the field theory is geometrically discretized using Discrete Exterior Calculus (DEC) Hirani (2003); Desbrun et al. , Whitney interpolating forms Whitney (1957); Squire et al. (2012a, b); Xiao et al. (2015a), and the powerful Hamiltonian splitting method for Vlasov-Maxwell systems Xiao et al. (2015a); He et al. (2015a, 2017). The resulting structure-preserving geometric PIC algorithm is able to preserve the non-canonical symplectic structure associated with the particle-field system, and numerical results show that the simulation error on the energy of the system is bounded by a small number for all time-steps. In addition, the algorithm is gauge independent and thus exactly complies with the discrete local charge conservation law. The knowledge of magnetic potentials is not needed, which is convenient in practical. Furthermore, the algorithm is locally explicit such that it is more efficient on parallel clusters compared with implicit schemes.

In the last ten years, structure-preserving geometric algorithm has become an active research topic in plasma physics. Since 1980s, symplectic integrators for solving Hamiltonian systems have been systematically studied Ruth (1983); Feng (1985, 1986); Feng and Qin (2010); Forest and Ruth (1990); Channell and Scovel (1990); Candy and Rozmus (1991); Hong and Qin (2002); Tang (1993); Shang (1994, 1999); Sanz-Serna and Calvo (1994); Marsden et al. (1998); Sun and Qin (2000); Marsden and West (2001); Hairer et al. (2006). The idea of geometric integrators is to find a discrete one-step iteration map that preserves the symplectic 2-form exactly as the analytical solution of the Hamiltonian system does. According to theoretical and numerical investigations, numerical errors of symplectic integrators on invariants of the systems, such as the total energy and momentum, can be bounded by small numbers for all time-steps Feng (1986); Sanz-Serna and Calvo (1994); Shang (1999); Hairer et al. (2006). In plasma physics, many fundamental models are canonical or non-canonical Hamiltonian systems, and corresponding structure-preserving geometric integrators were recently developed, including those for guiding centers Qin and Guan (2008); Qin et al. (2009); Zhang et al. (2014); Ellison et al. (2015a); Burby and Ellison (2017); Kraus ; Ellison et al. (2018); Ellison (2016), charged particles Qin et al. (2013); He et al. (2015b); Zhang et al. (2015); Ellison et al. (2015b); Zhang et al. (2016); Wang et al. (2016); He et al. (2017, 2016a); Tu et al. (2016); Zhou et al. (2017); Xiao and Qin ; Shi et al. (2019), Vlasov-Maxwell systems Squire et al. (2012a, b); Xiao et al. (2013); Evstatiev and Shadwick (2013); Shadwick et al. (2014); Xiao et al. (2015b, a); Qin et al. (2016); He et al. (2016b); Kraus et al. (2017); Morrison (2017); Xiao et al. (2017, 2018a, 2018b), ideal two-fluid systems Xiao et al. (2016), magnetohydrodynamics Zhou et al. (2014, 2016, 2017); Zhou (2017), Schrödinger-Maxwell system Chen et al. (2017) and Klein-Gorden-Maxwell Shi et al. (2016, 2018); Shi (2018) system. One of the defining characteristics of structure-preserving geometric algorithms is that they are all based on the underpinning field theories and the geometric discretization thereof. Structure-preserving geometric algorithms have demonstrated unparalleled long-term stability and conservative properties compared with conventional non-geometric methods.

The study reported here represents a new development in this research field. We customarily design a field theory for low frequency electrostatic perturbations with fully kinetic ions and adiabatic electrons, and geometrically discrete the field theory to build a structure-preserving geometric PIC algorithm for simulating drift wave instabilities and turbulence in magnetized plasmas.

The paper is organized as follows. In Sec. 2, the field theory for low frequency electrostatic perturbations with fully kinetic ions and adiabatic electrons is established, which is the starting point of our study. Section 3 constructs structure-preserving geometric PIC algorithm by geometrically discretizing the field theory. The algorithm is tested using the examples of Ion Bernstein Waves (IBWs) and drift waves in Sec. 4, and then applied to study the Ion Temperature Gradient (ITG) instability and turbulence in a 2D slab geometry.

2 Field theory for low frequency electrostatic perturbations

To build an effective geometric PIC algorithm for drift wave instabilities and turbulence, a field theory for low frequency electrostatic perturbations is required. With the assumptions of adiabatic electrons and quasi-neutrality condition, the key of establishing the field theory is to find an appropriate action integral whose Euler-Lagrange (EL) equations recover Eqs. (1) and (2). We have found such an action integral. It is

[TABLE]

where $T_{e}$ and $n_{e0}$ are functions of $\mathbf{x}$ and $t$, $\mathbf{A}_{0}$ is the external magnetic potential which gives $\mathbf{B}_{0}=\nabla\times\mathbf{A}_{0}$ and $\mathbf{E}_{0}=-\dot{\mathbf{A}}_{0}$. The system evolves according to the EL equations,

[TABLE]

It can be easily verified that Eqs. (4) and (5) are equivalent to Eqs. (2) and (1).

If we insert $\phi$ obtained from Eq. (5) to the action integral Eq. (3), then the resulting new action integral is

[TABLE]

where

[TABLE]

The action integral $S^{\prime}$ does not depend on $\phi$, and the corresponding EL equation,

[TABLE]

is equivalent to Eqs. (4) and (2).

By design, the field theory is only applicable to low frequency electrostatic perturbations with adiabatic electrons. However, it captures the fully kinetic dynamics of ions, and the corresponding geometric algorithm constructed in the next section possesses long-term accuracy and fidelity, a necessity for fully kinetic particle simulations guaranteed only by the structure-preserving nature of the algorithm.

3 Structure-preserving geometric PIC algorithm

In previous study, we have built explicit geometric algorithms for the Vlasov-Maxwell system and ideal two-fluid system Xiao et al. (2015a, 2016, 2018a). The action integral given by Eq. (3) or (6) is similar to those in previous work. We therefore apply the same techniques of DEC Hirani (2003); Desbrun et al. and high-order Whitney interpolating forms Xiao et al. (2015a, 2016) to perform the spatial discretization. The resulting spatially discretized action integral is

[TABLE]

where

[TABLE]

is the spatially discretized Lagrangian. Here, the subscript $I$ is the grid index, $T_{e,I}$ and $n_{e0,I}$ are electron temperature and density fields on the grid, $W_{\sigma_{0}I}\left(\mathbf{x}\right)$ is the Whitney interpolating map for 0-forms (scalar fields) Xiao et al. (2015a, 2016). The equations of motion for the discrete system are

[TABLE]

Equation (13) plays the role of Poisson’s equation, which links the charge density and the electrostatic potential, i.e.,

[TABLE]

Equation (12) is Newton’s equation with the Lorentz force for the $p$-th particle,

[TABLE]

where

[TABLE]

Using the property of Whitney interpolating map Xiao et al. (2015b), Eq. (16) can be rewritten as

[TABLE]

where

[TABLE]

Akin to the relationship between $S\left[\mathbf{x}_{p},\phi\right]$ and $S^{\prime}\left[\mathbf{x}_{p}\right]$, we can obtain a discrete Lagrangian independent of $\phi_{I}$ by inserting Eq. (14) into Eq. (11),

[TABLE]

From Eq. (21) we see that the system now involves only particles. It is not difficult to build symplectic algorithms for $L_{sd}^{\prime}\left[\mathbf{x}_{p}\right]$ using the techniques of variational integrators Marsden and West (2001); Hairer et al. (2006). However, directly applying these techniques will break the electromagnetic gauge symmetry of the system, which causes charge accumulation and results in numerical stability. To overcome this shortcoming, explicit Hamiltonian splitting method Xiao et al. (2015a); He et al. (2015a, 2017); Zhou et al. (2017) for charged particle dynamics and the Vlasov-Maxwell system have been developed. To solve for $L_{sd}^{\prime}\left[\mathbf{x}_{p}\right]$, here we adopt a similar but more general Hamiltonian splitting method Xiao and Qin .

First, we introduce a non-canonical Hamiltonian structure for the $p$-th charged particle by extending the phase space into 8-dimensional,

[TABLE]

The associated Poisson bracket is

[TABLE]

where

[TABLE]

The Hamiltonian $H$ and Poisson bracket $\left\{\cdot,\cdot\right\}$ for the extended system are

[TABLE]

where $V$ is defined in Eq. (22). Hamilton’s equation is

[TABLE]

Newton’s equation with the Lorentz force (16) is equivalent to $\dot{\mathbf{v}}_{p}=\left\{\mathbf{v}_{p},H\right\}$ and $\dot{\mathbf{x}}_{p}=\left\{\mathbf{x}_{p},H\right\}$.

Next, we split the Hamiltonian $H$ into 5 parts,

[TABLE]

where $H_{W}=-\sum_{p}W_{p}$, $H_{V}=V$, $H_{x}=\sum_{p}m_{i}v_{p,x}^{2}/2$, $H_{y}=\sum_{p}m_{i}v_{p,y}^{2}/2$, and $H_{z}=\sum_{p}m_{i}v_{p,z}^{2}/2$. Each part represents a sub-Hamiltonian system. For example, the equation of motion generated by $H_{V}$ is

[TABLE]

i.e.,

[TABLE]

Its exact solution map $\Theta_{V}\left(\Delta t\right)$ is

[TABLE]

We can exactly solve all other sub-systems$H_{W}$, $H_{x}$, $H_{y}$ and $H_{z}$ in a similar way. The exact solution maps are listed as follows.

[TABLE]

Using these exact solutions of the subsystems, we can compose symplectic iteration schemes of the entire system. For instance, a 1st-order symplectic scheme is

[TABLE]

and a symmetric 2nd-order symplectic scheme can be built as

[TABLE]

A $2(l+1)$-th order symplectic scheme can be constructed from a $2l$-th order symplectic scheme as Yoshida (1990)

[TABLE]

We should point out that there exist a discrete variational approach that generates the same explicit schemes as $\Theta_{1}$, $\Theta_{2}$ and $\Theta_{2l}$. See Refs. Xiao et al. (2018a); Xiao and Qin for details.

4 Simulations of ion Bernstein waves and drift wave instabilities

We have implemented the 2nd-order explicit structure-preserving geometric PIC algorithm given by Eq. (70) in the SymPIC code to simulate low-frequency electrostatic perturbations with fully kinetic ions and adiabatic electrons. The Whitney interpolating maps are chosen to be the same as those in Ref. Xiao et al. (2016). As a benchmark and test, the algorithm is applied to study the ion Bernstein waves. It is then used to simulate the drift wave instability and ion temperature gradient turbulence in a 2D slab geometry.

4.1 Dispersion relation of ion Bernstein waves

To simulate the IBWs in a homogeneous magnetized plasma, the follow system parameters are chosen. External magnetic field is in the $z$-direction, $\mathbf{B}_{0}=B_{0,z}\mathbf{e}_{z}$ with $B_{0,z}=2.5$T. Plasma density $n_{i,0}=1\times 10^{18}\mathrm{m}^{-3}$, and the thermal velocity of ions $v_{T,i,0}=7.23\times 10^{-4}\mathrm{c}$, where $\mathrm{c}$ is the speed of light in the vacuum. The mass and charge of ions are $m_{i}=3.342\times 10^{-27}$kg and $q_{i}=1.6\times 10^{-19}$C, respectively. The simulation domain is a $256\times 1\times 1$ grid and periodic boundaries are imposed for all 3 directions. On average there are 256 simulation particles (sampling points) per grid cell. The grid sizes are $\Delta x=4\times 10^{-4}\mathrm{m}=0.2208\rho_{T,i}$ and $\Delta t=120\Delta x/\mathrm{c}=0.0192/\omega_{c,i}\leavevmode\nobreak\ .$ Here the time-step $\Delta t$ is relatively small compared with the cyclone period $2\pi/\omega_{c,i}$, because it needs to satisfy the Courant condition for stability. Initially the perturbed electromagnetic fields is set to zero, and electrostatic waves are generated from noise. The total number of time-steps is 8192.

Theoretically the dispersion relation of the electrostatic IBWs in the $x$-direction is Miecnikowski et al. (2018); Fried and Conte (1968)

[TABLE]

where

[TABLE]

The spectra of the electric field in the $x$-direction is plotted in Fig. 1, which clearly shows that the simulated dispersion relation matches the theoretical result very well.

To test the energy conservation property, we performed a long-term simulation. The total number of time-steps is $1\times 10^{6}$. The simulation domain is a $32\times 32\times 32$ grid mesh, and the averaged number of simulation particles per cell is 16. During the simulation the total energy is recorded, and the result is shown in Fig. 2. It is evident that the error of total energy is bounded by a small number for all simulation time-steps.

4.2 Ion temperature gradient instability and turbulence in a slab geometry

In certain parameter regimes, the ion temperature gradient in a magnetized plasma can excite the drift wave instability, which often nonlinearly evolves into a turbulent stage to produce anomalous transport of energy and particles Krall and Rosenbluth (1965); Coppi (1967); Tang (1978); Romanelli (1989); Hammett and Perkins (1990); Cowley et al. (1991); Horton (1999); Dorland and Hammett (1993); Parker et al. (1999); Rogers et al. (2000); Dimits et al. (2007); Ku et al. (2009); Merz and Jenko (2010); Sturdevant et al. (2017); Miecnikowski et al. (2018); Hu et al. (2018). We demonstrate the simulations of the ITG instability and turbulence by the structure-preserving geometric PIC algorithm in a 2D slab geometry. System parameters are similar to those in Sec. 4.1, except that the temperatures for both ions and electrons are now functions of the $x$-coordinate,

[TABLE]

The simulation domain is a $N_{x}\times N_{y}\times N_{z}=512\times 512\times 1$ grid, and the total number of time-steps is $1.2\times 10^{6}$. To balance the pressure gradient for equilibrium, we use an external electric field, which is set to

[TABLE]

It can be checked that the local Maxwell distribution function

[TABLE]

is the steady state solution of the [math]th, 1st, and 2nd-order moment equations of the Vlasov equation, i.e.,

[TABLE]

To obtain a more precise kinetic equilibrium, we first use the $f_{0}$ specified by Eq. (78) to perform a 1-D simulation, i.e., $N_{y}=1$. After $10^{6}$ time-steps when the ion distribution function reaches a steady state, we take this numerically calculated distribution function as the equilibrium distribution function for the 2D simulation in the slab geometry. For the system parameters selected in this example, the ion temperature gradient excites unstable drift modes. According to the theory of drift wave, the phase velocity in the $y$-direction of modes is approximately

[TABLE]

Plotted in Fig. 3 are the phase velocity in the $y$-direction calculated from the electric field perturbations observed in the simulation and the theoretical drift velocity $v_{d,y}\left(x\right)$ as a function of $x$ given by Eq. (79). It is clear that the simulation agrees with the theoretical predication very well. We also plotted in Fig. 3 the averaged bulk velocity of the ions in the $y$-direction, which by comparison is smaller. This indicates that the space-time structure observed in the simulation is produced by the drift wave, instead of the bulk flow of the ions.

To illustrate the instability, the time history of the amplitude of density perturbation with different $k_{y}$ at $x=224\Delta x$ are plotted in Fig. 4. We observe that all modes displayed grow initially, and saturate after $t>3000/\omega_{c,i}$. From the simulation data, the dispersion relation of the instability at $x=224\Delta x$ can be calculated. It is plotted in Fig. 5.

After a sufficient long time, the unstable modes nonlinearly evolve into a turbulent state, as evident from the distribution of kinetic energy density and number density of ions at different times. Figure 6 shows that the kinetic energy diffuses as the instability grows, saturates, and becomes turbulent. Figure 7 shows that density blobs generated by the rupture of unstable modes are the prominent structures of the fully developed ITG turbulence. The details of the instability and turbulence, especially the formation of density blobs, can be observed from the video of the density evolution available at http://staff.ustc.edu.cn/~xiaojy/ditg.html.

When turbulence develops, the energy or particle diffusion of the plasma across the magnetic field is conjectured empirically to follow the scaling of Bohm diffusion or the gyro-Bohm diffusion Taylor (1961); Drummond and Rosenbluth (1962); Waltz et al. (1990); Petty et al. (1995). The corresponding diffusion coefficients are

[TABLE]

where $\rho^{*}=\rho_{T,i}/L$ is the gyro-radius of ions measured in $L$, the characteristic length of the plasma. Assuming that the diffusion coefficient varies slowly with $x$ and the energy density $E_{k}$ diffuses according to

[TABLE]

we can calculate the numerical diffusion coefficient of ions $\chi$ in the $x$-direction. The results are plotted in Fig. 8, where the local plasma characteristic length $L$ is estimated using

[TABLE]

We observe that at $t=3841/\omega_{c,i}$ the diffusion coefficient $\chi$ is between the Bohm scaling $\chi_{B}$ and gyro-Bohm scaling $\chi_{gB}$. Afterwards, $\chi$ decreases. When the ITG turbulence is fully developed at $t=19159/\omega_{c,i}$, the diffusion coefficient $\chi$ is closer to the gyro-Bohm scaling $\chi_{gB}$.

We emphasize again that simulating long-term dynamical and transport behavior of magnetized plasmas using fully kinetic particles demands a large number of simulation time-steps, $1.2\times 10^{6}$ in this case. A structure-preserving geometric algorithm with long-term accuracy and fidelity is desirable for this purpose. To verify the long-term conservative properties of the simulation, we plotted the energy error in Fig. 9. Clearly, the error is globally bounded by a small number for all simulation time-steps.

5 Discussion and Conclusion

In conclusion, we have customarily designed a field theory for low frequency electrostatic perturbations with fully kinetic ions and adiabatic electrons, and geometrically discretized the field theory to build a structure-preserving geometric PIC algorithm for simulating drift wave instabilities and turbulence in magnetized plasmas. The geometric discretization of the field theory is accomplished using DEC, high-order Whitney interpolation forms, and the non-canonical Hamiltonian splitting method. It preserves the non-canonical symplectic structure of the particle-field system, as well as the gauge symmetry. And as a result, the PIC algorithm is automatically charge-conserving and possesses long-term conservation properties that are indispensable for simulating the dynamics of fully kinetic particles. We have successfully implemented the algorithm in the SymPIC code. The algorithm was tested and benchmarked using the examples of ion Bernstein waves and drift waves, and applied to study the ion temperature gradient instability and turbulence in a 2D slab geometry. Simulation results show that at the early stage of the ITG turbulence, the energy diffusion is between the Bohm scaling and gyro-Bohm scaling. At later time, the observed diffusion is closer to the gyro-Bohm scaling, and density blobs generated by the rupture of unstable modes are the prominent structures of the fully developed ITG turbulence.

Acknowledgements.

This research was supported by the National Key Research and Development Program (2016YFA0400600, 2016YFA0400601 and 2016YFA0400602), the National Natural Science Foundation of China (NSFC-11775219 and NSFC-11575186), China Postdoctoral Science Foundation (2017LH002), the Fundamental Research Funds for the Central Universities (WK2030040096) and the U.S. Department of Energy (DE-AC02-09CH11466).