Exponential energy-preserving methods for charged-particle dynamics in a strong and constant magnetic field

TL;DR

This paper introduces exponential energy-preserving numerical methods for charged-particle dynamics in strong magnetic fields, ensuring exact energy conservation and near conservation of magnetic moment over long times, with proven properties and numerical validation.

Contribution

The paper develops and analyzes a novel exponential energy-preserving method that exactly conserves energy and nearly conserves magnetic moment for charged-particle simulations.

Findings

Exact energy preservation in the numerical method.

Near conservation of magnetic moment over long times.

Numerical experiments demonstrating long-term stability.

Abstract

In this paper, exponential energy-preserving methods are formulated and analysed for solving charged-particle dynamics in a strong and constant magnetic field. The resulting method can exactly preserve the energy of the dynamics. Moreover, it is shown that the magnetic moment of the considered system is nearly conserved over a long time along this exponential energy-preserving method, which is proved by using modulated Fourier expansions. Other properties of the method including symmetry and convergence are also studied. An illustrated numerical experiment is carried out to demonstrate the long-time behaviour of the method.

Equations

$$\begin{array}[c]{ll}\ddot{x}=\dot{x}\times\frac{1}{\epsilon}B+F(x),\quad x(t_{0})=x_{0},\quad\dot{x}(t_{0})=\dot{x}_{0},\ \ t\in[t_{0},T],\end{array}$$

$$E(x,v)=\frac{1}{2}{\left|v\right|}^2+U(x),$$

$$\tilde{B}=\left(\begin{array}[]{ccc}0&B_{3}&-B_{2}\\ -B_{3}&0&B_{1}\\ B_{2}&-B_{1}&0\\ \end{array}\right).$$

$$\begin{array}[c]{ll}\dot{x}=v,\ \dot{v}=\frac{1}{\epsilon}\tilde{B}v+F(x).\end{array}$$

$$I(x,v)=\frac{{\left|v_{\perp }\right|}^2}{2\left|B\right|}=\frac{1}{2{\left|\tilde{B}\right|}^3}{\left|\tilde{B}v\right|}^2$$

$$\left\{\begin{array}[c]{ll}x_{n+1}=x_{n}+h\varphi_{1}(\frac{h}{\epsilon}\tilde{B})v_{n}+h^{2}\varphi_{2}(\frac{h}{\epsilon}\tilde{B})\int_{0}^{1}F\big{(}x_{n}+\sigma(x_{n+1}-x_{n})\big{)}d\sigma,\\ v_{n+1}=e^{\frac{h}{\epsilon}\tilde{B}}v_{n}+h\varphi_{1}(\frac{h}{\epsilon}\tilde{B})\int_{0}^{1}F\big{(}x_{n}+\sigma(x_{n+1}-x_{n})\big{)}d\sigma,\end{array}\right.$$

$${\phi }_0(z)=e^z,{\phi }_k(z)={\int }_0^1{e}^{(1-\sigma )z}\frac{{\sigma }^{k-1}}{(k-1)!}d\sigma ,k=1,2.$$

$$\left\{\begin{array}[c]{ll}x_{n}=x_{n+1}-h\varphi_{1}(-\frac{h}{\epsilon}\tilde{B})v_{n+1}+h^{2}\varphi_{2}(-\frac{h}{\epsilon}\tilde{B})\int_{0}^{1}F\big{(}x_{n}+\sigma(x_{n+1}-x_{n})\big{)}d\sigma,\\ v_{n}=e^{-\frac{h}{\epsilon}\tilde{B}}v_{n+1}-h\varphi_{1}(-\frac{h}{\epsilon}\tilde{B})\int_{0}^{1}F\big{(}x_{n}+\sigma(x_{n+1}-x_{n})\big{)}d\sigma,\end{array}\right.$$

$$\int_{0}^{1}F\big{(}x_{n+1}+\sigma(x_{n}-x_{n+1})\big{)}d\sigma=\int_{0}^{1}F\big{(}x_{n}+\sigma(x_{n+1}-x_{n})\big{)}d\sigma.$$

$$v_{n+1}=e^{\frac{h}{ϵ}\tilde{B}}{v}_n+h{e}^{\frac{h}{ϵ}\tilde{B}}{\phi }_1(-\frac{h}{ϵ}\tilde{B})\mathcal{I}=e^{\frac{h}{ϵ}\tilde{B}}{v}_n+h{\phi }_1(\frac{h}{ϵ}\tilde{B})\mathcal{I}$$

$$\begin{array}[]{rl}x_{n+1}&=x_{n}+h\varphi_{1}(-\frac{h}{\epsilon}\tilde{B})v_{n+1}-h^{2}\varphi_{2}(-\frac{h}{\epsilon}\tilde{B})\mathcal{I}\\ &=x_{n}+h\varphi_{1}(-\frac{h}{\epsilon}\tilde{B})e^{\frac{h}{\epsilon}\tilde{B}}v_{n}+h^{2}\big{(}\varphi_{1}(-\frac{h}{\epsilon}\tilde{B})\varphi_{1}(\frac{h}{\epsilon}\tilde{B})-\varphi_{2}(-\frac{h}{\epsilon}\tilde{B})\big{)}\mathcal{I}\\ &=x_{n}+h\varphi_{1}(\frac{h}{\epsilon}\tilde{B})v_{n}+h^{2}\varphi_{2}(\frac{h}{\epsilon}\tilde{B})\mathcal{I}.\end{array}$$

$$x_0=\mathcal{O}(1),v_0=\mathcal{O}(1)$$

$$|\sin (\frac{1}{2}kh\tilde{\omega })|\ge c\sqrt{h}\mathrm{for}k=1,2,\dots ,N\mathrm{with}N\ge 2,$$

$$E(x_{n+1},v_{n+1})=E(x_n,v_n)\text{for}n=0,1,\dots .$$

$$\displaystyle I(x_{n},v_{n})=I(x_{0},v_{0})+\mathcal{O}\Big{(}\frac{1}{\left|\cos(\frac{1}{2}h\tilde{\omega})\right|}h\Big{)}+\mathcal{O}(h),$$

$$I(x_n,v_n)=I(x_0,v_0)+\mathcal{O}(h)$$

$$\begin{array}[c]{ll}E(x_{n+1},v_{n+1})=\dfrac{1}{2}v_{n+1}^{\intercal}v_{n+1}+U(x_{n+1}).\end{array}$$

$$(e^{\frac{h}{ϵ}\tilde{B}}{)}^{⊺}=e^{-\frac{h}{ϵ}\tilde{B}},({\phi }_1(\frac{h}{ϵ}\tilde{B}){)}^{⊺}={\phi }_1(-\frac{h}{ϵ}\tilde{B}),({\phi }_2(\frac{h}{ϵ}\tilde{B}){)}^{⊺}={\phi }_2(-\frac{h}{ϵ}\tilde{B}).$$

$$\displaystyle E(x_{n+1},v_{n+1})=\frac{1}{2}\big{(}e^{\frac{h}{\epsilon}\tilde{B}}v_{n}+h\varphi_{1}(\frac{h}{\epsilon}\tilde{B})\mathcal{I}\big{)}^{\intercal}\big{(}e^{\frac{h}{\epsilon}\tilde{B}}v_{n}+h\varphi_{1}(\frac{h}{\epsilon}\tilde{B})\mathcal{I}\big{)}+U(x_{n+1})$$

$$=$$

$$E(x_{n+1},v_{n+1})=\frac{1}{2}{v}_n^{⊺}{v}_n+h{\mathcal{I}}^{⊺}{\phi }_1(\frac{h}{ϵ}\tilde{B})v_n+\frac{1}{2}{h}^2{\mathcal{I}}^{⊺}{\phi }_1(-\frac{h}{ϵ}\tilde{B}){\phi }_1(\frac{h}{ϵ}\tilde{B})\mathcal{I}+U(x_{n+1}).$$

$$U(x_n)-U(x_{n+1})=-{\int }_0^{1}dU((1-\tau )x_n+\tau x_{n+1})$$

$$=$$

$$E(x_{n+1},v_{n+1})=$$

$$-h{\mathcal{I}}^{⊺}{\phi }_1(\frac{h}{ϵ}\tilde{B})v_n-h^2{\mathcal{I}}^{⊺}{\phi }_2(\frac{h}{ϵ}\tilde{B})\mathcal{I}$$

$$=$$

$$\displaystyle\Big{(}\varphi_{1}(-\frac{h}{\epsilon}\tilde{B})\varphi_{1}(\frac{h}{\epsilon}\tilde{B})-2\varphi_{2}(\frac{h}{\epsilon}\tilde{B})\Big{)}^{\intercal}=\sum\limits_{k=1}^{\infty}c_{k}\big{(}-\frac{h}{\epsilon}\tilde{B}\big{)}^{2k+1}$$

$$\displaystyle=-\sum\limits_{k=1}^{\infty}c_{k}\big{(}\frac{h}{\epsilon}\tilde{B}\big{)}^{2k+1}=-\Big{(}\varphi_{1}(-\frac{h}{\epsilon}\tilde{B})\varphi_{1}(\frac{h}{\epsilon}\tilde{B})-2\varphi_{2}(\frac{h}{\epsilon}\tilde{B})\Big{)}.$$

$$\displaystyle\mathcal{I}^{\intercal}\Big{(}\varphi_{1}(-\frac{h}{\epsilon}\tilde{B})\varphi_{1}(\frac{h}{\epsilon}\tilde{B})-2\varphi_{2}(\frac{h}{\epsilon}\tilde{B})\Big{)}\mathcal{I}$$

$$\displaystyle=-\mathcal{I}^{\intercal}\Big{(}\varphi_{1}(-\frac{h}{\epsilon}\tilde{B})\varphi_{1}(\frac{h}{\epsilon}\tilde{B})-2\varphi_{2}(\frac{h}{\epsilon}\tilde{B})\Big{)}\mathcal{I},$$

$$E(x_{n+1},v_{n+1})=\frac{1}{2}{v}_n^{⊺}{v}_n+U(x_n)=E(x_n,v_n).$$

$$\tilde{x}(t)=P^{\text{H}}x(t),\tilde{v}(t)=P^{\text{H}}v(t),$$

$$\begin{array}[c]{ll}\dot{\tilde{x}}=\tilde{v},&\tilde{x}_{0}=P^{\textup{H}}x_{0},\\ \dot{\tilde{v}}=\mathrm{i}\tilde{\Omega}\tilde{v}+\tilde{F}(\tilde{x}),&\tilde{v}_{0}=P^{\textup{H}}v_{0},\end{array}$$

Full text

Keywords: charged-particle dynamics, exponential energy-preserving methods, modulated Fourier expansions, long-time conservation, strong and constant magnetic field

MSC: 65L06, 65P10, 78A35, 78M25.

1 Introduction

In this paper, we derive and analyse energy-preserving methods for charged-particle dynamics in a strong and constant magnetic field

[TABLE]

where $x(t)\in{\mathbb{R}}^{3}$ describes the position of a particle, $B=\nabla_{x}\times A(x)$ is a constant magnetic field with the vector potential $A(x)=-\frac{1}{2}x\times B\in{\mathbb{R}}^{3}$ and $F(x)=-\nabla_{x}U(x)$ is an electric field with the scalar potential $U(x)$. Following [15], we are devoted to the situation of a small positive acaling parameter $0<\epsilon\ll 1$ and assume that $\left|B\right|\geq 1$ in the Euclidean norm. The energy of this dynamics is given by

[TABLE]

where $v=\dot{x}$ is the velocity of the particle. Denote the constant vector $B$ by $B=(B_{1},B_{2},B_{3})^{\intercal}$ with $B_{i}\in{\mathbb{R}}$ for $i=1,2,3.$ By the definition of the cross product, we obtain $\dot{x}\times B=\tilde{B}\dot{x},$ where $\tilde{B}$ is a skew symmetric matrix

[TABLE]

The system (1) as well as $v=\dot{x}$ can be rewritten as

[TABLE]

With the analysis given in [1, 4], it is well known that the magnetic moment

[TABLE]

is an adiabatic invariant, where $v_{\perp}=\frac{v\times B}{\left|B\right|}$ is orthogonal to $B$ and $|\tilde{B}|=\sqrt{B_{1}^{2}+B_{2}^{2}+B_{3}^{2}}.$ Recently, it has been shown in [15] that this quantity is nearly conserved over long time scales, which is proved by using the technique of modulated Fourier expansion with state-dependent frequencies and eigenvectors.

Charged-particle dynamics have been received much attention for a long time (see, e.g. [1, 4]) and many effective methods have been developed for solving this system. The Boris method was presented in [2] and it was researched further in [8, 13, 26]. Various other kinds of methods have also been researched for charged-particle dynamics, such as volume-preserving algorithms in [17], symplectic or K-symplectic algorithms in [18, 29, 33, 35], and symmetric multistep methods in [14]. It is worth mentioning that, more recently, a variational integrator has been studied in [15] for solving charged-particle dynamics in a strong magnetic field. In this paper we are interested in the formulation and analysis of energy-preserving (EP) methods when applied to charged-particle dynamics in a strong and constant magnetic field.

For the energy-preserving methods, there have been a lot of studies on this topic. Many different EP methods have been presented and analysed, such as the average vector field (AVF) method (see, e.g. [5, 27]), discrete gradient methods (see, e.g. [23]), Hamiltonian Boundary Value Methods (see, e.g. [3]), EP collocation methods (see, e.g. [10]) and trigonometric/exponential EP methods (see, e.g. [21, 31, 34]). Based on these work, we will derive and analyse novel EP integrators for charged-particle dynamics in a strong and constant magnetic field. The long time magnetic moment conservation of the new integrator will also be researched via its modulated Fourier expansion. The technique of modulated Fourier expansion was firstly given in [11] and then it has been successfully used in the study of long-time behaviour for numerical methods/differential equations (see, e.g. [6, 7, 9, 12, 16, 28]).

The rest of this paper is organized as follows. In Section 2, we first formulate the scheme of the method and then prove that it is symmetric. In Section 3, two main results concerning energy preservation and magnetic moment preservation are presented and a numerical experiment is reported to support these results. The proofs of the two main results are given in Sections 4-5, respectively. The concluding remarks of this paper are given in the last section.

2 Formulation of the method

In order to effectively solve the system (3), we consider the following method.

Definition 2.1

The exponential energy-preserving method for solving (3) is defined by:

[TABLE]

where $h$ is a stepsize and the $\varphi$-functions are defined by

[TABLE]

We denote this method by EEP.

Remark 2.2

It is noted that this kind of method belongs to exponential integrators, which have been widely developed and researched for solving highly oscillatory systems (see, e.g. [19, 20, 21, 24, 30, 32]).

Proposition 2.3

The EEP integrator (5) is symmetric.

Proof Exchanging $(x_{n},v_{n})\leftrightarrow(x_{n+1},v_{n+1})$ and $h\leftrightarrow-h$ in (5) yields

[TABLE]

where the following fact has been used

[TABLE]

We thus obtain

[TABLE]

by considering the second formula of (7). Inserting this into the first formula of (7), we get

[TABLE]

Therefore, the energy-preserving integrator (5) is symmetric.

3 Main results and numerical experiment

3.1 Main results

Before presenting the main results of this paper, we need the following assumptions, which have been considered in [11, 16].

Assumption 3.1

• •

We consider the initial values

[TABLE]

such that the energy $E$ is bounded independently of $\epsilon$ along the solution.

• •

It is assumed that the numerical solution stays in a compact set.

• •

We require a lower bound on the stepsize $h\tilde{\omega}\geq c_{0}>0$ with $\tilde{\omega}=\frac{|\tilde{B}|}{\epsilon}$.

• •

The following numerical non-resonance condition is assumed to be true

[TABLE]

which imposes a restriction on $N$ for a given $h$ and $\tilde{\omega}$.

We first give the result about the energy preservation of the EEP method.

Theorem 3.2

(Energy preservation.) The EEP integrator (5) preserves the energy $E$ in (2) exactly, i.e.,

[TABLE]

Now, we present the result about the long time magnetic moment conservation of the EEP method.

Theorem 3.3

(Magnetic moment conservation.) Under the conditions of Assumption 3.1, the long time magnetic moment is nearly conserved over long times by the EEP method

[TABLE]

where $0\leq nh\leq h^{-N+1}.$ The constants symbolized by $\mathcal{O}$ depend on $N$, the final time $T$ and the constants in the assumptions, but are independent of $n,\ h,\ \epsilon$.

Remark 3.4

We note that it is easy to make a choice of the stepsize $h$ such that $\left|\cos(\frac{1}{2}h\tilde{\omega})\right|$ is not small. In this case, an improved result concerning the long time magnetic moment conservation is obtained

[TABLE]

for $0\leq nh\leq h^{-N+1}.$

3.2 Numerical experiment

As an illustrative numerical experiment, we consider the charged particle system of [14] with a constant magnetic field and an additional factor $1/\epsilon$. The system can be given by (1) with the potential $U(x)=\frac{1}{100\sqrt{x_{1}^{2}+x_{2}^{2}}}$ and the constant magnetic field $B=(0,0,1)^{\intercal}.$ The initial values are chosen as $x(0)=(0.7,1,0.1)^{\intercal}$ and $v(0)=(0.9,0.5,0.4)^{\intercal}.$

We consider applying four-point Gauss-Legendre’s rule to the integral of the EEP integrator (5). The fixed-point iteration is chosen here and we set $10^{-16}$ as the error tolerance and $50$ as the maximum number of each iteration. We choose $\epsilon=0.01,0.0001$. This problem is integrated on $[0,10000]$ with $h=0.01$ and see Figures 1-2 for the relative errors $(H(x_{n},v_{n})-H(x^{0},v^{0}))/H(x^{0},v^{0})$ of the energy and $(I(x_{n},v_{n})-I(x^{0},v^{0}))/I(x^{0},v^{0})$ of the magnetic moment, respectively. In order to show the performance of our EEP method, we choose Boris method for comparison. We consider $\epsilon=0.0001$ and solve this system on $[0,10000]$ with $h=0.01$. The results of Boris method are shown in Figure 3. Finally, the problem is solved with $\epsilon=0.05,$ $T=10,100,1000$ and $h=1/(50\times 2^{i})$ for $i=0,\ldots,3$. The global errors are shown in Figure 4.

From the results, it can be observed that our EEP method shows an excellent energy-preserving property, a prominent long-term behavior in the numerical magnetic moment conservation and a good accuracy. All these observations support the theoretical results given in Theorems 3.2- 3.3.

4 Proof of energy preservation

In this section, we give the proof of Theorem 3.2.

Proof  In this paper, we let \mathcal{I}:=\int_{0}^{1}F\big{(}x_{n}+\sigma(x_{n+1}-x_{n})\big{)}d\sigma for brevity. We firstly compute

[TABLE]

Keeping the fact in mind that $\tilde{B}$ is skew-symmetric, one obtains that

[TABLE]

Inserting the second formula of (5) into (9) with some manipulation yields

[TABLE]

It follows from (6) that $\varphi_{1}(-\frac{h}{\epsilon}\tilde{B})e^{\frac{h}{\epsilon}\tilde{B}}=\varphi_{1}(\frac{h}{\epsilon}\tilde{B}).$ Therefore, we obtain

[TABLE]

On the other hand, it can be checked that

[TABLE]

where the first equation of (5) has been used. Inserting this result into (10) implies

[TABLE]

It follows from the definition (6) that \varphi_{1}(-\frac{h}{\epsilon}\tilde{B})\varphi_{1}(\frac{h}{\epsilon}\tilde{B})-2\varphi_{2}(\frac{h}{\epsilon}\tilde{B})=\sum\limits_{k=1}^{\infty}c_{k}\big{(}\frac{h}{\epsilon}\tilde{B}\big{)}^{2k+1} with the coefficients $c_{k}$ for $k=1,2,\ldots$, which gives that

[TABLE]

Based on this result, one arrives at

[TABLE]

which shows that \mathcal{I}^{\intercal}\Big{(}\varphi_{1}(-\frac{h}{\epsilon}\tilde{B})\varphi_{1}(\frac{h}{\epsilon}\tilde{B})-2\varphi_{2}(\frac{h}{\epsilon}\tilde{B})\Big{)}\mathcal{I}=0. Therefore, (11) becomes

[TABLE]

5 Proof of the magnetic moment conservation

This section is devoted to the proof of Theorem 3.3. Modulated Fourier expansion will be used here. We will first derive the modulated Fourier expansion for the EEP method in Subsection 5.1 and then show an almost-invariant of the expansion in Subsection 5.2. Based on these analysis, Theorem 3.3 will be proved immediately.

Since the matrix $\tilde{B}$ is skew-symmetric, there exists a unitary matrix $P$ and a diagonal matrix $\Lambda$ such that $\tilde{B}=P\Lambda P^{\textup{H}}$ with $\Lambda=\textmd{diag}(-|\tilde{B}|\mathrm{i},0,|\tilde{B}|\mathrm{i})$. By the linear change of variable

[TABLE]

we rewrite the system (3) as

[TABLE]

where $\tilde{\Omega}=\textmd{diag}(-\tilde{\omega},0,\tilde{\omega})$ with $\tilde{\omega}=\frac{|\tilde{B}|}{\epsilon}$ and $\tilde{F}(\tilde{x})=P^{\textup{H}}F(P\tilde{x})=-\nabla_{\tilde{x}}U(P\tilde{x})$. In this paper, a vector $x$ in ${\mathbb{R}}^{3}$ or ${\mathbb{C}}^{3}$ is denoted by $x=(x_{-1},x_{0},x_{1})^{\intercal}$. For the transformed system (13), its magnetic moment has the following form

[TABLE]

The EEP integrator (5) for solving this transformed system is defined as

[TABLE]

In this section, the following five operators will be used:

[TABLE]

where $D$ is the differential operator (see [16]). We have the following properties of these operators.

Proposition 5.1

The Taylor expansions of the operator $L(hD)$ are expressed by

[TABLE]

where

[TABLE]

The operator $(L_{3}L_{4}^{-1})(hD)$ has the following Taylor expansions

[TABLE]

Moreover, for the operator $L_{5}(hD,\tau,k)$ and $\left|k\right|>0$, it is true that

[TABLE]

5.1 Modulated Fourier expansion

A modulated Fourier expansion of the EEP method (15) for solving the transformed system (13) is given as follows.

Theorem 5.2

Assume that all the conditions of Assumption 3.1 hold. The numerical solution given by the EEP integrator (15) admits a modulated Fourier expansion

[TABLE]

where $0\leq t=nh\leq T$, $N$ is a fixed integer such that (8) holds and the remainder terms are bounded by

[TABLE]

The bounds of the coefficient functions of $\tilde{\zeta}^{k}$ as well as all their derivatives are given by

[TABLE]

and further by

[TABLE]

The bounds of the coefficient functions of $\tilde{\eta}^{k}$ as well as all their derivatives are

[TABLE]

By noticing $P\tilde{\zeta}^{-k}=\bar{P}\overline{\tilde{\zeta}^{k}},$ $P\tilde{\eta}^{-k}=\bar{P}\overline{\tilde{\eta}^{k}}$ and P^{\textup{H}}\bar{P}=\left(\begin{array}[]{ccc}0&0&1\\ 0&1&0\\ 1&0&0\\ \end{array}\right), one obtains that $\tilde{\zeta}_{-l}^{-k}=\overline{\tilde{\zeta}_{l}^{k}}$ and $\tilde{\eta}_{-l}^{-k}=\overline{\tilde{\eta}_{l}^{k}}$. The constants symbolized by the notation depend on the constants from Assumption 3.1 and the final time $T$, but are independent of $h$ and $\tilde{\omega}$.

Proof  We will construct the functions

[TABLE]

with smooth coefficient functions and prove that there is only a small defect when $\tilde{x}_{h}$ and $\tilde{v}_{h}$ are inserted into the numerical scheme (15).

$\bullet$ Construction of the coefficients functions. For the term $(1-\tau)\tilde{x}_{n}+\tau\tilde{x}_{n+1}$, we consider the function

[TABLE]

as its modulated Fourier expansion. Then in the light of (22), we get

[TABLE]

which yields

[TABLE]

By eliminating $\mathcal{I}$ in (15), one gets

[TABLE]

Inserting (22) into (24) and comparing the coefficients of $\mathrm{e}^{\mathrm{i}k\tilde{\omega}t}$ yields

[TABLE]

On the basis of this formula, we can obtain

[TABLE]

which will be used in the analysis of the next subsection.

Based on the second formula of (15) and (25), one has

[TABLE]

where the sum ranges over $m\geq 0$, $\alpha=(\alpha_{1},\ldots,\alpha_{m})$ with integer $\alpha_{i}$ satisfying $0<|\alpha_{i}|<N$, $s(\alpha)=\sum\limits_{j=1}^{m}\alpha_{j},$ and $(\tilde{\xi}(t,\tau))^{\alpha}$ is an abbreviation for $(\tilde{\xi}^{\alpha_{1}}(t,\tau),\ldots,\tilde{\xi}^{\alpha_{m}}(t,\tau))$. This formula as well as (23) presents the modulation system for the coefficients $\tilde{\zeta}^{k}(t)$ of the modulated Fourier expansion $q_{n}$. By using the results given in Proposition 5.1 and choosing the dominate terms in the relations yields the ansatz of $\tilde{\zeta}^{k}(t)$:

[TABLE]

Similarly, the ansatz of the modulated Fourier functions $\tilde{\eta}^{k}(t)$ can be obtained by considering the dominating terms of (25) and (27).

$\bullet$ Initial values.

According to the conditions $x_{h}(0)=\tilde{x}_{0}$ and $v_{h}(0)=\tilde{v}_{0}$, we obtain

[TABLE]

Thus the initial values $\tilde{\zeta}_{0}^{0}(0)=\mathcal{O}(1)$ and $\dot{\tilde{\zeta}}_{0}^{0}(0)=\mathcal{O}(1)$ can be arrived by considering the first and third formulae. It follows from the fourth formula that \tilde{\zeta}_{1}^{1}(0)=\frac{1}{\mathrm{i}\tilde{\omega}}\big{(}\tilde{v}_{1}^{0}-\dot{\tilde{\zeta}}_{1}^{0}(0)+\mathcal{O}(\epsilon)\big{)}=\mathcal{O}(\epsilon) and similarly one has $\tilde{\zeta}_{-1}^{-1}(0)=\mathcal{O}(\epsilon)$. Then one gets the initial value $\tilde{\zeta}_{\pm 1}^{0}(0)=\mathcal{O}(1)$ in the light of (28).

$\bullet$ Bounds of the coefficients functions.

The bounds (19) can be obtained by considering the initial values obtained above, the ansatz and Assumption 3.1. On the basis of the initial values obtained above, we get the bounds (20). The bounds given in (21) are true in the light of (26).

$\bullet$ Defect.

The defect (18) can be obtained by using the Lipschitz continuous of the nonlinearity and the standard convergence estimates.

Theorem 5.3

The numerical solution given by the EEP integrator (5) for solving (3) admits the following modulated Fourier expansion

[TABLE]

where $\zeta^{k}(t)=P\tilde{\zeta}^{k}(t)$ and $\eta^{k}(t)=P\tilde{\eta}^{k}(t).$ This relation yields $\zeta^{-k}=\overline{\zeta^{k}}$ and $\eta^{-k}=\overline{\eta^{k}}$. The bounds of the remainders $R_{h,N}$ and $S_{h,N}$ are the same as those given in Theorem 5.2.

5.2 An almost-invariant

With the coefficient functions constructed in last subsection, we let

[TABLE]

and

[TABLE]

An almost-invariant of the modulated Fourier expansion (29) is given as follows.

Theorem 5.4

Under the conditions of Theorem 5.2, there exists a function $\widehat{\mathcal{M}}[\vec{\tilde{\zeta}},\vec{\tilde{\eta}}]$ such that

[TABLE]

for $0\leq t\leq T,$ where

[TABLE]

Proof  With Theorem 5.2, we obtain that

[TABLE]

where we have used the denotations $\tilde{x}_{h}(t)=\sum\limits_{|k|<N}\tilde{x}^{k}_{h}(t),\ \tilde{q}_{h}(t,\tau)=\sum\limits_{|k|<N}\tilde{q}^{k}_{h}(t,\tau)$ with $\tilde{x}^{k}_{h}(t)=\mathrm{e}^{\mathrm{i}k\tilde{\omega}t}\tilde{\zeta}^{k}(t)$ and $\tilde{q}^{k}_{h}(t,\tau)=\mathrm{e}^{\mathrm{i}k\tilde{\omega}t}\tilde{\xi}^{k}(t,\tau).$ Multiplication of this result with $P$ yields

[TABLE]

where $x_{h}(t)=\sum\limits_{|k|<N}x^{k}_{h}(t)$ with $x^{k}_{h}(t)=\mathrm{e}^{\mathrm{i}k\tilde{\omega}t}\zeta^{k}(t)$ and $q_{h}(t,\tau)=\sum\limits_{|k|<N}q_{h}^{k}(t,\tau)$ with $q_{h}^{k}(t,\tau)=\mathrm{e}^{\mathrm{i}k\tilde{\omega}t}\xi^{k}(t,\tau).$ Rewrite the equation in terms of $x^{k}_{h}$ and then we get

[TABLE]

where $\mathcal{U}(\vec{q}(t,\tau))$ is defined as

[TABLE]

with

[TABLE]

Define the vector function

[TABLE]

We have the invariance property that $\mathcal{U}(\vec{q}(\lambda,t,\tau))$ is independent of $\lambda$ and $\tau$. Thus, differentiation with respect $\lambda$ implies

[TABLE]

By letting $\lambda=0$ and $\tau=\frac{1}{2}$, we obtain $\sum\limits_{|k|<N}\mathrm{i}k\tilde{\omega}(q^{k}_{h}(t,\frac{1}{2}))^{\intercal}\nabla_{k}\mathcal{U}(\vec{q}(t,\frac{1}{2}))=0.$ Consequently, we have

[TABLE]

According to the relation of the coefficient functions given in Theorem 5.3, we obtain

[TABLE]

By Proposition 5.1 and the following the “magic formulas” on p. 508 of [16], it can be verified that the right-hand side of (33) is a total derivative. This proves that there exists a function $\widehat{\mathcal{M}}$ such that $\frac{d}{dt}\widehat{\mathcal{M}}[\vec{\tilde{\zeta}},\vec{\tilde{\eta}}]=\mathcal{O}(h^{N})$. The first statement follows by integration this result.

Based on the bounds of the coefficients functions given in Theorem 5.2, the $\widehat{\mathcal{H}}$ can be expressed as

[TABLE]

where the second formula of (26) was used. This implies the second statement of the theorem.

5.3 Long time magnetic moment conservation

From the bounds presented in Theorem 5.2 and from the expression (14) it follows that

[TABLE]

Looking closing at this formula and (30), we get the following relationship between the magnetic moment and the almost-invariant $\widehat{\mathcal{M}}$:

[TABLE]

Based on the above analysis and following the way used in Chapter XIII of [16], Theorem 3.3 is easily proved by patching together the local near-conservation result.

6 Conclusions

In this paper, we have studied exponential energy-preserving methods for charged-particle dynamics in a strong and constant magnetic field. It was shown that this method can exactly preserve the energy of the charged-particle dynamics. It is worth mentioning that the long-time magnetic moment conservation was also been researched by deriving a modulated Fourier expansion of the method and showing an almost-invariant of the modulation system. In this paper, we have also studied other properties of the method. The effectiveness of the method is emphasized by carrying out a numerical experiment.

Acknowledgements

The author is grateful to Professor Christian Lubich for drawing my attention to the charged-particle dynamics and the long-term analysis of energy-preserving methods.