Quasilinear diffusion coefficients in a finite Larmor radius expansion for ion cyclotron heated plasmas

The increase of the kinetic energy density of plasmas due to the energy transfer from RF waves can be described by Brambilla and Krucken (1988)

where $q$ is the species charge, $\widetilde{f}$ is the perturbed distribution function of the species due to the RF waves and $\left\langle...\right\rangle_{w}$ indicates the average over a number of wave periods in time and space. Here, $\omega$ is the RF wave frequency, and $\mathbf{v}^{\prime}$ and $\mathbf{r}^{\prime}$ denote the velocity and space vector at the past time $t^{\prime}$, respectively. Using the solution of the linearized Vlasov equation for $\widetilde{f}$ and the time harmonic form for the frequency $\omega$, the energy increase rate (so-called Wdot) is

where $m$ is the mass of the species, $\gamma$ is the damping rate of the wave, and $f$ is the background distribution function. This general expression of the wave damping can be reformulated differently in each numerical code depending on its specific numerical formation and assumption. Using the Fourier spectral representation of the space coordinate, it may be described by

where $\mathbf{W_{l}}$ is the resonance kernel and, for example, it is defined in Eq. (11) of Jaeger et al. (2006a) for AORSA.

In this paper, we represent the quasilinear diffusion tensor in spherical coordinates, $(v,\vartheta,\phi)$, where $v=\sqrt{v_{\perp}^{2}+v_{\|}^{2}}$ is the speed, $\vartheta=\arctan(v_{\perp}/v_{\|})$ is the pitchangle, and $\phi$ is the gyroangle. Here, $v_{\perp}$ and $v_{\|}$ are the velocity perpendicular and parallel to the static magnetic field, respectively. This coordinate is beneficial to evaluate the quasilinear diffusion tensor because it can reduce the computation time due to the adiabatic invariant, $v$, as will be explained in Section IV. Then, the quasilinear diffusion coefficients (B, C, E and F) determine the divergence of the flux in velocity space by Kerbel and McCoy (1985)

Because the gyro-averaged quasilinear diffusion coefficients (averaged in $\phi$) are sufficient to model the energy transfer, we define the four coefficients in $Q(f)$ and do not retain the flux in the gyro-phase direction. For the perpendicular momentum transfer, the flux along the gyro-phase direction should be included Lee et al. (2012), but it is not interest of this paper.

The energy transfer using the quasilinear diffusion coefficient is

Using the resonance kernel in Eqs. (3) and (5), we can evaluate the quasilinear diffusion coefficients for any numerical formulation. Note that the Kennel-Engelmann quasilinear diffusion operator is derived in uniform plasmas with a homogenous static magnetic field, as shown in Appendix A. In a toroidal geometry, we redefine $\dot{W}$ using the local space coordinate $(x,y,z)$ where $x$ and $y$ coordinates are orthogonal in the plane that is perpendicular to the static magnetic field along the $z$ coordinate. The unit vector for $(x,y,z)$ coordinates are $\mathbf{e_{x}}$, $\mathbf{e_{y}}$, and $\mathbf{e_{\|}}$, respectively. For convenience, we define the rotating coordinate $\mathbf{e_{+}}=(\mathbf{e_{x}}+i\mathbf{e_{y}})/\sqrt{2}$ and $\mathbf{e_{-}}=(\mathbf{e_{x}}-i\mathbf{e_{y}})/\sqrt{2}$ and the electric field $E_{+}=\mathbf{E}\cdot\mathbf{e_{+}}$ and $E_{-}=\mathbf{E}\cdot\mathbf{e_{-}}$ Brambilla and Krucken (1988). For example, the diagonal components of the quasilinear tensor in the speed direction $v$ is

where $\epsilon$ is the vacuum permittivity, $\Omega$, $\omega_{p}$ and $n_{s}$ are the gyrofrequency, the plasma frequency and the density of the species, respectively, and $k_{\|}=\mathbf{k}\cdot\mathbf{e_{\|}}$ is the parallel wavevector. The dirac-delta function is obtained by Plemelj relation Stix (1992) for $(\omega-n\Omega-v_{\|}k_{\|})^{-1}$ term in the resonance kernel in Eq. (3) and the contribution of the wave polarization is determined by the effective potential $\chi_{\mathbf{k},n}$ (see Appendix A),

where $J_{n}=J_{n}(k_{\perp}\rho_{i})$ is the first kind Bessel function for the order n. The relations between the quasilinear coefficients are determined by the diffusion direction. The diffusion direction is fixed by $G(f)=0$ Kennel and Engelmann (1966); Stix (1992), where the operator $G$ in Eq. (67) can be also described in spherical coordinates,

The operator $G$ results in the Onsager relations of the coefficients Kerbel and McCoy (1985):

These relations are important when we derive the expansion of the quasilinear diffusion coefficients in terms of small Larmor radius in the following sections. The relations should be preserved to fix the diffusion direction in any approximation. Because the relations do not depend on the perpendicular wavevector, $k_{\perp}$, it can be used in the FLR expansion.

When the ion Larmor radius is much smaller than the perpendicular wavelength, $\dot{W}$ can be expanded by Brambilla and Krucken (1988),

where the SCK approximation is used as in Eq. (23) of Brambilla and Krucken (1988). In Brambilla and Krucken (1988) the electromagnetic contribution is ignored for the Maxwellian distribution function $f=f_{M}(v)$ (i.e. $(\mathbf{v}\times\mathbf{B})\cdot\nabla_{v}f_{M}=0$), but it is included in Eq. (12) for an arbitrary distribution function $f(v,\vartheta)$. For the small Larmor radius compared to the perpendicular wavelength (i.e. $k_{\perp}\rho_{i}\ll 1$), we use the Taylor series for the position vector $\mathbf{r}$ from the guiding center position $\mathbf{r_{g}}$. In this approximation, we expand it by

where the number in the parenthesis of the superscript denotes the order in the small parameter $k_{\perp}\rho_{i}$. The velocity in the lowest order is

where $\phi$ is the gyroangle between $v_{\perp}$ and $v_{x}$ in the perpendicular velocity plane at the time $t$. The position vector from the gyro-center in the first order is

In Eq. (12), the electrostatic contribution is associated with both quasilinear coefficients B and C by

where $E_{\perp}(\mathbf{r}^{\prime})=\mathbf{E}\cdot(\mathbf{v}^{\prime}-v_{\|}^{\prime}\mathbf{e_{\|}})/v_{\perp}^{\prime}=\{E_{+}\exp({i(\phi+\int_{t}^{t^{\prime}}\Omega^{\prime}d\tau)})+E_{-}\exp({-i(\phi+\int_{t}^{t^{\prime}}\Omega^{\prime}d\tau)})\}/\sqrt{2}$. The electromagnetic contribution depends on only the coefficient C that is associated with the pitch-angle direction variation of the distribution function ($\partial f/{\partial\vartheta}$) by

where $\nabla_{\perp}=\{\exp({i(\phi+\int_{t}^{t^{\prime}}\Omega^{\prime}d\tau)})\partial_{+}+\exp({-i(\phi+\int_{t}^{t^{\prime}}\Omega^{\prime}d\tau)})\partial_{-}\}/{\sqrt{2}}$, $\partial_{+}=\mathbf{e_{+}}\cdot\nabla$, $\partial_{-}=\mathbf{e_{-}}\cdot\nabla$, and $\nabla_{\|}=\partial_{\|}=\mathbf{e_{\|}}\cdot\nabla$.

For simplicity of the representation and the connection to the coefficients, we can separate the contributions on $\dot{W}$ by each harmonic gyrofrequency ($n$) and by the diffusion directions associated with B and C. The B part that is associated with $\partial f/\partial v$ for the ion fundamental damping ($n=1$) is Brambilla (1989)

where the term in the first line on the right hand side of Eq. (18) is for the lowest order $O(1)$, the term in the second line is in the first order of $O(k_{\perp}\rho_{i})$, and the remaining terms are in the second order of $O((k_{\perp}\rho_{i})^{2})$.

The C part that is associated with $\partial f/\partial\vartheta$ for ion fundamental damping ($n=1$) is determined by both electrostatic and electromagnetic parts (i.e. $\dot{W}^{n=1}_{C}=\dot{W}^{n=1}_{C,ES}+\dot{W}^{n=1}_{C,EM}$). The electrostatic part is

and the electromagnetic part is

where the first line is for the lowest order $O(1)$, the second line is for the first order of $O(k_{\perp}\rho_{i})$, and the remaining is for the second order of $O((k_{\perp}\rho_{i})^{2})$ in both Eq. (19) and (20).

Before formulating the quasilinear diffusion coefficients using the $\dot{W}^{1}_{B}$ and $\dot{W}^{1}_{C}$, we note three important problems that occur in the expansion of $\dot{W}$:
(1) The coefficient B needs to be positive-definite so that the energy is transferred from the waves to plasmas for the Maxwellian equilibrium distribution and the H-theorem by K-E coefficients is guaranteed Kennel and Engelmann (1966). However, only the lowest order term in $\dot{W}^{1}_{B}$ guarantees the positive definite property because of the square form. This problem is also explained in Brambilla (1989).
(2) The zero order in $\dot{W}^{1}_{C}$ in the first line has the both contribution from $E_{+}$ and $E_{\|}$, while it is only determined by $E_{+}$ in the Kennel-Engelmann form in the lowest order Kennel and Engelmann (1966). It results in the different contribution from the wave polarization even in the lowest order. In $\chi_{\mathbf{k},n=1}$ of Eq. (7) for the K-E coefficient, only $E_{\mathbf{k},+}J_{0}$ results in the zero order contribution.
(3) The relation between B and C is different from the relation between K-E form of Eq. (9), which results in the different direction of diffusion.

In the next section, we will resolve these three problems to derive the quasilinear diffusion coefficients correctly in the small Larmor radius limit.

A solution of the first problem described in the Section II is to retain only the lowest order in the quasilinear diffusion. The lowest order for the fundamental damping ($n=1$) is the zero order in the FLR expansion $O(1)$, while the lowest order for the second harmonic damping ($n=2$) is the second order $O((k_{\perp}\rho_{i})^{2})$. Keeping the lowest order is sufficient unless the parameter $k_{\perp}\rho_{i}$ of the fast ions is too large, as will be shown in Section V. The parameter $k_{\perp}\rho_{i}$ is determined by the wave power density, the ion density and mass, and the static magnetic field strength.

Another advantage of selecting only the lowest order term is the compatibility with the dielectric tensor for the current. The current can be derived in FLR expansion Brambilla (1989) by

where $\sum_{s}$ is the summation over the species. In the finite Larmor radius approximation,

the zero Larmor radius current is

The $n=1$ fundamental damping is determined by the $\mathbf{e_{+}}$ in the lowest order. After gyro-average in $\phi$, it can be described by the non-local operators $\hat{L}$, $\Delta\hat{L}_{1}$ and $\Delta\hat{L}_{2}$,

where $\mu_{0}$ is the vacuum permeability and the integral operator $\hat{L}$ Brambilla (1999) is

The electromagnetic contribution to the integral operator are $\Delta\hat{L}_{1}$ and $\Delta\hat{L}_{2}$,

Because of the cancelation by full FLR terms as will be shown in the next subsection, the lowest order of $\dot{W}^{n=1}_{C,EM}$ has only the term with $E_{+}$ in the first line of Eq. (20), and the $\Delta\hat{L}_{2}E_{\|}$ only contributes to the higher order $O(k_{\perp}\rho_{i})$. Then, to the lowest order, one can prove that

where we used the SCK approximation, giving the electric field for the trajectory along the static magnetic field by $\mathbf{E}^{\prime}=\mathbf{E}\left(\mathbf{r}-\int^{t}_{t^{\prime}}v_{\|}\mathbf{e_{\|}}\prime d\tau\right)$ and

In the general relation, there exists the contribution of kinetic flux making the difference between $\dot{W}$ and $\langle\mathbf{E}\cdot\mathbf{J}\rangle_{w}$,

where the kinetic flux $\mathbf{T}$ has been derived in many studies Brambilla and Krucken (1988); Smithe (1989). However, to the lowest order, this kinetic flux contribution vanishes in the $n=1$ damping. Also, the $\dot{W}$ and $\langle\mathbf{E}\cdot\mathbf{J}\rangle_{w}$ in Eq. (28) are equivalent to the $\dot{W}$ in Eq. (14) of Jaeger et al. (2006a) by making $J_{0}(k_{\perp}\rho_{i})=1$ and $J_{n>0}(k_{\perp}\rho_{i})=0$ in the FLR approximation. For the iteration between TORIC-CQL3D, this equivalence between $\dot{W}$ and $\langle\mathbf{E}\cdot\mathbf{J}\rangle_{w}$ in the lowest order makes the numerical implementation much simpler, and we can use the existing implementation for non-Maxwellian dielectric tensor in the lowest order Phillips et al. (2003); Bertelli et al. (2017).

The second and third problems described in the Section II can be resolved by considering the cancelations of the terms using the summation of all FLR expansions, as described in Chapters 10 and 17 of Stix (1992). In this subsection, we explain the cancellation, which can be applicable for both forms in $\dot{W}^{n=1(0)}_{B}+\dot{W}^{n=1(0)}_{C}$ and $(\hat{L}+\Delta\hat{L_{1}})E_{+}+\Delta\hat{L}_{2}E_{\|}$. The methods of cancellation in both forms are the same, so we only derive it for the latter form.

We note two facts in the form of $(\hat{L}+\Delta\hat{L_{1}})E_{+}+\Delta\hat{L}_{2}E_{\|}$. First, the operator $(\hat{L}+\Delta\hat{L_{1}})$ has the operator $U$,

where $k_{\|}$ comes from $\partial_{\|}$ in Eq. (26). The operator $U$ guarantees the diffusion direction of K-E coefficients towards $G=0$. Secondly, the $\Delta\hat{L}_{2}E_{\|}$ term cancels with other contributions of $E_{\|}$ (on $P$ operator in Stix notation Stix (1992); Brambilla (1999)) by

where $\lambda=k_{\perp}\rho_{i}$. Because the term for the resonance condition $\omega-k_{\|}v_{\|}-n\Omega$ also exists in the denominator of the operators $\Delta\hat{L}_{2}$ and $P$, the second term in the last line of Eq. (32) vanishes for any distribution function $f$. The remaining term depends on the $U$ operator and it is in the higher order of $O(k_{\perp}\rho_{i})$ for n=1 because of $J_{1}(\lambda)\sim k_{\perp}\rho_{i}/2$. Here, when we derive Eq. (32), we utilize the Bessel function expansion for the sinusoid phase in Eq. (65) and the following Bessel function identities are used for the cancellation

Here, note that the full summation over $n$ number is required in the cancellation. The small contributions from the higher orders $n>1$ accumulate for the cancelation in the lower order.

Because the operator $U$ only remains for both $(\hat{L}+\Delta\hat{L_{1}})$ and $\Delta\hat{L}_{2}E_{\|}$, the diffusion direction in the relations of Eq. (9-11) still holds. Also, to the lowest order, the coefficient B for n=1 is only determined by $\hat{L}$, giving