Solving Poisson equation with slowing-down equilibrium distribution for global gyrokinetic simulation

1.   Introduction

In gyrokinetic particle simulation, the difference between particle distribution and gyro-center distribution leads to the double gyro-average of potential field, which is manifested in polarization density [1, 2]. The polarization density, which is essential for the gyrokinetic Poisson equation, depends critically on the phase structure of the equilibrium particle distribution, especially for short wavelength modes. Global gyrokinetic simulations are crucial for studying many important physics issues in magnetic fusion plasmas, such as turbulent transport scaling and turbulence spreading [3–7]. With the advent of burning plasmas, alpha particles inevitably excite Alfvénic turbulence with a slowing-down distribution via electron–alpha collisions [8–14]. A more accurate global gyrokinetic Poisson solver with slowing-down background distribution is desirable for simulating alpha particle physics in burning plasmas and neutral beam injection heating scenarios.

The multi-point average method [2, 15] is a key invention that enables global simulation of microturbulence, which has been widely used in many global gyrokinetic particle-in-cell codes [16, 17]. However, these global gyrokinetic particle simulations are based on Maxwellian background distribution. Some researchers have recently used so-called equivalent Maxwellian distribution, whose temperature or second order velocity moment stays the same as the slowing-down distribution [18], to simulate the turbulent transport of alpha particles, which may be valid for a number of physics scenarios that depend weakly on the phase space structure of the equilibrium distribution. Some other works [19], though correctly considering velocity space derivatives for alpha particles, have used the local flux-tube approximation to avoid tackling the global spatial dependence of the gyrokinetic Poisson equation.

In this work, a novel global method based on the multi-point average [2, 15] to solve the gyrokinetic Poisson equation is developed to include slowing-down equilibrium distribution for alpha particles, and the accuracy of this new method is verified in the long and short wavelength limits. This method is useful for global gyrokinetic simulation to accurately investigate alpha particle physics in burning plasmas with the advent of International Thermonuclear Experimental Reactor (ITER) operation when the alpha particles accumulate to a non-negligible portion of the ion species.

The rest of this article is structured as follows: section 2 describes how to derive polarization density in gyrokinetics with slowing-down particle distribution; section 3 shows our numerical scheme to solve the gyrokinetic Poisson equation with the slowing-down equilibrium distribution; numerical verification is provided in section 4 for the accuracy of this new scheme; and the numerical results are summarized in section 5.

2.   Gyrokinetic simulation with α particles

The gyrokinetic equation describes the time evolution of the gyrocenter distribution function $\overline{f}\left(\mathit{X},v_{||},\mu ,t\right)$ in 5D phase space [1, 20, 21], where $\mathit{X}$ is the spatial gyrocenter coordinates, $v_{||}$ is parallel velocity and $\mu$ is the magnetic moment with $\mu =m{v}_{\perp }^2/2B.$ Suppose that the gyrocenter distribution function $\overline{f}$ can be decomposed into an equilibrium component ${\overline{f}}_0$ and a perturbed component $\delta \overline{f},$ i.e. $\overline{f}={\overline{f}}_0+\delta \overline{f},$ where the overbar denotes a function of the gyrocenter coordinates. Following the usual gyrokinetic ordering, the perturbed distribution function $\delta \overline{f}$ is determined by the perturbed gyrokinetic equation [22–24]:

where the subscript $s=\mathrm{i},\mathrm{e},\alpha$ stands for particle species, i.e. ion, electron and alpha particle respectively; ${\mathit{B}}_0$ is equilibrium magnetic field with magnitude $B_0$ and direction $\mathit{b}={\mathbf{B}}_0/B_0.$ ${\mathit{B}}^{★}={\mathit{B}}_0+m_s{v}_{||}\mathrm{\nabla }\times \mathit{b}/q_s+\delta \mathit{B};$ $q_s$ is particle charge; and $m_s$ is particle mass.

The equations of motion for the gyrocenter are given by

where is magnetic drift for the guiding centers with representing magnetic field curvature, $\boldsymbol{v}_{\mathrm{E}}=-\nabla\left\langle\delta \phi_{\mathrm{gc}}\right\rangle / B_0$ is the gyro-averaged drift, and $\langle\ldots\rangle$ represents gyrophase average. The gyrocenter perturbed electromagnetic fields are given by pushing forward transform from particle coordinates to guiding center coordinates

where the gyroradius $\mathit{\rho }$ is defined as $\mathit{\rho }=\mathit{b}\times {\mathit{v}}_{\perp }/{\mathrm{\Omega }}_s$ with the gyrofrequency ${\mathrm{\Omega }}_s=q_{s}B/m$ and $\zeta$ is the gyrophase. The perturbed vector potential can be calculated by the parallel Ampère's law:

where the perturbed parallel fluid velocity $\delta u_{||s}$ can be calculated by $\delta u_{||s}={\int }^{}{\mathrm{d}}^{3}v{v}_{||}\exp \left(-\mathit{\rho }·\mathrm{\nabla }\right)\delta {\overline{f}}_s.$ The perturbed electrostatic potential $\delta \varphi$ is calculated by the quasi-neutrality condition: $q_{\mathrm{i}}\delta n_{\mathrm{i}}+q_{\alpha }\delta n_{\alpha }=e\delta n_{\mathrm{e}}.$ Due to the finite Larmor radius, the perturbed response of ions or alpha particles can be divided into two parts: the contribution from perturbed gyrocenter distribution and the contribution from particle polarization due to fluctuating electric fields:

where the guiding center density $\delta n_{s,\mathrm{gc}}$ is defined as $\mathrm{\delta }{n}_{\mathrm{i},\mathrm{gc}}={\int }^{}{\mathrm{d}}^{3}v\exp \left(-\mathit{\rho }·\mathrm{\nabla }\right)\delta {\overline{f}}_s.$ In many cases, the adiabatic response can be conveniently assumed for electrons due to their fast parallel motion, i.e. $\delta n_{\mathrm{e}}=e{n}_0\delta \varphi /T_{\mathrm{e}}$ with $T_{\mathrm{e}}$ the electron temperature. However, the perturbed density $\delta n_{\mathrm{e}}$ actually contains contributions not only from adiabatic electrons but also from non-adiabatic electrons, i.e. $\delta n_{\mathrm{e}}=\frac{e{n}_0\delta \varphi }{T_{\mathrm{e}}}+\delta n_{\mathrm{e}}^{\mathrm{NA}}.$ We also note that simply employing equations (1), (6), (7) does not provide a stable numerical algorithm for nonlinear electromagnetic simulation. In [23], a fluid–kinetic hybrid electron simulation model is invented to resolve this difficulty.

In the preceding equation, the polarization density for thermal ions is given by [1]

where $T_{\mathrm{i}}$ is the ion temperature and the double gyro-average of electric potential is defined as

where Maxwellian background distribution $f_{\mathrm{M}}$ is assumed for the thermal ions, and $\widetilde{\delta \varphi }$ has a complicated form in real space but a neat form in Fourier space:

with the function ${\mathrm{\Gamma }}_0=I_0\left(b\right){\mathrm{e}}^{-b},$ $b=k_{\perp }^2{\rho }_{\mathrm{i}}^2,$ where ${\rho }_{\mathrm{i}}=\sqrt{T_{\mathrm{i}}/m_{\mathrm{i}}}$ and $I_0$ is a modified Bessel function of the first kind. In the gyrokinetic particle simulation, $\widetilde{\delta \varphi }$ can be calculated by a multi-point average method in the real space, which enables global gyrokinetic simulation [15, 25].

However, the calculation of polarization density for alpha particles with the slowing-down background distribution remains to be unsolved for gyrokinetic particles simulation. The slowing-down distribution $f_{\mathrm{sld}}$ is the steady-state solution to the collisional scattering for an isotropic particle source with a large birth speed $v_0,$ i.e. the alpha particle source produced by the thermal nuclear fusion, where the birth energy of alpha particle $\frac{1}{2}{m}_{\alpha }{v}_0^2=3.5\mathrm{MeV}$ for typical D–T fusion. It is discovered that $f_{\mathrm{sld}}$ can be explicitly expressed as [26, 27]:

where $v_{\mathrm{c}}\equiv {\left(\frac{3\sqrt{\pi }{m}_{\mathrm{e}}}{4n_{\mathrm{e}}}{\displaystyle {\sum }_i\frac{q_{\alpha }^2{n}_{0\alpha }}{e^2{m}_{\alpha }}}\right)}^{1/3}{v}_{\mathrm{th},\mathrm{e}}$ is the critical speed at which the electron drag is comparable to the thermal ion drag, and $I_1=\frac{1}{3}\ln \left(1+v_0^3/v_{\mathrm{c}}^3\right)$ is an auxiliary function for the purpose of normalization. The alpha temperature $T_{\alpha }$ can be defined as the second velocity moment of the slowing-down distribution function, which is similar to that of the Maxwellian:

with $I_2=\frac{v_0^2}{2v_{\mathrm{c}}^2}-\frac{1}{6}\left(\frac{\pi }{\sqrt{3}}-2\sqrt{3}\arctan \frac{1-2v_0/v_{\mathrm{c}}}{\sqrt{3}}-\ln \frac{{\left(1+v_0/v_{\mathrm{c}}\right)}^3}{1-v_0^3/v_{\mathrm{c}}^3}\right).$ We note that for the fusion-born alpha particles in a 10 keV 50%–50% D-T plasma, $v_{\mathrm{c}}/v_0\approx 0.3,$ $T_{\mathrm{\alpha }}\approx 1.28m_s{v}_{\mathrm{c}}^2.$ The integrations $I_1,I_2$ have also been defined in [18].

Next, we show how to derive the form of polarization density for the alpha particles, i.e. $\delta n_{\mathrm{pol},\alpha }$ in equation (7), based on a gyrokinetic approach. Employing the Lie transform method [28], we can calculate the alpha density $n_{\alpha }$ by integrating over the velocity space on the distribution function, generated by pulling back the gyrocenter alpha distribution function $\overline{f_{\alpha }}=f_{\mathrm{sld}}+\delta {\overline{f}}_{\alpha }$,

where $\delta \tilde{\phi}_{\mathrm{gc}} \equiv \delta \phi_{\mathrm{gc}}-\left\langle\delta \phi_{\mathrm{gc}}\right\rangle$ is the gyrophase-dependent part of As is discussed, the perturbed alpha density can be separated into a perturbed gyrocenter center density and a polarization density, i.e. and the polarization density is associated with the fluctuating electric field through

We note that in this expression, the nonlinear term in $\delta n_{\alpha ,\mathrm{pol}}$ is ignored for simplicity.

To solve the gyrokinetic Poisson equation, a proper numerical algorithm is needed to deal with the guiding center transformation $\exp \left(-\mathit{\rho }·\mathrm{\nabla }\right)$ and gyrophase average in equation (11). In the local gyrokinetic simulation, this operation is carried out in Fourier space, and the associated phase angle and velocity space integration can be done analytically, which greatly simplifies numerical calculation. For global gyrokinetic simulation, the four-point or multi-point gyro-average method has been invented to calculate the ion polarization density in real space when the equilibrium distribution $f_{0\mathrm{i}}$ is Maxwellian [15]. Here we modify the original four-point average method to accommodate a slowing-down equilibrium distribution in order to investigate the self-consistent turbulent transport physics involving alpha particles. Considering the Fourier representation of the perturbed potential $\delta \varphi ={\displaystyle {\sum }_k\delta {\varphi }_k}\exp \left(\mathrm{i}\mathit{k}·\mathit{x}\right)$ in equation (11), then we can obtain:

where and $J_0= J_0\left(k_{\perp} \rho_{\mathrm{i}}\right)=\langle\exp (\mathrm{i} \boldsymbol{k} \cdot \boldsymbol{\rho})\rangle$ is the zeroth-order Bessel function, and is defined as

which can be considered as the expectation of $\lambda$ weighted by the equilibrium distribution $f_{\mathrm{sld}}$ after double gyro-averaging due to the back and forth transformation between particle position and gyrocenter position, where ${\tilde{f}}_{\mathrm{sld}}\left(v_{\perp }\right)={\int }^{}\mathrm{d}{v}_{||}\lambda f_{\mathrm{sld}}.$ Unlike the Maxwellian case, ${\mathrm{\Gamma }}_{\alpha }$ does not have a simple analytic expression in Fourier space and has to be evaluated numerically. In principle, one can solve the quasi-neutrality equation, i.e. equation (7), using equations (12) and (13) in Fourier space. Figure 1 shows the polarization term $c_{\alpha }-{\mathrm{\Gamma }}_{\alpha }\left(k_{\perp }{\rho }_{\mathrm{c}}\right)$ in Fourier space between Maxwellian (red solid line) and slowing-down (blue dashed line) with identical temperature but different $v_0/v_{\mathrm{c}}.$ From bottom to top, we have $v_0/v_{\mathrm{c}}$ ranging from 2 to 5 with a 0.5 step size. The polarization calculated from slowing-down and that from the equivalent Maxwellian are the same for long wavelengths but differ while $k_{\perp }{\rho }_{\mathrm{c}}\gtrsim 1.$

Figure  1.  Polarization comparison between slowing-down distribution and equivalent Maxwellian in k-space. The blue dashed line is for the slowing-down distribution and the red solid line is for the equivalent Maxwellian distribution. The lines from bottom to top correspond to various ratios of $v_0/v_{\mathrm{c}}$ ranging from 2 to 5 with a step size of 0.5.

DownLoad: Full-Size Img PowerPoint

With these newly defined functions, the gyrokinetic Poisson equation can be written in a dimensionless form

where $\widetilde{\delta {\varphi }_{\alpha }}$ is defined by the following Fourier representation:

For the typical D–T plasmas parameters mentioned before in the definition in temperature of slowing-down distribution, the polarization density of alpha particles is enhanced by the mass ratio $m_{\alpha }/m_{\mathrm{i}}$ for long wavelengths if the equilibrium density ratio is given. However, the ratio $n_{\alpha ,\mathrm{pol}}/n_{\mathrm{i},\mathrm{pol}}$ goes down as $k_{\perp }$ become larger, as shown in figure 2.

Figure  2.  Polarization density ratio between alpha particles and thermal fuel ions varying in $k$ space, normalized by their equilibrium densities respectively.

DownLoad: Full-Size Img PowerPoint

This Fourier representation is neat in form but it mixes up the configuration space and velocity space dependences through the $J_0$ term. In reality, the background magnetic field and perpendicular temperature can vary in real space, and then ${\mathrm{\Gamma }}_{\alpha }$ will gain global spatial dependence. Besides, the Fourier transform approach makes it more difficult to deal with realistic tokamak geometry, where no periodicity exists in the radial direction and on many occasions the global effects have to be considered seriously [3]. For the Maxwellian background distribution, the multi-point/four-point gyro-average method has been developed to solve this gyrokinetic Poisson equation in real space [1, 2]. Here we modify this method by including the slowing-down background distribution $f_{\mathrm{sld}}$ as the equilibrium distribution in the gyrokinetic Poisson equation.

3.   Multi-point method for gyrokinetic Poisson solver with slowing-down distribution

The crucial part of implementing this multi-point average method for a gyrokinetic Poisson solver is to represent $\widetilde{\delta {\varphi }_{\alpha }}$ in equation (15) by the values of $\delta \varphi$ at various field points in real space. By numerical interpolation, we note that $\widetilde{\delta \varphi }$ can be expressed as a linear combination of the $\delta \varphi$ values on a number of nearby grid points and consequently equation (14) is transformed into a discrete matrix form such as $\mathit{A}·\mathit{x}=\mathit{b},$ which can then be solved by many known matrix inversion algorithms, such as the Krylov subspace solver, provided by an existing parallel library like PETSc [29].

Starting from the integral form of $\widetilde{\delta {\varphi }_{\alpha }}$ instead of the Fourier form in equation (15), one finds that

To calculate at a grid point for a specific one needs to evaluate the gyro-averaged function $\left\langle\exp (-\rho \cdot \nabla)\left\langle\delta \phi_{\mathrm{gc}}\right\rangle\right\rangle,$, which is the average value of $\left\langle\delta \phi_{\mathrm{gc}}\right\rangle$ on a ring with radius around as shown by the dotted circle in figure 3. The gyro-averaged quantity $\left\langle\delta \phi_{\mathrm{gc}}\right\rangle$ can also be calculated by this ring average method, e.g., the value of $\left\langle\delta \phi_{\mathrm{gc}}\right\rangle$ at the black triangle in figure 3 can be calculated by the average value on a solid circle. It is not necessary to actually integrate numerically along the whole ring to compute the gyrophase average, which would make the gyro-average process rather time-consuming and expensive. According to [1, 2], a selection of four points uniformly distributed on the ring (four-point average method) is sufficient to compute the gyro-average for wavelengths up to Thus, nine neighboring points are required to compute on the grid point, as shown by the eight red points and the central blue diamond in figure 3. This scheme can be modified to adopt 4N (N = 1, 2 ···) points on the ring straightforwardly. In more general geometry, these points required for the gyro-average computation may not lie exactly on the grids, but their values can be acquired by a linear interpolation of the nearby grid points. Finally, a few rings are summed with different values with the weight function and the relationship between and on each grid point is found.

Figure  3.  Scheme for calculating $\widetilde{\delta \varphi }$ at grid point $i,j.$

DownLoad: Full-Size Img PowerPoint

The remaining issue in evaluating equation (16) is how to discretize the $v_{\perp }$ integral with the weight function ${\tilde{f}}_{\mathrm{sld}}\left(v_{\perp }\right).$ Here we approximate the integral by a weighted summation by choosing a few sampling grid points along the $v_{\perp }$ coordinate. From the definition of $\widetilde{\delta \varphi },$ one can tell that it is equivalent to the approximate equation (13) by:

where $c_j$ are the summing weights due to $\tilde{f}\left(v_{\perp }\right)$ and $v_{\perp j}$ are the sampling grid points. The value pairs of $\left(c_j,v_{\perp j}\right)$ are chosen by minimizing the following error function:

Here $a$ is the maximum value of $k_{\perp }{\rho }_{\mathrm{c}}$ that we are interested in. Since low-frequency micro-turbulence usually peaks around $k_{\perp }{\rho }_{\mathrm{c}}<1,$ it is required that this approximation has a better accuracy for long wavelengths or $k_{\perp }{\rho }_{\mathrm{c}}\ll 1.$ Considering the Taylor expansion for $J_0\left(x\right)$ and ${\mathrm{\Gamma }}_{\alpha }\left(x\right)$ around $x~0,$ one finds that $J_0\left(x\right)=1-x^2/4+O\left(x^4\right)$ and ${\mathrm{\Gamma }}_{\alpha }\left(x\right)=c_0-\frac{T_{\alpha }}{m{v}_{\mathrm{c}}^2}{x}^2+O\left(x^4\right).$ Let the first two terms equal each other:

These two constraints are then used to reduce the degrees of freedom. In order to minimize $\mathit{ϵ}$ with respect to $\left(c_j,v_{\perp j}\right),$ we use the Nelder–Mead method [30], which is a gradient-free iterative optimization algorithm. $I_{1,2,3}$ are functions of $v_{\mathrm{c}}/v_0,$ which is chosen to be 0.3 here to show the numerical result of $\left(c_j,v_{\perp j}\right).$ In the one-velocity-node case, we find that $c=1.226$ with the velocity node $v_{\perp }/v_{\mathrm{c}}=1.443$ and the relative error is $3.6\%$ for $k_{\perp }{\rho }_{\mathrm{c}}<0.5.$ When using two velocity nodes, we find that $c_1=0.9347$ and $c_2=0.2910$ with the velocity nodes $v_{\perp 1}/v_{\mathrm{c}}=0.8778$ and $v_{\perp 2}/v_{\mathrm{c}}=2.510,$ and the relative error is about $3.6\%$ for $k_{\perp }{\rho }_{\mathrm{c}}<1.5.$ In the three-velocity-node case, we find that $\left(c_1,c_2,c_3\right)=\left(0.1186,0.3881,0.7190\right)$ with $\left(v_{\perp 1}/v_{\mathrm{c}},v_{\perp 2}/v_{\mathrm{c}},v_{\perp 3}/v_{\mathrm{c}}\right)=\left(0.7016,1.716,2.984\right),$ and the relative error is only $0.46\%$ for $k_{\perp }{\rho }_{\mathrm{c}}<2.$ The three-velocity-node approximation is compared with the exact value from direct numerical integration, as shown in figure 4. Satisfactory accuracy is achieved with a relative error less than $0.46\%$ for $k_{\perp }{\rho }_{\mathrm{c}}<2,$ which is sufficient to include the most interesting finite Larmor radius effects due to the slowing-down alpha particles. We find that using more velocity nodes can achieve better accuracy for larger $k_{\perp }{\rho }_{\mathrm{i}},$ e.g., using four velocity nodes can make the ${\mathrm{\Gamma }}_{\alpha }$ function from the four-point-average method sufficiently accurate up to $k_{\perp }{\rho }_{\mathrm{i}}~4.$ Notice that $v_{\mathrm{c}}$ can be varied at equilibrium scale and the weights $c_{\mathrm{i}}$ and velocity nodes $v_{\perp \mathrm{i}}$ will change as $v_{\mathrm{c}}/v_0$ changes. We can compute a series of weights and velocity nodes of different $v_{\mathrm{c}}/v_0,$ then generate a spline with these discrete values. When constructing the matrix corresponding to ${\mathrm{\Gamma }}_{\alpha },$ the local values of weights and velocity nodes can be obtained from the pre-calculated spline. In this way, the four-point average method can capture the non-uniformity of the equilibrium quantities. As a contrast, this non-uniformity affects only $\rho$ but not the coefficients $\left(c_j,v_{\perp j}\right)$ for the Maxwellian case [15].

Figure  4.  Exact ${\mathrm{\Gamma }}_{\alpha }$ function (blue solid line) and its numerical approximations vary with perpendicular wavelength $k_{\perp }{\rho }_{\mathrm{c}}:$ 4-point average method with three velocity nodes in the integration (red dotted line), and Padé approximation (black dashed line).

DownLoad: Full-Size Img PowerPoint

We also test the widely used Padé approximation for the thermal ions, and find that it can introduce a $10\%$ relative error near $k_{\perp }{\rho }_{\mathrm{c}}~1.5$ comparing to the exact solution. Figure 5 shows the comparison between the four-point average method and the Padé approximation with the following form:

Figure  5.  Comparison of analytic expression and numeric solutions using four-point average approximation and Padé approximation for the gyrokinetic Poisson equation in the long wavelength limit with $k_r{\rho }_{\mathrm{i}}=0.11$ and $m=6.$

DownLoad: Full-Size Img PowerPoint

We note that the Padé approximation with the Maxwellian distribution has a similar form with $c_{\alpha }=1.$ In the long wavelength limit, these approximations are both very close to the exact value, as shown in figure 5, and it can be further verified by the numerical benchmarks shown in the next section.

Although the method demonstrated here is applied to the isotropic background distribution, it can also be generalized for more complex anisotropic background distributions. To obtain the ${\mathrm{\Gamma }}_{\alpha }$ function in equation (13) for the polarization density term, we just need to perform $\partial /\partial \mu$ and integrate it over velocity space with the anisotropic distribution function, which can be performed even for the numerical distribution functions. This will be carried out in our future work.

4.   Numerical verification

To verify our global algorithm for the slowing-down background distribution, we shall solve the gyrokinetic Poisson equation without electron response in a large-aspect-ratio toroidal geometry tokamak with circular cross section as a sample problem. We have implemented this multi-point method for the gyrokinetic Poisson solver in Gyrokinetic Toroidal Code (GTC) by loading slowing-down marker distribution, modifying weighting coefficients and using its embedded PETSc Library for the gyrokinetic Poisson solver. The simulation domain is chosen to be $r\in \left[a_0,a_1\right]$ on the poloidal plane, with $a_0=0.142a,$ $a_1=0.938a$ and ${\rho }_{\alpha }/a\approx 2.51\times {10}^{-2}.$ In the GTC, the accuracy of gyrokinetic Poisson solver has been verified for thermal ions and electrons [31]. Thus, we can verify the accuracy of gyrokinetic Poisson solver for alpha particles by ignoring contributions from the thermal ion and electron in equation (14). In the long wavelength limit equation (14) can be reduced to the following form using the Taylor expansion of ${\mathrm{\Gamma }}_{\alpha }:$

In the large safety factor limit, the toroidal effect can be ignored and ${\mathrm{\nabla }}_{\perp }^2=\frac{1}{r}\frac{\partial }{\partial r}r\frac{\partial }{\partial r}+\frac{1}{r^2}\frac{{\partial }^2}{\partial {\theta }^2}$ in polar coordinates on the poloidal plane since the difference between the perpendicular plane and poloidal plane can be ignored in the large-aspect-ratio limit. After these simplifications, equation (22) is just a normal Poisson equation and we can choose $\delta n_{\alpha ,\mathrm{gc}}$ to be the eigenfunction of the Laplacian operator to ensure an analytic solution. Let $\delta n_{\mathrm{i},\mathrm{gc}}/n_{\mathrm{i}}=\left(J_m\left(k_{0}r\right)-Y_m\left(k_{0}r\right)J_m\left(k_0{a}_1\right)/Y_m\left(k_0{a}_1\right)\right)\cos m\theta$ in which $k_0$ satisfies $J_m\left(k_0{a}_0\right)Y_m\left(k_0{a}_1\right)-Y_m\left(k_0{a}_0\right)J_m\left(k_0{a}_1\right)=0.$ Here $J_m$ and $Y_m$ are mth-order Bessel functions of the first and second kinds, respectively. Then the solution of this Poisson equation is just $\delta \varphi =k_0^{-2}\left(J_m\left(k_{0}r\right)-Y_m\left(k_{0}r\right)J_m\left(k_0{a}_1\right)/Y_m\left(k_0{a}_1\right)\right)\cos m\theta$ with zero boundary condition on $r=a_0,a_1.$ With $m=6$ as an example, the comparison between the analytic solution and numeric solution along the line $\theta =0$ is shown in figure 5. The relative differences are very small at 0.6% between the analytic and four-point average and 1.3% between the analytic and Padé approximation. More generally, the poloidal cross section contour for the solution is shown in figure 6. The difference is negligibly small between the numerical solution from the four-point average and that from the Padé approximation, as shown in figure 6(c). Thus, in the long wavelength limit, our four-point average method works perfectly for the slowing-down equilibrium distribution.

Figure  6.  2D poloidal contours for solutions to different operators to the gyrokinetic Poisson equation in the long wavelength limit: (a) four-point average operator, (b) Padé approximation operator, where $k_r{\rho }_{\mathrm{i}}=0.11$ and $m=6.$ The differences between these two solutions are shown in (c).

DownLoad: Full-Size Img PowerPoint

In order to simulate short wavelength modes, we need to verify the validity of our algorithm in the short wavelength limit. We solve equation (14) with only slowing-down alpha particles in a natural unit $T_{\alpha }=T_{\mathrm{e}}=B_0=q_{\alpha }=e=1,$ with $T_{\alpha },T_{\mathrm{e}}$ and $B_0$ being the quantity at the magnetic axis. Using eigenfunctions of the Laplacian operator, $\delta \varphi \left(r,\theta \right)$ can be expressed as a series:

Here, is the ith positive value satisfying and the boundary condition is fulfilled. The gyro-averaged potential $\left\langle\delta \phi_{\mathrm{gc}}\right\rangle$ can now be expanded as

Then, following the same procedure of how equation (12) is derived, the double gyro-averaged term is

In order to compare this analytic solution with that from the four-point average method, we keep the magnetic field and temperature profile constant and thus ${\rho }_{\mathrm{c}}$ is spatially uniform. If the density is given by a summation of the aforementioned eigenfunctions, then equation (14) can be solved analytically. Here we sum up three eigenfunctions with $m=6,$ and different radial wavelengths at $k_r{\rho }_{\mathrm{c}}=1,1.5,2,$ and corresponding amplitudes to $1,0.5$ and $0.2$ respectively. Figures 7 and 8 show the 1D and 2D plots for the numerical solutions from the four-point average method and Padé approximation, respectively. The globally averaged differences between the numerical methods and the analytic solution are about 1% and can be ignored. Figure 7 also shows the solution from the original four-point average method with an equivalent Maxwellian with the same temperature as the slowing-down distribution, i.e. the orange dotted-and-dashed line. This numerical solution deviates from the analytic solution with an overall error of more than 6%. The amplitude of the solution using the four-point average method is slightly larger, which can be ascribed to the fact that the operator of the four-point average is larger than the Padé approximation in $k$ space, as shown in figure 4. Though not so precise in $k$ space, the Padé approximation still gives a solution as good as that from the four-point average method in this case. This suggests that we may not need to be very accurate at large $k_{\perp }{\rho }_{\alpha }$ since the polarization density is nearly identical to the adiabatic response in the short wavelength limit, regardless of the particle's velocity distribution.

Figure  7.  Comparison of analytic solution (blue solid line) and numeric solutions using the new four-point average method with slowing-down equilibrium (red X marker), original four-point average method with equivalent Maxwellian equilibrium (orange dotted-and-dashed line) and Padé approximation (black long dashed line) for the gyrokinetic Poisson equation in the short wavelength limit.

DownLoad: Full-Size Img PowerPoint

Figure  8.  2D contour of the solution to gyrokinetic Poisson equation in the short wavelength limit on the poloidal plane. The numerical methods used in solving the Poisson equation are the (a) four-point average method and (b) Padé approximation. The difference between them is shown in (c).

DownLoad: Full-Size Img PowerPoint

Finally, we investigate how the fraction of energetic particle affects the gyrokinetic Poisson solver. Here we assume that the main ions are a 50%–50% mixture of D and T with a temperature of 10 keV, and the energetic particles are the fusion-born alpha particles with a slowing-down distribution. The density term is the same as that in the above short wavelength verification where $k_r{\rho }_{\mathrm{c}}~1$ and in this case $k_r{\rho }_{\mathrm{i}}\ll 1.$ The solutions are shown in figure 9 with different energetic particle fractions from 5% to 20%. The difference could be significant when $n_{\alpha }/n_{\mathrm{e}}$ increases above 5%, as demonstrated in the figure.

Figure  9.  Solutions to an arbitrary density source with a dominant wavelength $k_r{\rho }_{\mathrm{c}}~1$ with different energetic alpha particle fraction from 5% to 20%.

DownLoad: Full-Size Img PowerPoint

5.   Conclusion

A real space gyrokinetic Poisson solver for slowing-down equilibrium distribution has been developed based on the multi-point average method [2, 15] and verified for its accuracy in the long and short wavelength limits. The discovery process for this method is shown in detail and it can be further modified to accommodate more complicated equilibrium distributions. This method can be incorporated into the global gyrokinetic particle simulation to study the crucial alpha particle physics in the burning plasmas, i.e. to simulate the drift Alfvénic turbulence accurately in the presence of slowing-down alpha particle distribution.