On the validity of the guiding-center approximation in the presence of strong magnetic gradients

TL;DR

This paper investigates the accuracy of the guiding-center approximation for charged particle motion in nonuniform magnetic fields with gradients, demonstrating its validity even beyond traditional small-parameter limits through exact solutions.

Contribution

It provides an exact solution for particle motion in a magnetic field with a gradient and compares it to guiding-center theory, showing the approximation's robustness beyond standard assumptions.

Findings

Guiding-center theory agrees well with exact solutions in this scenario.

The approximation remains valid even when the small parameter condition is relaxed.

The work bridges the gap between exact particle motion and approximate guiding-center predictions.

Abstract

The motion of a charged particle in a nonuniform straight magnetic field with a uniform magnetic-field gradient is solved exactly in terms of elliptic functions. The connection between this problem and the guiding-center approximation is discussed. It is shown that, for this problem, the predictions of guiding-center theory agree very well with the orbit-averaged particle motion and hold well beyond the standard guiding-center limit $\epsilon \equiv \rho/L \ll 1$, where $\rho$ is the gyromotion length scale and $L$ is the magnetic-field gradient length scale.

Equations

$$\mathbf{B}(y)=B_0\left(1-\frac{}{}y/L\right)\hat{\mathsf{z}},$$

$$x^{\prime \prime }$$

$$y^{\prime \prime }$$

$$\left(x_0,x_0^{\prime };\frac{}{}{y}_0,y_0^{\prime }\right)=(0,1;1,0).$$

$$\overline{y}(ϵ)\equiv \frac{1}{4\mathsf{K}({ϵ}^2)}{\int }_0^{4\mathsf{K}({ϵ}^2)}y(\phi ,ϵ)d\phi ,$$

$$\mathsf{K}(m)\equiv {\int }_0^{\pi /2}\frac{d\phi }{\sqrt{1-m{\sin }^2\phi }}.$$

$${\mathsf{T}}_{\mathrm{gc}}\mathbf{X}\equiv \mathbf{X}+G_1^{\mathbf{X}}+G_2^{\mathbf{X}}+\frac{1}{2}{\mathsf{G}}_1\cdot \mathsf{d}{G}_1^{\mathbf{X}}+\cdots$$

$$⟨\mathbf{x}{⟩}_{\mathrm{gc}}-\mathbf{X}$$

$$\frac{1}{\rho}\,(\langle y\rangle_{\rm gc}-Y)=\widehat{\sf y}\,\mbox{\boldmath$\cdot$}\,\left(-\;\frac{\mu B_{0}}{2\,m\rho\Omega^{2}}\;\nabla\ln B\right)=\frac{\epsilon}{4}.$$

$$\widehat{{\sf b}}\,\mbox{\boldmath$\times$}\,\left(\frac{\mu\,B_{0}}{m\rho\Omega^{2}}\,\nabla\ln B\right)\;=\;\widehat{\sf z}\,\mbox{\boldmath$\times$}\,\left(-\,\frac{\epsilon}{2}\,\widehat{\sf y}\right)\;=\;\frac{\epsilon}{2}\,\widehat{\sf x}.$$

$$x^{\prime }=\frac{ϵ}{2}\left(1-\frac{}{}{y}^2\right)+y,$$

$$y^{\prime \prime }=-\frac{ϵ}{2}-\left(1-\frac{{ϵ}^2}{2}\right)y+\frac{3}{2}ϵy^2-\frac{{ϵ}^2}{2}{y}^3.$$

$$(y^{\prime }{)}^2$$

$$y=\left(\frac{2}{ϵ}-1\right)-\frac{2}{ϵ}\left(\frac{1-ϵ}{1-ϵw^2}\right),$$

$$y(\phi ,ϵ)=\left(\frac{2}{ϵ}-1\right)-\frac{2}{ϵ}\left[\frac{1-ϵ}{1-ϵ{\mathrm{sn}}^2(\phi /2|{ϵ}^2)}\right]$$

$${\int }_0^s\frac{du}{1-ϵ{\mathrm{sn}}^2(u|{ϵ}^2)}=\Pi \left(ϵ,\mathrm{am}(s|{ϵ}^2)|{ϵ}^2\right),$$

$$\overline{y}(ϵ)$$

$$\frac{\Pi \left(ϵ,\pi |{ϵ}^2\right)}{\mathsf{K}({ϵ}^2)}=1+\frac{\pi }{2(1-ϵ)\mathsf{K}({ϵ}^2)}.$$

$$\theta (\phi ,ϵ)=\frac{\pi \phi }{2\mathsf{K}({ϵ}^2)}+\Delta \theta (\phi ,ϵ),$$

$$\Delta \theta (\phi ,ϵ)=(1-ϵ)\left[4\Pi \left(ϵ,\mathrm{am}(\phi /2)\right)-\frac{\Pi \left(ϵ,\pi \right)}{\mathsf{K}({ϵ}^2)}\phi \right]$$

$$y(\phi ,ϵ)=\frac{1}{ϵ}\left(1-\frac{\pi }{2\mathsf{K}({ϵ}^2)}\right)-\frac{1}{ϵ}\frac{\partial \Delta \theta }{\partial \phi },$$

$$x(\phi ,ϵ)={\int }_0^{\phi }\left[\frac{ϵ}{2}\left(1-\frac{}{}{y}^2(t,ϵ)\right)+y(t,ϵ)\right]dt,$$

$$x(\phi ,ϵ)$$

$$\mathcal{Z}\left(\mathrm{am}(\phi /2|{ϵ}^2)\frac{}{}|{ϵ}^2\right)\equiv \frac{\partial }{\partial \phi }\ln \left[{\vartheta }_4^2(\zeta ,q)\right],$$

$$\frac{2\mathrm{sn}\mathrm{cn}\mathrm{dn}}{1-ϵ{\mathrm{sn}}^2}=-\frac{2}{ϵ}\frac{\partial }{\partial \phi }\ln \left[1-ϵ{\mathrm{sn}}^2(\phi /2|{ϵ}^2)\right],$$

$$1-ϵ{\mathrm{sn}}^2(\phi /2|{ϵ}^2)\equiv 1-{\vartheta }_1^2(\zeta ,q)/{\vartheta }_4^2(\zeta ,q),$$

$$x(\phi ,ϵ)$$

$$v_{\mathrm{D}}(ϵ)\equiv \frac{x(4\mathsf{K},ϵ)}{4\mathsf{K}}=\frac{1}{ϵ}\left(1-\frac{\mathsf{E}({ϵ}^2)}{\mathsf{K}({ϵ}^2)}\right).$$

Full text

I Introduction

The guiding-center dynamics of charged particles moving in a nonuniform magnetic field plays a crucial role in our understanding of magnetically-confined laboratory and space plasmas Cary_Brizard_2009 ; Tronko_Brizard_2015 . The standard guiding-center ordering used in deriving these guiding-center equations of motion is based on the dimensionless small parameter $\epsilon\equiv\rho/L\ll 1$, where $\rho$ is the gyromotion length scale and $L$ is the magnetic-field gradient length scale Note . Because of its important extensions to the gyrokinetic self-consistent treatment of low-frequency fluctuations in magnetized plasmas Brizard_Hahm_2007 , it is necessary to gain a full understanding of the validity of the guiding-center approximation, especially when plasma gradients are strong (see, e.g., Refs. Dimits_2012 and Coury_2016 ). For example, in the pedestal region of advanced tokamak plasmas Coury_2016 , the gradient length scale can be as small as $L\simeq 1-2$ cm, which means that a 10 keV proton confined by a 5 T magnetic field (with a thermal gyroadius $\rho=2$ mm) is represented by $\epsilon\simeq 0.1-0.2$, which falls well outside the standard guiding-center limit $\epsilon\ll 1$.

In order to investigate the validity of the guiding-center approximation in the presence of strong gradients, we consider the simplest non-trivial problem of a charged particle (with mass $m$ and charge $e$) moving in a straight magnetic field with a uniform magnetic-field gradient:

[TABLE]

where $B_{0}$ is a constant and $L\equiv|\nabla\ln B|^{-1}$ is the constant gradient length scale.

The Lorentz-force equations for the perpendicular motion in the $(x,y)$-plane are

[TABLE]

where we introduced the dimensionless coordinates $(x,y)\equiv(x/\rho,y/\rho)$ and a prime denotes a derivative with respect to $\varphi\equiv\Omega\,t$ (with $\Omega=eB_{0}/mc$). Since the magnetic field is straight $(\widehat{{\sf b}}=\widehat{\sf z})$, we ignore parallel motion along the $z$-axis. We note that these equations satisfy the energy conservation law $x^{\prime 2}+y^{\prime 2}\equiv 2\mu\,B_{0}/(m\rho^{2}\Omega^{2})=1$, where $\mu$ denotes the lowest-order magnetic moment.

In the present paper, we will study the solution for Eqs. (2)-(3) for $0\leq\epsilon\leq 1$ (instead of the standard guiding-center limit $\epsilon\ll 1$ Cary_Brizard_2009 ), subject to the initial conditions

[TABLE]

Figure 1 shows the numerical solution of Eqs. (2)-(3), subject to the initial conditions (4), for the case $\epsilon=0.5$ (i.e., the gyromotion length scale $\rho=L/2$ is half of the magnetic-gradient length scale L). Here, we clearly see the standard grad-B drift motion along the $x$-axis due to the magnetic-field gradient (which is exaggerated, here, by choosing $\epsilon=0.5$).

In Fig. 1, we also note that, while the transverse motion along the $y$-axis is periodic in the range $-1\leq y\leq 1$, the average $y$-position is not zero for $\epsilon\neq 0$. Here, the orbital average (dashed horizontal line) of the particle’s $y$ position is evaluated as

[TABLE]

where the complete elliptic integral of the first kind ${\sf K}(m)$ is defined as Mathematica

[TABLE]

In the uniform limit (i.e., $\epsilon=0$), we find $\overline{y}(0)=0$, as expected for the trigonometric solution $y(\varphi,0)=\cos\varphi$, with period $4{\sf K}(0)=2\pi$.

The purpose of the present work is to compare the orbit-averaged properties of the exact analytical solution of Eqs. (2)-(3) for the perpendicular motion in the $(x,y)$-plane, which will be explicitly expressed in terms of Jacobi elliptic functions NIST_Jacobi and elliptic integrals, with the predictions of higher-order guiding-center theory Tronko_Brizard_2015 . With these exact orbit averages, we will show that the results of guiding-center theory are valid well outside the guiding-center limit $\epsilon\ll 1$, which is consistent with Mynick’s work Mynick_1979 ; Mynick_1980 , where guiding-center theory was extended to $\epsilon$ comparable to unity within the context of bounce-center theory.

We note that the Jacobi elliptic functions have also recently been used in the description of guiding-center and bounce-center dynamics in axisymmetric tokamak geometry Brizard_2011 ; Brizard_2014 ; Duthoit_2014 . They have also been used in obtaining exact analytical solutions for charged-particle motion in nonuniform magnetic fields represented by power laws Headland_Seymour_1975 ; Repko_1991 , which includes the case of particle orbits crossing a line where $B=0$ (case 3 in Ref. Headland_Seymour_1975 ), and stationary current distributions Essen_Nordmark_2016 . The case of a discontinuous magnetic field (with a gradient length scale $L=0$) is also solved exactly in Ref. Dodin_Fisch_2001 .

II Guiding-center Approximation

We begin by deriving the lowest-order guiding-center results to be quoted in the remainder of this paper. Because magnetic-curvature and parallel-gradient effects vanish for the magnetic field (1), we ignore the parallel motion along field lines.

Two features associated with the averaged properties of the particle motion are seen in Fig. 1: an averaged displacement $\overline{y}(\epsilon)$ along the $y$-axis (in the opposite direction to the magnetic-field gradient) and a drift motion $v_{\rm D}(\epsilon)\equiv x(4{\sf K},\epsilon)/(4{\sf K})$ along the $x$-axis. Both features can easily be explained by higher-order guiding-center theory Cary_Brizard_2009 ; Tronko_Brizard_2015 .

First, for the guiding-center approximation of the averaged transverse displacement $\overline{y}(\epsilon)$, we compute the guiding-center averaged particle position $\langle{\bf x}\rangle_{\rm gc}\equiv\langle{\sf T}_{\rm gc}{\bf X}\rangle_{\rm gc}$ according to high-order guiding-center theory Tronko_Brizard_2015 , where

[TABLE]

is defined as the guiding-center pull-back of the guiding-center position ${\bf X}=(X,Y)$ to particle phase space, where the phase-space vector fields $({\sf G}_{1},{\sf G}_{2},\cdots)$ generate the guiding-center transformation Tronko_Brizard_2015 . From this expression, we compute the guiding-center averaged-particle displacement

[TABLE]

where the gyroradius vector \mbox{\boldmath\rho}_{0}(Y,\mu,\theta)\equiv-\,G_{1}^{\bf X} depends on the gyroangle $\theta$, the guiding-center magnetic moment $\mu$, and the guiding-center coordinate $Y$ (because its magnitude depends on the strength of the magnetic field), and we used \langle\mbox{\boldmath\rho}_{0}\rangle_{\rm gc}=0 in Eq. (7). From Eq. (7), we see that the guiding-center-averaged $x$-position of the particle $\langle x\rangle_{\rm gc}=X$ is equal to the guiding-center position $X$. On the other hand, because of the magnetic-field gradient (with $\nabla\ln B=-\,\widehat{\sf y}/L$), the normalized guiding-center-averaged particle displacement along the $y$-axis is

[TABLE]

This guiding-center prediction is, in fact, quite close (see Fig. 3) to the orbit-averaged position (5) shown by the solid horizontal line in Fig. 1.

Next, the dimensionless guiding-center drift velocity $\dot{\bf X}/(\rho\Omega)$ in a straight magnetic field with a constant gradient (normalized to the perpendicular particle velocity scale $\rho\Omega$) is Cary_Brizard_2009

[TABLE]

This guiding-center result is clearly seen as the orbit-averaged drift motion along the $x$-axis in Fig. 1. The guiding-center prediction (9) is also quite close to the orbit-averaged drift velocity along the $x$-axis (see Fig. 5). We note that the case of a discontinuous magnetic field discussed in Ref. Dodin_Fisch_2001 exhibits a drift motion that can be calculated exactly along lines similar to the derivation of Eq. (9).

While these guiding-center predictions (8)-(9) are derived in the standard guiding-center limit $\epsilon\ll 1$ (i.e., when the magnetic-gradient length scale $L$ is assumed to be much longer than the gyromotion length scale $\rho$), we will show that these guiding-center predictions are, in fact, also valid well outside the limit $\epsilon\ll 1$.

III Exact Transverse Motion

We now derive the exact solution for the transverse motion along the $y$-axis, which is expressed in terms of the Jacobi elliptic functions NIST_Jacobi . With this solution, we will obtain an explicit expression for the orbit average (5), which can then be compared with the guiding-center predictions (8).

III.1 Exact solution

We begin by noting that, since the right side of Eq. (2) can be expressed as an exact derivative with respect to $\varphi$, we integrate it to obtain

[TABLE]

where we used the initial conditions (4). By inserting Eq. (10) into Eq. (3), we obtain the nonlinear second-order ODE

[TABLE]

Next, we multiply this equation by $2\,y^{\prime}$, and, then using the initial conditions (4), we integrate it to finally obtain

[TABLE]

We will now seek an exact solution $y(\varphi,\epsilon)$ of Eq. (11) for the transverse motion in terms of the Jacobi elliptic functions.

For this purpose, we transform Eq. (11) by writing

[TABLE]

where $w(\varphi,\epsilon)$ is a solution of the differential equation $(w^{\prime})^{2}=\frac{1}{4}\,(1-w^{2})\,(1-\epsilon^{2}\,w^{2})$. By using the initial condition $w(0,\epsilon)=0$, the solution $w(\varphi,\epsilon)={\rm sn}(\varphi/2|\epsilon^{2})$ is expressed in terms of the doubly-periodic Jacobi elliptic function ${\rm sn}(z|m)$, with a real period of $4{\sf K}(m)$ and an imaginary period $2i\,{\sf K}(1-m)$. Hence, the transverse motion

[TABLE]

corresponds exactly to the numerical solution shown in Fig. 1, which is periodic in the range $-1\leq y\leq 1$ with a period of $4\,{\sf K}(\epsilon^{2})$, as used in Eq. (5). Note that at the half-period $2\,{\sf K}$, we have the exact result $y(2\,{\sf K},\epsilon)=-1$, which holds for all values of $\epsilon$. This exact solution is shown in Fig. 2 for $\epsilon=0.1$ (solid) and $\epsilon=0.99$ (dashed). We note that, as $\epsilon$ increases toward 1, the solution spends an increasingly larger portion of its orbit in the range $y>0$, although it must still satisfy $y(2\,{\sf K},\epsilon)=-1$. This feature explains why the orbit average (5) is positive for finite magnetic gradients $\epsilon>0$.

III.2 Averaged transverse motion

We can now use the Jacobi solution (12) to obtain an explicit expression for the orbit average (5). For this purpose, we introduce the definition

[TABLE]

expressed in terms of the incomplete elliptic integral of the third kind, where ${\rm am}(s|m)\equiv\arcsin[{\rm sn}(s|m)]$ is the Jacobi amplitude function NIST_Jacobi . Hence, we compute the orbit average (5) of Eq. (12):

[TABLE]

where we used ${\rm am}(2{\sf K}|\epsilon^{2})=\pi$ and

[TABLE]

The plots of Eq. (14) (solid curve) and the lowest-order guiding-center approximation $\epsilon/4$ (dashed curve), given by Eq. (8), are shown in Fig. 3. We note that the lowest-order guiding-center approximation (8) for the averaged particle $y$-position is remarkably close to the orbit averaged position (14) for values of $\epsilon$ well beyond the guiding-center limit $\epsilon\ll 1$. In Fig. 1, for example, we find $\overline{y}(0.5)=0.136...$, and the guiding-center result $\epsilon/4=0.125$, which is just 8% below the orbit-averaged value $\overline{y}(\epsilon)$.

Lastly, we note that an alternative way to write the transverse solution (12) involves the gyroangle $\theta(\varphi,\epsilon)$, which satisfies the differential equation $\theta^{\prime}=1-\epsilon\,y(\varphi,\epsilon)$. Using Eq. (12), its solution is expressed as

[TABLE]

where the periodic function

[TABLE]

vanishes at $\varphi=0,2{\sf K}$, and $4{\sf K}$. The transverse solution (12) can, therefore, be expressed as

[TABLE]

where the orbit average of the second term is explicitly zero because of the periodicity of Eq. (16).

IV Drift Motion

The solution for $x(\varphi,\epsilon)$ can now be obtained by integrating Eq. (10):

[TABLE]

subject to the initial conditions (4). By inserting Eq. (12), we find the solution

[TABLE]

which is expressed in terms of the Jacobi elliptic functions $({\rm sn},{\rm cn},{\rm dn})$, the complete elliptic integral of the second kind, ${\sf E}(\epsilon^{2})$, and the Jacobi zeta function NIST_Jacobi

[TABLE]

which is expressed in terms of the logarithmic derivative of the elliptic theta function $\vartheta_{4}(\zeta,q)$ NIST_Theta , where $\zeta=\pi\varphi/(4{\sf K}(\epsilon^{2}))$ and $q(\epsilon)=\exp[-\pi\;{\sf K}(1-\epsilon^{2})/{\sf K}(\epsilon^{2})]$. Since the last term in Eq. (19) can also be expressed as a logarithmic derivative

[TABLE]

we use the identity NIST_Theta

[TABLE]

so that Eq. (19) can also be expressed as

[TABLE]

where the second term is a periodic function of $\zeta$. Once again, this solution, which is shown in Fig. 4, agrees exactly with the numerical solution shown in Fig. 1. It is interesting to note that an expression similar to the second term in Eq. (IV) has previously appeared in Ref. Brizard_2014 , where the generating functions for the canonical transformation for trapped/passing guiding-center orbits in axisymmetric tokamak geometry were derived.

The dashed line in Fig. 4 represents the drift motion along the $x$-axis, with a slope equal to the orbit-averaged particle drift velocity obtained from Eq. (IV) :

[TABLE]

Equation (21) and the lowest-order guiding-center prediction $\epsilon/2$, derived in Eq. (9), are shown in Fig. 5. Once again, the lowest-order guiding-center prediction (9) yields an excellent agreement with the exact result (21) well outside the standard guiding-center limit $\epsilon\ll 1$. For example, for the case $\epsilon=0.5$, the guiding-center prediction $\epsilon/2$ yields $0.5/2=0.25$, which is just 3% below the orbit-averaged particle drift velocity $v_{\rm D}(0.5)=0.259...$

V Summary

By solving exactly the motion of a charged particle in a straight magnetic field (1) with a constant gradient, we have been able to investigate the validity of the guiding-center approximation. In the present work, we have shown that the guiding-center predictions (8) and (9) agree very well with the orbit-averaged particle displacement (14) and the orbit-averaged particle drift velocity (21) for values of $\epsilon$ well outside the standard guiding-center limit $\epsilon\ll 1$.

These results can, thus, be used to justify the applications of gyrokinetic theory for magnetized plasmas with strong gradients Dimits_2012 . For example, in the pedestal region of advanced tokamak plasmas Coury_2016 , where $\epsilon\simeq 0.2$, the guiding-center predictions (8) and (9) are 1.3% and 0.5% below the orbit-averaged particle displacement (14) and the orbit-averaged particle drift velocity (21), respectively.

Acknowledgements.

The author wishes to thank H. Essén for pointing out Refs. Headland_Seymour_1975 and Essen_Nordmark_2016 . This work was supported by the U.S. Dept. of Energy under contract No. DE-SC0014032.