Continuous-stage adapted exponential methods for charged-particle dynamics with arbitrary magnetic fields

Abstract

This paper is devoted to the numerical symplectic approximation of the charged-particle dynamics (CPD) with a homogeneous magnetic field and its extension to a non-homogeneous magnetic field. By utilizing continuous-stage methods and exponential integrators, a general class of symplectic methods is formulated for CPD under a homogeneous magnetic field. Based on the derived symplectic conditions, two practical symplectic methods up to order four are constructed where the error estimates show that the proposed second order scheme has a uniform accuracy in the position w.r.t. the strength of the magnetic field. Moreover, the symplectic methods are extended to CPD under a non-homogeneous magnetic field and three algorithms are formulated. Rigorous error estimates are investigated for the proposed methods and one method is proved to have a uniform accuracy in the position w.r.t. the strength of the magnetic field. Numerical experiments are provided for CPD under homogeneous and non-homogeneous magnetic fields, and the numerical results support the theoretical analysis and demonstrate the remarkable numerical behavior of our methods.

Appendix

Firstly, consider a Taylor expansion of the function $\vartheta (\nu )$ at the origin with real coefficients $c_{n}$, and formulate $\vartheta (iy)$ as the following form

$$\begin{aligned} \vartheta (iy)=\vartheta (0)+iy\vartheta _{1}(y)-y^{2}\vartheta _{2}(y),\end{aligned}$$

(1.1)

where $\vartheta _{1}(y)=\sum \limits _{j \ge 0} c_{2j+1}(-y^2)^{j}$ and $\vartheta _{2}(y)=\sum \limits _{j \ge 0} c_{2j+2}(-y^2)^{j}$.

Then for a vector $B=\big (B_1,B_2,B_3\big )^{\intercal }$ and a skew symmetric matrix $M=\frac{1}{\epsilon }\left( \begin{array}{ccc} 0 &{} B_3 &{} -B_2 \\ -B_3 &{} 0 &{} B_1 \\ B_2 &{} -B_1 &{} 0 \\ \end{array} \right)$, it can be checked that $M^3=-m^2 M$ with $m= \frac{\Vert B\Vert }{\epsilon }$. Based on (5.1), one obtains that

$$\vartheta (M)=\vartheta (0) I+\vartheta _{1}(m) M+\vartheta _{2}(m) M^2,$$

which leads to the calculation of $\vartheta (M) v$ by estimating the scalars $\vartheta (0)$, $\vartheta _{1}(m)$, $\vartheta _{2}(m)$, and by forming twice a product of M with a vector.

Finally, we present the computation of $\varphi _{0}$ and $\varphi _{1}$ via the above results. Before going further, it is noted that $\vartheta _{1}(y)=0$ if $\vartheta (\nu )$ has only even powers of $\nu$ and $\vartheta (0)=0,\ \vartheta _{2}(y)=0$ if there are only odd powers of $\nu$. Now we are in a position to get

$$\begin{aligned} \begin{aligned}&\varphi _{0}(hM)=I+\frac{sin(hm)}{m}M+\frac{1-cos(hm)}{m^2}M^2,\\&\varphi _{1}(hM)=I+\frac{1-cos(hm)}{hm^2}M+\frac{1-sinc(hm)}{m^2}M^2. \end{aligned} \end{aligned}$$

The expressions of the coefficients $\alpha _{\tau \sigma }(hM)$, $\beta _{\tau }( hM )$ and $\gamma _{\tau }( hM )$ can be derived similarly. Moreover, from the above analysis, it follows that $\bar{\varphi }(W)=\varphi (\bar{W}),$ where $W=h\tilde{\Lambda }\textrm{i}$ with $\tilde{\Lambda }=diag (-m,0,m)$.