3.1. Composing a first-order scheme with its adjoint?

A common method for constructing a second-order-accurate integrator starting from a first-order-accurate integrator is to compose it with its adjoint (see Hairer et al. Reference Hairer, Lubich and Wanner2006), discussed above for the special case of the symplectic Euler scheme. To show what can go wrong with composing a first-order scheme with its adjoint in the context of variational integration, let us introduce the discrete phase-space Lagrangian for a canonical system in one degree of freedom:

(3.2)\begin{equation} L_{d}=p_{k}\frac{q_{k+1}-q_{k}}{h}-H(q_{k},p_{k}). \end{equation}

Except for the change $H(q_{k+1},p_{k})\to H(q_{k},p_{k})$, this discretization is identical to that of the symplectic Euler scheme of (2.20) and (2.21a,b). (The reason for this modification of the symplectic Euler Lagrangian is that it has an $O(h)$ term in its preserved two-form, similar to the term in the magnetic field line and guiding centre examples. Ramifications of this $O(h)$ term are discussed below.) Variations with respect to $p_{k}$ and $q_{k}$ yield the map

(3.3a,b)\begin{equation} q_{k+1}=q_{k}+hH_{2}(q_{k},p_{k}),\quad p_{k+1}=p_{k}-hH_{1}(q_{k+1},p_{k+1}). \end{equation}

The first of these is explicit with respect to $q$; the second update is implicit in $p$ and leapfrogged in $q$, so that this scheme is slightly different from the updating in symplectic Euler.Footnote ² The adjoint of this scheme follows from $k\leftrightarrow k+1, h\rightarrow -h$. We find

(3.4)\begin{equation} L_{d}=p_{k+1}\frac{(q_{k+1}-q_{k})}{h}-H(q_{k+1},p_{k+1}), \end{equation}

leading to

(3.5a,b)\begin{equation} p_{k+1}=p_{k}-hH_{1}(q_{k},p_{k}),\quad q_{k+1}=q_{k}+hH_{2}(q_{k+1},p_{k+1}). \end{equation}

This scheme is explicit in $p$, implicit and leapfrogged in $q$.

By direct substitution, we find that

(3.6)\begin{equation} \left(1+hH_{12}(q,p)\right)\,{\rm d}q\wedge {\rm d}p\end{equation}

is preserved by this scheme. As discussed in Ellison et al. (Reference Ellison, Finn, Burby, Kraus, Qin and Tang2018), the preservation of this two-form can also be shown by inspecting the first and last terms in ${\rm d}S={\rm d}\sum _{k}hL_{d}$. This is the discrete analogue to the variation of $S=\int _{0}^{T}(p\dot {q}-H)\,{\rm d}t$, and $\delta S=0$ gives $\int _{0}^{T}[(\dot {q}-H_{p})\,{\rm d}p+(-\dot {p}-H_{q})\,{\rm d}q]\,{\rm d}t+p\,{\rm d}q\vert _{0}^{T}$. Here, the vanishing of the integral provides Hamiltonian equations, and the endpoints, with $d^{2}S=0$, imply the conservation of the canonical two-form $\omega ={\rm d}q\wedge {\rm d}p$, showing the canonical symplectic property of the equations. Similarly, for a general non-canonical system, the endpoint terms in $S=\int (\vartheta _i(z)\dot {z}^i-H(z))\,{\rm d}t$ show the preservation of $\omega _{ij}\,{\rm d}z^i \wedge {\rm d}z^j$, where $\omega _{ij}=\vartheta _{i,j}-\vartheta _{j,i}$. The preservation of this form shows the non-canonical symplectic nature of the discrete map. This property is consistent with the fact that this scheme and the adjoint symplectic Euler scheme are equivalent under a non-canonical change of variables. Unlike the symplectic Euler scheme, which preserves the canonical two-form $\omega ={\rm d}q\wedge {\rm d}p$, this form has $\omega =\omega (h)=(1+O(h))\,{\rm d}q\wedge {\rm d}p$. That is, in this scheme the phase-space coordinates $(q,p)$ are not canonical. By analogous arguments, or by direct inspection, the adjoint scheme preserves

(3.7)\begin{equation} \omega^{{{\dagger}}}(h)=\omega({-}h)=\left(1-hH_{12}(q,p)\right)\,{\rm d}q\wedge {\rm d}p, \end{equation}

and because the correction is $O(h)$, the two-forms differ by $O(h)$. From these observations, when $\omega (h)\neq \omega (-h)$ it is not clear how to determine which two-form, if any, is preserved by the composition of the scheme and its adjoint. In § 5.1 we show numerical evidence that such a composed form does not in general preserve any two-form. Also, in Appendix C we show the direct analogy to the symplectic Euler scheme in (2.21a,b) (not (3.3a,b)), and show that its preserved two-form also has $O(h)$ corrections. While we cannot rule out the existence of a first-order-accurate scheme in our class of non-canonical variables (which would allow composition with its adjoint to obtain second-order accuracy), we proceed in the next section to show two separate second-order-accurate schemes for such systems.

3.2. Centred schemes for second-order accuracy

Because composing a first-order DVI with its adjoint does not reliably produce a scheme preserving a symplectic form (although such a scheme is second-order accurate), we are naturally led to consider the problem of proceeding to higher order by identifying improved properly degenerate discrete Lagrangians. We first illustrate for a one-degree-of-freedom case ($m=2$) in canonical variables, introducing two different staggered, centred schemes to discretize the action $S=\int L\,{\rm d}t$ for the phase-space Lagrangian $L(q,p,\dot {q},\dot {p})=p\dot {q}-H(q ,p)$. The first scheme has

(3.8)\begin{equation} L_{d}(q_{k},q_{k+1},p_{k+1/2})=p_{k+1/2}\frac{(q_{k+1}-q_{k})}{h}-H\left(\frac{q _{k}+q_{k+1}}{2},p_{k+1/2}\right). \end{equation}

This scheme includes a staggered nature of $p$ and a midpoint nature with respect to $q$. Taking variations with respect to $q_{k}$ and $p_{k+1/2}$ we find

(3.9)\begin{gather} q_{k+1} =q_k+h H_2\left(\frac{q_k+q_{k+1}}{2},p_{k+1/2}\right), \end{gather}

(3.10)\begin{gather} p_{k+1/2} =p_{k-1/2}-\frac{h}{2} H_1\left(\frac{q_k+q_{k+1}}{2},p_{k+1/2}\right) -\frac{h}{2} H_1\left(\frac{q_{k-1}+q_{k}}{2},p_{k-1/2}\right). \end{gather}

We call this the midpoint DVI (MDVI) scheme. This scheme is clearly time-centred, leading to second-order accuracy, as we discuss below. This scheme must be advanced implicitly in both variables. It can easily be shown to be properly degenerate, by showing that the rank of the $2\times 2$ discrete Hessian is one, essentially because, as for the symplectic Euler scheme, its adjoint and the schemes of (3.2) and (3.4), the discrete Lagrangian $L_d$ depends on $p$ at only one time level. This degeneracy is in spite of the fact that the scheme appears to be a two-step scheme, connecting $q_{k-1}$, $q_{k}$ and $q_{k+1}$ (but only $p_{k-1/2}$ and $p_{k+1/2}$). As in Ellison et al. (Reference Ellison, Finn, Burby, Kraus, Qin and Tang2018), $q_{k-1}$ can be expressed in terms of $q_{k}$ and $p_{k-1/2}$ by the implicit equation $q_{k-1} =q_{k}-h H_2(({q_{k-1}+q_{k}})/{2},p_{k-1/2})$ obtained from (3.9), retarded by one step. With this substitution, we define a time-advance map from (3.9) and (3.10) at the current step of the form $(q_k, p_{k-1/2}) \mapsto (q_{k+1}, p_{k+1/2})$, obtaining a one-step method. This property of appearing to be of two steps but showing a single-step nature after some substitutions, as discussed at length in Ellison et al. (Reference Ellison, Finn, Burby, Kraus, Qin and Tang2018), shows the importance of the discrete Hessian test; indeed, without the assurance of the discrete Hessian test, it would be easy to miss the possibility of this substitution.

The second scheme uses

(3.11)\begin{equation} L_{d}(q_{k},q_{k+1},p_{k+1/2})=p_{k+1/2}\frac{(q_{k+1}-q_{k})}{h}- \frac{1}{2}H\left(q_{k},p_{k+1/2}\right)-\frac{1}{2}H\left(q_{k+1},p_{k+1/2} \right), \end{equation}

and variations lead to

(3.12)\begin{gather} q_{k+1}=q_k+\frac{h}{2}H_2(q_k,p_{k+1/2})+\frac{h}{2}H_2(q_{k+1},p_{k+1/2}), \end{gather}

(3.13)\begin{gather} p_{k + 1/2}=p_{k - 1/2}-\frac{h}{2} H_1\left(q_{k},p_{k + 1/2}\right) -\frac{h}{2} H_1\left(q_{k},p_{k-1/2}\right). \end{gather}

This second scheme is also time-centred, again showing second-order accuracy. This scheme must be advanced implicitly in both variables. For this scheme, the discrete Lagrangian is obtained using a trapezoidal quadrature scheme, hence we call this scheme the trapezoidal DVI (TDVI) scheme. In this scheme, the single-step nature is evident from the DEL equations; of course, the Hessian test confirms the single-step character. Finally, backward error analysis for both the MDVI and the TDVI schemes, with $q_{k}=q(t=kh)$ and $p_{k+1/2}=p((k+1/2)h)$, shows second-order accuracy. This derivation follows from a Taylor expansion about $t_k=kh$, e.g. $q_{k+1}=q(t_k)+h\dot {q}(t_k)+(h^2/2)\ddot {q}(t_k)+O(h^3),$ $p_{k+1/2}=p(t_k)+(h/2)\dot {p}(t_k)+(h^2/8)\ddot {p}(t_k)+O(h^3)$; upon substituting in either (3.9) and (3.10) or (3.12) and (3.13) and dividing by $h$, we find that the $O(h)$ terms all cancel, proving second-order accuracy.

In both schemes, if initial conditions $q_{0}$ and $p_{0}$ are both given at $t=0$, the momentum variable needs to be regressed to $p_{-1/2}$ by a processing scheme, which we discuss shortly. Numerical trials using the non-reversible (cf. Appendix B) Hamiltonian $H=(p^{2}+q^{2})/2+\alpha q p^{3}/3$ indicate second-order accuracy and the good long-time behaviour of a scheme with a preserved two-form.

For the non-canonical (but not completely general) phase-space Lagrangians of the form of (3.1), the MDVI scheme has

(3.14)\begin{equation} L_d(x_k,y_{k + 1/2}, x_{k+1}) = f_i\left( \frac{x_k+x_{k+1}}{2}, y_{k+1/2}\right ) \frac{x^i_{k+1}-x^i_k}{h} - H\left ( \frac{x_k+x_{k+1}}{2},y_{k+1/2}\right). \end{equation}

Again, the centredness suggests, and backward error analysis indeed shows, second-order accuracy. Also, the staggered-grid discrete Lagrangian (3.14) preserves a symplectic form determined, as before, by examining the endpoint terms in the variation of the action. Again, the discrete Hessian has half-rank because $L_d$ depends on $y$ at only one time level.

The DEL equations stemming from (3.14) are given by

(3.15a)\begin{align} & \tfrac{1}{2} \left( f_{i,j}\left(k + \tfrac{1}{2}\right)(x_{k+1}^i - x_{k}^i) + f_{i,j}\left(k-\tfrac{1}{2}\right) (x_k^i - x_{k-1}^i) \right) \nonumber\\ & \quad - \left(f_j\left(k + \tfrac{1}{2}\right) - f_j\left(k - \tfrac{1}{2}\right)\right) - \frac{h}{2} \left(H_{,j}\left(k + \tfrac{1}{2}\right) + H_{,j}\left(k-\tfrac{1}{2} \right) \right) = 0, \end{align}

(3.15b)\begin{align} & f_{i,\alpha}\left(k + \tfrac{1}{2} \right) \left(x_{k+1}^i - x_k^i\right) - h H_{,\alpha}\left(k + \tfrac{1}{2}\right) = 0, \end{align}

where $i,j \in \{1,\ldots, m/2\}$, $m$ is even, while $\alpha \in \{m/2+1,\ldots, m\}$ and $(k+1/2)$ is shorthand for the arguments $(({x_{k+1} + x_k})/{2}, y_{k+1/2} )$ and $(k-1/2)$ is analogous. Like the canonical MDVI, the apparent dependence on $x_{k-1}$ may be eliminated in favour of a function depending on $(x_k, y_{k-1/2})$. In practice one has access to $x_{k-1}$ due to the state at prior iterations (except for on the first step) and the stored value may be used instead of directly solving for $x_{k-1}$ at each new iteration.

The TDVI scheme for this class of non-canonical cases arises from

(3.16)\begin{gather} L_d(x_k,y_{k+1/2},x_{k+1}) = \tfrac{1}{2}\left(f_i(x_k, y_{k+1/2})+f_i(x_{k+1}, y_{k+1/2})\right) \frac{x^i_{k+1}-x^i_k}{h}\nonumber\\ - \tfrac{1}{2}\left(H(x_{k},y_{k+1/2}) + H(x_{k+1},y_{k+1/2}) \right). \end{gather}

Again, the time-centred property leads to second-order accuracy, and its discrete Hessian shows that it is properly degenerate, again because $L_d$ depends on $y$ at only one time level.

Performing variations with respect to $x_{k}, y_{k+1/2}$ yields the TDVI scheme:

(3.17a)\begin{align} & \tfrac{1}{2} \left(f_{i,j}(x_k, y_{k+1/2}) \left( x_{k+1}^i - x_{k}^i \right) + f_{i,j}(x_k, y_{k-1/2}) \left( x_{k}^i - x_{k-1}^i \right) \right) \nonumber\\ & \quad-\tfrac{1}{2} \left(f_j(x_{k+1}, y_{k+1/2}) + f_{j}(x_{k}, y_{k+1/2}) - f_j(x_k, y_{k-1/2}) - f_j(x_{k-1}, y_{k-1/2} \right) \nonumber\\ & \quad- \tfrac{1}{2} \left( H_{,j}(x_k, y_{k+1/2}) + H_{,j}(x_{k+1}, y_{k+1/2}) \right) = 0, \end{align}

(3.17b)\begin{align} & \tfrac{1}{2} \left( f_{i,\alpha} (x_k, y_{k+1/2}) + f_{i,\alpha}(x_{k+1}, y_{k+1/2}) \right) \left(x_{k+1}^i - x_{k}^i\right) \nonumber\\ & \quad- \tfrac{h}{2} \left( H_{,\alpha}(x_k, y_{k+1/2}) + H_{,\alpha}(x_{k+1}, y_{k+1/2}) \right) = 0. \end{align}

In contrast to the canonical setting, the TDVI scheme introduces dependence on $x_{k-1}$ in general. Of course, this dependence is superficial, and can be eliminated in favour of a function of $(x_k, y_{k-1/2})$ as previously discussed. Relative to the MDVI scheme, the TDVI scheme requires more function evaluations (i.e. additional evaluations of $f$ and $H$) in the update rule; this has the potential to increase the computational expense of TDVI relative to MDVI.

Of course, initial conditions for the time-marching scheme defined by (3.14) or (3.16) will be supplied at an integer time step instead of directly on the staggered grid. In order to transform integer-time-step initial conditions $(x_0,y_0)$ to staggered-grid initial conditions $(x_0,y_{-1/2})$, it is sufficient to advance $y$ backward in time by a half-step $h/2$ using any scheme that is accurate to first (or higher) order. This transformation is encoded in the mapping $\varphi _{-h/2}:(x_0,y_0)\mapsto (x_0,y_{-1/2})$. Similarly, at the end of a simulation, the staggered-grid data $(x_N,y_{N-1/2})$ must be collocated to the integer-grid data $(x_N,y_N)$. The natural way to do this is simply to apply the inverse of $\varphi _{-h/2}$, i.e. set $(x_N,y_N) = \varphi _{-h/2}^{-1}(x_N,y_{N-1/2})$, and second-order accuracy is preserved. These two conditions are special cases of enforcing the relationship

(3.18)\begin{equation} \varphi_{{-}h/2}(x_k,y_{k}) = (x_k,y_{k-1/2}) \end{equation}

for all $k$. Indeed, for displaying results during a computation, e.g. where $k$ equals a multiple of a fundamental period $M$, typically with $M>1$, these results will retain second-order accuracy if the points are collocated according to (3.18). Furthermore, this process preserves the symplectic nature of the scheme. This is because applying the processing scheme at the initial step and its inverse at the final step essentially defines the step in terms of a non-symplectic change of variables. Writing the processing step in (3.18) as $\varPhi$ and the time step as $T_h$, this new scheme is written as $\tilde {T}_h=\varPhi ^{-1}\circ T_h \circ \varPhi$. For a symplectic scheme in canonical variables, this defines an equivalent map in non-canonical variables; for the case under consideration here, where the symplectic scheme is already in non-canonical variables, this transformation just specifies an $O(h)$ transformation to other non-canonical variables. Regarding invariants, the original symplectic scheme may not preserve invariants exactly, but typically such a scheme will preserve invariants such as energy approximately for sufficiently small $h$, preventing, for example, secular or exponential growth in the energy. The same properties hold for the modified scheme $\tilde {T}_h$. In Blanes et al. (Reference Blanes, Casas and Murua2004), the idea of increasing the order of a low-order scheme using a map and its inverse as pre- and post-processors is explored in greater detail. It would be interesting to determine if even higher-order DVIs may be derived using more elaborate processing than the more-or-less obvious processors described here.

4.1. Magnetic field line

For the problem of tracing magnetic field lines, we take the action to be equal to the flux:

(4.1)\begin{equation} \varPhi= \int L\left(\boldsymbol{x}, \frac{{\rm d}\boldsymbol{x}}{{\rm d}t}\right) \,{\rm d}t = \int \boldsymbol{A}(\boldsymbol{x}) \boldsymbol{\cdot} \frac {{\rm d}\boldsymbol{x}}{{\rm d}t}\,{\rm d}t,\end{equation}

where $\boldsymbol {A}$ is the magnetic vector potential. The invariance with respect to reparameterization of time (see Ellison et al. Reference Ellison, Finn, Burby, Kraus, Qin and Tang2018) in (4.1) is consistent with its Euler–Lagrange equation, namely $({\rm d}\boldsymbol {x}/{\rm d}t)\times \boldsymbol {B}=0$, where $\boldsymbol {B}=\boldsymbol {\nabla } \times \boldsymbol {A}$ is the magnetic field: this equation determines the direction of the flow but not the speed. This time invariance is dealt with by parameterizing the field line trajectory in terms of one of the coordinates (without loss of generality, we choose $x^3$) instead of ‘time’ $t$:

(4.2)\begin{align} \varPhi=\int L(x^1,x^2,x^3,{{\rm d} x}^2/{{\rm d} x}^3,x^3)\,{{\rm d} x}^3 =\int \left( A_2(x^1,x^2,x^3)\,{{\rm d} x}^2/{{\rm d} x}^3+A_3(x^1,x^2,x^3)\right)\,{{\rm d} x}^3.\end{align}

This Lagrangian is of the form in (3.1), with $f = [0 \quad A_2]$ and $H=-A_3$; the Euler–Lagrange equations for (4.2) are

(4.3a)\begin{gather} {{\rm d} x}^1/{{\rm d} x}^3 =(A_{3,2}-A_{2,3})/A_{2,1} = B^1/B^3, \end{gather}

(4.3b)\begin{gather}{{\rm d} x}^2/{{\rm d} x}^3 ={-}A_{3,1}/A_{2,1} = B^2/B^3. \end{gather}

The MDVI then follows from (3.15), which when expressed in terms of the magnetic vector potential becomes

(4.4a)\begin{align} & A_{2,1}(k + \tfrac{1}{2}) \left(x^2_{k+1} - x^2_k \right) + h A_{3,1}(k+\tfrac{1}{2}) = 0, \end{align}

(4.4b)\begin{align} & \tfrac{1}{2}\left(A_{2,2}(k + \tfrac{1}{2}) \left(x^2_{k+1} - x^2_k \right) + A_{2,2}(k - \tfrac{1}{2}) \left(x^2_{k} - x^2_{k-1} \right) \right) \nonumber\\ & \quad -\left( A_2(k+\tfrac{1}{2}) - A_2(k-\tfrac{1}{2}) \right) + \tfrac{h}{2}\left( A_{3,2}(k+\tfrac{1}{2}) + A_{3,2}(k-\tfrac{1}{2}) \right) = 0, \end{align}

where $(k+\tfrac {1}{2})$ denotes evaluation at $(x^1_{k+1/2}, (x^2_{k+1} + x_k^2)/2, x^3_{k+1/2})$. Similar to the process in Ellison (Reference Ellison2016) and Ellison et al. (Reference Ellison, Finn, Burby, Kraus, Qin and Tang2018), these equations are made into an explicitly single-step process by writing the first with $k\to k-1$, solving for $x_{k-1}$ and substituting for $x_{k-1}$ in the remaining equations without the $k\to k-1$ incrementation.

Similarly, the TDVI algorithm for the magnetic field line problem follows from (3.17):

(4.5a)\begin{align} & \tfrac{1}{2} \left(A_{2,1}(x_{k+1/2}^1, x_{k+1}^2, x_{k+1/2}^3) + A_{2,1}(x_{k+1/2}^1, x_{k}^2, x_{k+1/2}^3) \right) \left(x_{k+1}^2 - x_k^2 \right) \nonumber\\ & \quad+ \tfrac{h}{2} \left(A_{3,1}(x_{k+1/2}^1, x_{k}^2, x_{k+1/2}^3) + A_{3,1}(x_{k+1/2}^1, x_{k+1}^2, x_{k+1/2}^3) \right) = 0,\end{align}

(4.5b)\begin{align} & \tfrac{1}{2} \left(A_{2,2}(x^1_{k+1/2}, x^2_k, x^3_{k+1/2})(x^2_{k+1} - x^2_{k}) + A_{2,2}(x^1_{k-1/2}, x^2_k, x^3_{k-1/2})(x^2_k - x^2_{k-1}) \right) \nonumber\\ & \quad+ \tfrac{1}{2} \left( A_2(x^1_{k-1/2}, x^2_{k-1}, x^3_{k-1/2}) + A_2(x^1_{k-1/2}, x^2_{k}, x^3_{k-1/2})\right) \nonumber\\ & \quad- \tfrac{1}{2} \left( A_2(x^1_{k+1/2}, x^2_{k}, x^3_{k+1/2}) + A_2(x^1_{k+1/2}, x^2_{k+1}, x^3_{k +1/2})\right) \nonumber\\ & \quad+ \frac{h}{2} \left( A_{3,2}(x^1_{k-1/2}, x^2_{k}, x^3_{k-1/2}) + A_{3,2}(x^1_{k+1/2}, x^2_{k}, x^3_{k+1/2}) \right) = 0. \end{align}

Numerical tests of the MDVI and TDVI schemes are presented in § 5.

4.2. Guiding centre equations

We treat the example of the guiding centre equations based on discretizing the toroidally regularized (compare with Burby & Ellison Reference Burby and Ellison2017) phase-space Lagrangian

(4.6)\begin{align} L(z,\dot{z})& = A_2(\boldsymbol{x})\dot{x}^2+(A_3(\boldsymbol{x})+u)\dot{x}^3 - H_{\text{gc}} \nonumber\\ & \equiv A^{\dagger}_2(\boldsymbol{x})\,\dot{x}^2 + A^{\dagger}_3(\boldsymbol{x},u)\,\dot{x}^3 - H_{\text{gc}}, \end{align}

where $H_{{\rm gc}}=\tfrac {1}{2}(B^3/B)^2\,u^2+\mu B - \tfrac {1}{2}E_\perp ^2/B^2 - u\boldsymbol {b}\boldsymbol {\cdot } (\boldsymbol {E}\times \boldsymbol {\nabla } x^3)/B+\varphi$ is the guiding centre Hamiltonian, in toroidally regularized non-canonical variables, $u$ is the toroidally regularized parallel velocity, $\mu$ is the magnetic moment, $B$ is the magnitude of the magnetic field and $\varphi$ is the scalar potential for the electric field $\boldsymbol {E} = -\boldsymbol {\nabla } \varphi$. For simplicity, we neglect time dependence in the electromagnetic fields. Toroidal regularization requires $|B^3|>0$ everywhere. In toroidal geometries relevant to magnetic fusion energy, possible choices for coordinates include $(x^1,x^2,x^3) = (R,Z,\phi )$ or $(x^1,x^2,x^3) = (r,\theta,\phi )$, where $(R,Z,\phi )$ denote cylindrical coordinates and $(r,\theta,\phi )$ denote toroidal coordinates. As in the magnetic field example, the gauge condition $A_1=0$ was imposed to lead to proper degeneracy. (The original form introduced in Ellison et al. (Reference Ellison, Finn, Burby, Kraus, Qin and Tang2018) and Ellison (Reference Ellison2016) required that the covariant component $B_1$ of the magnetic field also vanish, but Burby & Ellison (Reference Burby and Ellison2017) subsequently relaxed that condition.)

The second-order-accurate DVIs follow by establishing the correspondence with (3.1); in this case, $f = [0 \quad A_2^{\dagger} \quad A_3^{\dagger} \quad 0]$ and $H = H_{{\rm gc}}$. We discretize this phase-space Lagrangian by the MDVI and TDVI schemes for second-order accuracy. For the MDVI discretization of the guiding centre system, we use the discrete Lagrangian

(4.7)\begin{align} & L_d(r_{k+1/2}, \theta_{k}, \theta_{k+1}, \phi_k, \phi_{k+1}, u_{k+1/2}) =A_\theta(k+\tfrac{1}{2}) (\theta_{k+1} - \theta_k)\nonumber\\ & \quad + \left( A_\phi(k+\tfrac{1}{2}) + u_{k+1/2} \right) (\phi_{k+1} - \phi_k) - h H(k+\tfrac{1}{2}). \end{align}

The MDVI scheme, following from (3.15), is

(4.8a)\begin{align} & A_{\theta, r}(k+\tfrac{1}{2}) (\theta_{k+1} - \theta_k) + A_{\phi, r}(k+\tfrac{1}{2}) (\phi_{k+1} - \phi_k) - h H_{,r}(k+\tfrac{1}{2}) = 0, \end{align}

(4.8b)\begin{align} & \tfrac{1}{2} A_{\theta,\theta}(k-\tfrac{1}{2}) (\theta_k - \theta_{k-1}) + \tfrac{1}{2} A_{\theta, \theta}(k+\tfrac{1}{2}) (\theta_{k+1} - \theta_k) \nonumber\\ & \quad+ \tfrac{1}{2} A_{\phi,\theta}(k-\tfrac{1}{2}) (\phi_k - \phi_{k-1}) + \tfrac{1}{2} A_{\phi, \theta} (k+\tfrac{1}{2}) (\phi_{k+1} - \phi_k) \nonumber\\ & \quad+ A_\theta(k-\tfrac{1}{2}) + u_{k-1/2} - A_\theta(k+\tfrac{1}{2}) - u_{k+1/2} - \tfrac{h}{2} \left( H_{,\theta}(k-\tfrac{1}{2}) + H_{,\theta}(k+\tfrac{1}{2}) \right) = 0, \end{align}

(4.8c)\begin{align} & \tfrac{1}{2} A_{\theta,\phi}(k-\tfrac{1}{2}) (\theta_k - \theta_{k-1}) + \tfrac{1}{2} A_{\theta, \phi}(k+\tfrac{1}{2}) (\theta_{k+1} - \theta_k) \nonumber\\ & \quad+ \tfrac{1}{2} A_{\phi,\phi}(k-\tfrac{1}{2}) (\phi_k - \phi_{k-1}) + \tfrac{1}{2} A_{\phi, \phi} (k+\tfrac{1}{2}) (\phi_{k+1} - \phi_k) \nonumber\\ & \quad+ A_\phi(k-\tfrac{1}{2}) + u_{k-1/2} - A_\phi(k+\tfrac{1}{2}) - u_{k+1/2} - \tfrac{h}{2} \left( H_{,\phi}(k-\tfrac{1}{2}) + H_{,\phi}(k+\tfrac{1}{2}) \right) = 0, \end{align}

(4.8d)\begin{align} \phi_{k+1} - \phi_k - \tfrac{h}{2} H_{,u}(k+\tfrac{1}{2}) = 0, \end{align}

where $(k+1/2)$ refers to $(r_{k+1/2},(\theta _{k}+\theta _{k+1})/2,(\phi _{k}+\phi _{k+1})/2)$ and similarly for $(k-1/2)$. The incrementation–substitution process for obtaining a scheme that is explicitly single step, introduced in Ellison (Reference Ellison2016) and Ellison et al. (Reference Ellison, Finn, Burby, Kraus, Qin and Tang2018) and used in the magnetic field lines case above, is also used here.

For the TDVI scheme, we use the discrete Lagrangian in (3.16), namely

(4.9)\begin{align} & L_d(r_{k+1/2}, \theta_{k}, \theta_{k+1}, \phi_k, \phi_{k+1}, u_{k+1/2}) \nonumber\\ & \quad= \tfrac{1}{2} \left( A_\theta(r_{k+1/2}, \theta_k, \phi_k) + A_\theta(r_{k+1/2}, \theta_{k+1}, \phi_{k+1}) \right) (\theta_{k+1} - \theta_k) \nonumber\\ & \qquad+ \tfrac{1}{2} \left( A_\phi(r_{k+1/2}, \theta_k , \phi_k ) + A_\phi(r_{k+1/2}, \theta_{k+1} , \phi_{k+1} ) + 2 u_{k+1/2} \right) (\phi_{k+1} - \phi_k) \nonumber\\ & \qquad- \tfrac{h}{2} \left( H(r_{k+1/2}, \theta_k , \phi_k, u_{k+1/2} ) + H(r_{k+1/2}, \theta_{k+1} , \phi_{k+1}, u_{k+1/2})\right]. \end{align}

The guiding centre TDVI update equations follow from straightforward variations of this discrete Lagrangian (or from (3.17)) and are omitted here for brevity.

5.1. Failure of composing with the adjoint

In § 3.1 we noted an example of a scheme with a preserved two-form with $\omega (h)\neq \omega (-h)$, suggesting complications if the scheme and its adjoint are composed in an attempt to obtain second-order accuracy. In this section we give a concrete example for which the composed scheme appears not to have a preserved two-form. We consider the non-reversible Hamiltonian

(5.1)\begin{equation} H=\frac{p^{2}+q^{2}}{2}+\frac{\alpha q p^{3}}{3}. \end{equation}

Numerical tests show that both of the schemes of (3.3a,b) and (3.5a,b) have good long-time behaviour for this Hamiltonian. The orbits of the two first-order schemes have bounded behaviour but, with first-order accuracy, have noticeable oscillations in the value of $H$ between bounds. The composed scheme, which has smaller oscillations in $H$, otherwise performs poorly: the points that should stay near the constant $H$ surfaces spiral out, as shown in figure 1. The behaviour of $H$ in time shows exponential increase, with $\gamma =O(h^2)$. Appendix B shows an example for which composition does appear to give useful results, but this particular example is a reversible system, and this reversibility by itself appears to be responsible for the positive results.

Figure 1.

Phase portrait (a) and energy versus time (b) for two different methods of integrating the non-reversible Hamiltonian $H=(p^{2}+q^{2})/2+\alpha q p^{3}/3$. The ‘direct’ scheme of (3.3a,b), and shown in green, is first-order-accurate and preserves a discrete symplectic structure that depends on the numerical time-step size $h$, with similar results for the adjoint scheme in (3.5a,b). Composing the direct method with its adjoint achieves second-order accuracy by virtue of centring in time, giving smaller oscillations in $H$. However, the numerical results indicate that the composed scheme leads to growth in $H$ and therefore does not preserve a two-form.

5.2. Magnetic field line

To test the proposed algorithms in a magnetic configuration representative of those of interest to the magnetic fusion community, we use the simple analytic expression for an axisymmetric, toroidal magnetic field presented in Qin, Guan & Tang (Reference Qin, Guan and Tang2009):

(5.2)\begin{equation} \boldsymbol{ A}(r, \theta, \phi) = \frac{B_0 R_0}{\cos^2 \theta} \left(r \cos \theta -R_0 \log \left( 1 + \frac{r \cos \theta}{R_0}\right) \right)\boldsymbol{\nabla} \theta - \frac{B_0 r^2}{2 q_0} \boldsymbol{\nabla} \phi, \end{equation}

where $B_0$ is a magnetic field amplitude, $R_0$ is the major radius and $q_0=\sqrt {2}$ is the on-axis safety factor. The variables $(x^1,x^2,x^3)$ of (4.2) are replaced by simple toroidal coordinates $(r,\theta,\phi )$, and $\phi$ takes the place of the time variable, as discussed in § 4.1.

First, we demonstrate numerically that the MDVI and TDVI achieve the anticipated second-order accuracy. Next, we demonstrate that the proposed algorithms exhibit the expected qualitative behaviour of symplectic integrators. To the axisymmetric magnetic field of (5.2) we add a perturbation of the form

(5.3)\begin{equation} \boldsymbol{ A}(r, \theta, \phi) = \boldsymbol{ A}_0(r, \theta, \phi) - \frac{B_0 r^2}{2 q_0}\sum_i \delta_i \sin(m_i \theta - n_i \phi)\boldsymbol{\nabla} \phi, \end{equation}

where $\boldsymbol { A_0}$ is given by (5.2). We choose two perturbative harmonics, $m_0=3, n_0=2$ and $m_1=7, n_1=5$, with amplitudes $\delta _0 = \delta _1 = 10^{-4}$. These perturbations lead to magnetic islands at the resonant magnetic surfaces, and small stochastic field line regions in the $(m_1,n_1)=(3,2)$ and $(m_2,n_2)=(7,5)$ resonant regions.

We test the order of accuracy by varying the step $\Delta \phi$ in factors of two across a range of $2^5$ and comparing with a fourth-order Runge–Kutta scheme with an extremely small value of $\Delta \phi$. See figure 2. We compare the MDVI and TDVI schemes with each other and with the first-order variant described in Ellison et al. (Reference Ellison, Finn, Burby, Kraus, Qin and Tang2018) in which the three variables $(r,\theta,\phi )$ are collocated, i.e. not staggered. We integrate over a large number of toroidal transits for this comparison. Second-order accuracy is confirmed for the MDVI and TDVI schemes, as is first-order accuracy for the non-staggered scheme. The MDVI and TDVI schemes exhibit relatively similar accuracy, with the MDVI scheme being more accurate for this particular example.

Figure 2.

Comparison of the order of accuracy of the MDVI scheme, the TDVI scheme and the collocated first-order-accuratescheme of Ellison et al. (Reference Ellison, Finn, Burby, Kraus, Qin and Tang2018) for the magnetic field line problem, over a range of $2^5$ in $\Delta \phi$. The dashed curves show Error$\sim h$ (blue) and Error$\sim h^2$ (green) for comparison.

The Poincaré surface of section at $\phi =0$ is shown in figure 3. Figures 3(a) and 3(b) show the second-order Runge–Kutta (RK2) scheme, with $\Delta \phi =0.1$ and $\Delta \phi =5\times 10^{-4}$, respectively. Figures 3(c) and 3(d) show the MDVI scheme for the same two values of $\Delta \phi$. The RK2 results in figure 3(a) show blurriness of the Kolmogorov–Arnold–Moser (KAM) surfaces, falsely indicating a higher degree of magnetic stochasticity. The results in figure 3(b) are greatly improved. The two MDVI cases in figure 3(c,d), for very different steps $\Delta \phi$, look almost identical, showing very good preservation of KAM tori. The TDVI scheme leads to results essentially indistinguishable from those of the MDVI scheme for comparable time steps.

Figure 3.

Poincaré section $\phi =0$ for the magnetic field line problem. For comparison we show results from (a) the RK2 scheme with $\Delta \phi = 0.1$, (b) the RK2 scheme with $\Delta \phi = 0.005$, (c) the MDVI scheme with $\Delta \phi = 0.1$ and (d) the MDVI scheme with $\Delta \phi = 0.005$. In all cases we have $\phi _\text {final} = 3\times 10^5$. The $(m,n)=(3,2)$ and $(7,5)$ island chains are evident, and in the RK2 scheme in (a) the KAM surfaces are blurred.

5.3. Guiding centre

In this section we show numerical results for the guiding centre example, with time-independent potentials and scalar potential $\varphi =0$. That is, physically, the electric field is zero. See figure 4, confirming second-order accuracy in $h=\Delta t$ for the MDVI and TDVI schemes. Interestingly, the error in the MDVI scheme is a factor of about $3$ larger, in contrast with the results shown in figure 2, where the MDVI and TDVI results are reversed. As for the magnetic field line results, this difference is due to the very long run times.

Figure 4.

Error in the second-order MDVI and TDVI schemes and the first-order collocated scheme of Ellison et al. (Reference Ellison, Finn, Burby, Kraus, Qin and Tang2018), compared with results obtained with a fourth-order Runge–Kutta method with extremely small time step $h=\Delta t$. Results are shown over almost two orders of magnitude in $\Delta t =h$. First-order and second-order reference lines, in blue dashed and green dashed, respectively, are shown for comparison. Second-order accuracy for the MDVI and TDVI schemes is confirmed.

6.1. Canonical systems

We first review the extended phase-space method in canonical variables, in one degree of freedom for transparency. The time stepping is defined by a time-step density $\rho (q,p)$, with $\rho (q,p)\,{\rm d}t={\rm d}\zeta$. We prescribe uniform steps in the time-like variable $\zeta$, with $\rho \Delta t\approx \Delta \zeta =h=\text {const}$. We extend the action $S=\int (p\,{\rm d}q/\,{\rm d}t-H)\,{\rm d}t$ for the Hamiltonian $H=H(q,p)$ to

(6.1)\begin{equation} S_{e}=\int\left[p\frac{{\rm d}q}{{\rm d}\zeta}-\frac{H(q,p)}{\rho(q,p)}+{\rm \pi}\left(\frac{{\rm d}w}{{\rm d}\zeta}-\frac{1}{\rho(q,p)}\right)\right]\,{\rm d}\zeta. \end{equation}

Here we have written time as the dependent variable $w$, and the Lagrange multiplier ${\rm \pi}$ enforces the time-step condition as a constraint. We rewrite this as

(6.2)\begin{equation} S_{e}=\int\left[p\frac{{\rm d}q}{{\rm d}\zeta}+{\rm \pi}\frac{{\rm d}w}{{\rm d}\zeta}-\frac{H(q,p)+{\rm \pi}}{\rho(q,p)}\right]\,{\rm d}\zeta, \end{equation}

leading to the extended phase-space Hamiltonian $K(q,p,{\rm \pi} )=(H(q,p)+{\rm \pi} )/\rho (q,p)$, with an added canonically conjugate pair, $(q,p)\to (q,p,w,{\rm \pi} )$. It is clear that, since $K$ does not depend on $\zeta$ explicitly, $K$ is exactly conserved. That is, if we set ${\rm \pi} =-H$ initially, then $H+{\rm \pi}$ remains exactly zero. Thus, the equations for $(q,p,w,{\rm \pi} )$, using $H+{\rm \pi} =0$, are

(6.3a,b)\begin{equation} \frac{{\rm d}q}{{\rm d}\zeta}=\frac{H_{2}(q,p,w)}{\rho(q,p)},\quad\frac{{\rm d}p}{{\rm d}\zeta}={-}\frac{H_{1}(q,p,w)}{\rho(q,p)}, \end{equation}

(6.4a,b)\begin{equation} \frac{{\rm d}w}{{\rm d}\zeta}=\frac{1}{\rho(q,p)},\quad \frac{{\rm d}{\rm \pi}}{{\rm d}\zeta}=0. \end{equation}

Note that the imposed time-step requirement is satisfied. The extension to a non-autonomous system, with $H(q,p,t)$, is straightforward.

For discrete integration, it is important to keep the terms (not shown) derived from $H+{\rm \pi}$ in (6.2): this quantity is not exactly zero for the discrete equations and if this quantity is set exactly equal to zero, the symplectic nature is lost, as discussed in Richardson & Finn (Reference Richardson and Finn2011).

To illustrate, we discretize $S_{e}$ as in the symplectic Euler scheme, here for two degrees of freedom, for uniform stepping in $\zeta$, $\Delta \zeta =h$. We have

(6.5)\begin{align} S_{e1}& =\sum_{k=0}^{N-1}\left[p_{k}(q_{k+1}-q_{k})+{\rm \pi}_{k}(w_{k+1}-w_{k})\right]\nonumber\\ & \quad -\sum_{k=0}^{N-1} \left[h\frac{H(q_{k+1},p_{k},w_{k+1})+{\rm \pi}_{k}}{\rho(q_{k+1},p_{k})}\right]. \end{align}

The symplectic Euler scheme derived from (6.5) preserves the canonical two-degree-of- freedom two-form

(6.6)\begin{equation} \omega={\rm d}q\wedge {\rm d}p+{\rm d}w\wedge {\rm d}{\rm \pi}, \end{equation}

by inspection or by taking the endpoint values of ${\rm d}S_{e1}$. This scheme is a single-step scheme, shown either by inspection or by computing the discrete Hessian.

We now consider the special class of systems in non-canonical variables with action of the form $\int (f_{i}\dot {x}_{i}-H(\boldsymbol {x},\boldsymbol {y}))\,{\rm d}t$, as in (3.1). With time-step condition $\rho (\boldsymbol {x},\boldsymbol {y})\,{\rm d}t={\rm d}\zeta$ and $t\to w$ we can write the analogue to (6.1) and (6.2):

(6.7)\begin{align} S& =\int\left[f_{i}(\boldsymbol{x},\boldsymbol{y})\frac{{{\rm d} x}^{i}}{{\rm d}\zeta}-\frac{H(\boldsymbol{x},\boldsymbol{y})}{\rho(\boldsymbol{x},\boldsymbol{y})} +{\rm \pi}\left(\frac{{\rm d}w}{{\rm d}\zeta}-\frac{1}{\rho(\boldsymbol{x},\boldsymbol{y})}\right) \right]\,{\rm d}\zeta,\nonumber\\ & =\int\left[f_{i}(\boldsymbol{x},\boldsymbol{y})\frac{{{\rm d} x}^{i}}{{\rm d}\zeta}+{\rm \pi}\frac{{\rm d}w}{{\rm d}\zeta}-\frac{H(\boldsymbol{x},\boldsymbol{y},w)+{\rm \pi}}{\rho(\boldsymbol{x},\boldsymbol{y})} \right]\,{\rm d}\zeta. \end{align}

Again, it is evident in the first form that ${\rm \pi}$ is a Lagrange multiplier enforcing the time-step condition. We can apply uniform stepping in the new time-like variable $\zeta$, with $\Delta \zeta =h={\rm const.}$, to give the required non-uniform time stepping. In (6.7) a term involving the canonical pair $(w,{\rm \pi} )$, namely ${\rm \pi} {\rm d}w/{\rm d}\zeta$, is added to $f_{i}\,{{\rm d} x}^{i}/{\rm d}\zeta$. Furthermore, the new Hamiltonian is $K=(H+{\rm \pi} )/\rho$, and the resulting action is of the same restricted non-canonical class of (3.1). Therefore any discretization that can be applied to the action with uniform time step can be applied to this non-canonical extended phase-space version.

6.2. Magnetic field line integration

For the magnetic field line integration problem, we use a gauge with $A_{1}=0$, as in § 4.1. For non-uniform time stepping, we first go back to considering $x^3$ to be a coordinate and put the action in the form in (6.7), leading to

(6.8)\begin{equation} \varPhi_{e}=\int\left[A_{2}(\boldsymbol{x})\frac{{{\rm d} x}^2}{{\rm d}\zeta}+\frac{A_{3}(\boldsymbol{x})}{\rho(\boldsymbol{x})}+{\rm \pi}\left(\frac{{{\rm d} x}^3}{{\rm d}\zeta}-\frac{1}{\rho(\boldsymbol{x})}\right)\right]\,{\rm d}\zeta, \end{equation}

where ${\rm \pi}$ is a Lagrange multiplier enforcing the time-step restriction and $\boldsymbol {x}=(x^1,x^2,x^3)$. This can again be put into the form

(6.9)\begin{equation} \varPhi_{e}=\int\left[A_{2}(\boldsymbol{x})\frac{{{\rm d} x}^2}{{\rm d}\zeta}+{\rm \pi}\frac{{{\rm d} x}^3}{{\rm d}\zeta}-\frac{-A_{3}(\boldsymbol{x})+{\rm \pi}}{\rho(\boldsymbol{x})}\right]\,{\rm d}\zeta. \end{equation}

This action is of the form in (6.7), with two-degree-of-freedom Hamiltonian equal to $(-A_{3}+{\rm \pi} )/\rho$.

We discretize this in a manner similar to the modified (adjoint) symplectic Euler scheme described in (3.3a,b) by forming $\varPhi _{e1}=\sum _{k}hL_{d}(\boldsymbol {x}_{k},{\rm \pi} _{k},\boldsymbol {x}_{k+1},{\rm \pi} _{k+1})$, or

(6.10)\begin{gather} \varPhi_{e1}=\sum_{k}\left[A_{2}(\boldsymbol{x}_{k+1})(x^2_{k+1}-x^2_{k})+{\rm \pi}_{k+1}(x^3_{k+1}-x^3_{k})\right]\nonumber\\ +\sum_k \left[h\frac{A_{3}(\boldsymbol{x}_{k+1})-{\rm \pi}_{k+1}}{\rho(\boldsymbol{x}_{k+1})}\right]. \end{gather}

Compared with (3.3a,b), $p_{k+1}\to A_{2}(\boldsymbol {x}_{k+1})$ and $H(q_{k+1},p_{k+1})\to A_{3}(\boldsymbol {x}_{k+1})$ with more direct substitutions for $x^3_{k}$ and ${\rm \pi} _{k}$.

The DEL equations from $d\varPhi _{e1}=0$ lead to the fourth-order system

(6.11)\begin{align} & A_{2,1}(\boldsymbol{x}_{k})(x^2_{k}-x^2_{k-1})+h\frac{A_{3,1}(\boldsymbol{x}_{k})}{\rho(\boldsymbol{x}_{k})}\nonumber\\ & \quad +h\left[A_{3}(\boldsymbol{x}_{k})-{\rm \pi}_{k}\right]\frac{\partial}{\partial x^1_{k}}\frac{1}{\rho(\boldsymbol{x}_{k})}=0, \end{align}

(6.12)\begin{align} & A_{2,2}(\boldsymbol{x}_{k})(x^2_{k}-x^2_{k-1})+A_{2}(\boldsymbol{x}_{k})-A_{2}(\boldsymbol{x}_{k+1})\nonumber\\ & \quad +h\frac{A_{3,2}(\boldsymbol{x}_{k})}{\rho(\boldsymbol{x}_{k})}+h\left[A_{3}(\boldsymbol{x}_{k})-{\rm \pi}_{k}\right]\frac{\partial}{\partial x^2_{k}}\frac{1}{\rho(\boldsymbol{x}_{k})}=0, \end{align}

(6.13)\begin{align} & A_{2,3}(\boldsymbol{x}_{k})(x^3_{k}-x^3_{k-1})+{\rm \pi}_{k}-{\rm \pi}_{k+1}+h\frac{A_{3,3}(\boldsymbol{x}_{k})}{\rho(\boldsymbol{x}_{k})}\nonumber\\ & \quad +h\left[A_{3}(\boldsymbol{x}_{k})-{\rm \pi}_{k}\right]\frac{\partial}{\partial x^3_{k}}\frac{1}{\rho(\boldsymbol{x}_{k})}=0, \end{align}

(6.14)\begin{align} & x^3_{k}-x^3_{k-1}-\frac{h}{\rho(\boldsymbol{x}_{k})}=0. \end{align}

As in earlier sections, the discrete Hessian test indicates that these four equations define a single-step scheme, and therefore a mapping on a four-dimensional space. The process of incrementing and substitution leads to equations for which this single-step property is explicit. For the uniform time-stepping case of § 4.1 the assumption $\rho =\rho (\boldsymbol {x})$ means that the time step depends on $x^3$, which takes the place of time. This suggests the possibility that complications such as parametric instabilities related to having a time-step density explicitly dependent on time might arise, as in Richardson & Finn (Reference Richardson and Finn2011). In the case treated in this subsection, on the other hand, $\zeta$ rather than $z$ is the independent (time-like) variable.

As occurred in Ellison et al. (Reference Ellison, Finn, Burby, Kraus, Qin and Tang2018), (6.12) and (6.13) appear to involve indices $k-1, k,k+1$ and therefore appear to be two-step equations, suggesting that (6.11)–(6.14) are difference equations of order higher than $4$. However, the discrete Hessian can be shown to have rank $4$, consistent with a first-order system in $(\boldsymbol {x},{\rm \pi} )$; indeed similar substitutions to those of Ellison et al. (Reference Ellison, Finn, Burby, Kraus, Qin and Tang2018) lead to a fourth-order system, i.e. a single-step method. That is, writing (6.12), (6.11) and (6.13) in the compact form

(6.15)\begin{gather} A_{2,1}(\boldsymbol{x}_{k})(x^2_{k}-x^2_{k-1})+hP(\boldsymbol{x}_{k},{\rm \pi}_{k})=0, \end{gather}

(6.16)\begin{gather}A_{2,2}(\boldsymbol{x}_{k})(x^2_{k}-x^2_{k-1})+A_{2}(\boldsymbol{x}_{k})-A_{2}(\boldsymbol{x}_{k+1})+hQ(\boldsymbol{x}_{k},{\rm \pi}_{k})=0, \end{gather}

(6.17)\begin{gather}A_{2,3}(\boldsymbol{x}_{k})(x^2_{k}-x^2_{k-1})+{\rm \pi}_{k}-{\rm \pi}_{k+1}+hR(\boldsymbol{x}_{k},{\rm \pi}_{k})=0, \end{gather}

we find $x^2_{k}-x^2_{k-1}$ from (6.15) and substitute into (6.16) and (6.17). When the indices are incremented $k\rightarrow k+1$ in (6.15) and (6.14), we find

(6.18)\begin{equation} \left. \begin{aligned} & A_{2,1}(\boldsymbol{x}_{k+1})(x^2_{k+1}-x^2_{k})+hP(\boldsymbol{x}_{k+1},{\rm \pi}_{k+1})=0,\\ & -hA_{2,2}(\boldsymbol{x}_{k})\frac{P(\boldsymbol{x}_{k},{\rm \pi}_{k})}{A_{2,1}(\boldsymbol{x}_{k})}+A_{2}(\boldsymbol{x}_{k})-A_{2}(\boldsymbol{x}_{k+1})+hQ(\boldsymbol{x}_{k},{\rm \pi}_{k})=0,\\ & -hA_{2,3}(\boldsymbol{x}_{k})\frac{P(\boldsymbol{x}_{k},{\rm \pi}_{k})}{A_{2,1}(\boldsymbol{x}_{k})}+{\rm \pi}_{k}-{\rm \pi}_{k+1}+hR(\boldsymbol{x}_{k},{\rm \pi}_{k})=0,\\ & x^3_{k+1}-x^3_{k}-\frac{h}{\rho(\boldsymbol{x}_{k+1})}=0. \end{aligned} \right\} \end{equation}

Because only indices $k$ and $k+1$ are involved, the single-step property predicted by the rank of the discrete Hessian is evident. Note the solvability condition $A_{2,1}=B^{3}\neq 0$, necessary for $z$ to parameterize the length along the field line. We have performed numerical tests with non-uniform time stepping, showing comparable properties with uniform time stepping of the field line equations (§ 5.2). We defer computations based on a time step adapted to an error estimator obtained by backward error analysis, as in Richardson & Finn (Reference Richardson and Finn2011), to a future publication.

6.3. Guiding centre equations

Guided by the results for the magnetic field line equations, it was shown in Ellison et al. (Reference Ellison, Finn, Burby, Kraus, Qin and Tang2018) that if the term $A_{1}^{{\dagger} }\dot {x}^1$ is zero, substitutions can be made to lead to a single-step scheme. Because the requirement $A_{1}^{{\dagger} }(\boldsymbol {x},u)=A_{1}(\boldsymbol {x})+ub_{1}(\boldsymbol {x})=0$ must hold for arbitrary values of $u$, it requires both the gauge condition $A_{1}=0$ and the physical condition $b_{1}=0$. In Ellison et al. (Reference Ellison, Finn, Burby, Kraus, Qin and Tang2018), numerical tests were performed for axisymmetric fields and for coordinates such that the covariant component $b_{1}$ is zero. In Burby & Ellison (Reference Burby and Ellison2017), it was shown that it is possible to find coordinates such that this condition is satisfied for arbitrary magnetic fields, provided one component, e.g. the toroidal component, does not change sign.

For applying non-uniform time steps to the guiding centre equations, we again introduce a new time-like variable $\zeta$ such that $\rho (\boldsymbol {x},u)\,\textrm {d}t=\textrm {d}\zeta$ for the time-step density $\rho (\boldsymbol {x},u)$. Then we make the substitution, with $w=t$ and

(6.19)\begin{equation} S_{{\rm gc}}=\int\left[A_{2}^{{{\dagger}}}(\boldsymbol{x})\frac{{{\rm d} x}^2}{{\rm d}t}+A_{3}^{{{\dagger}}}(\boldsymbol{x})\frac{{{\rm d} x}^3}{{\rm d}t}-H_{{\rm gc}}(\boldsymbol{x},u)\right]\,{\rm d}t \end{equation}

going to

(6.20)\begin{gather} S_{{\rm gc}}=\int\left[A_{2}^{{{\dagger}}}(\boldsymbol{x})\frac{{{\rm d} x}^2}{{\rm d}\zeta}+A_{3}^{{{\dagger}}}(\boldsymbol{x})\frac{{{\rm d} x}^3}{{\rm d}\zeta}-\frac{H_{{\rm gc}}(\boldsymbol{x},u)}{\rho(\boldsymbol{x},u)}\right] \,{\rm d}\zeta\nonumber\\ +\int \left[\frac{{\rm d}w}{{\rm d}\zeta}-\frac{1}{\rho(\boldsymbol{x},u)}\right]\,{\rm d}\zeta. \end{gather}

Again this can be put in the form

(6.21)\begin{equation} S_{{\rm gc}}=\int\left[A_{2}^{{{\dagger}}}(\boldsymbol{x})\frac{{{\rm d} x}^2}{{\rm d}\zeta}+A_{3}^{{{\dagger}}}(\boldsymbol{x})\frac{{{\rm d} x}^3}{{\rm d}\zeta} +{\rm \pi}\frac{{\rm d}w}{{\rm d}\zeta}-K(\boldsymbol{x},u,{\rm \pi})\right]\,{\rm d}\zeta, \end{equation}

where $K(\boldsymbol {x},u,p_{0})=(H_{\textrm {gc}}(\boldsymbol {x},u)+{\rm \pi} )/\rho (\boldsymbol {x},u)$. This is an action in the form of (6.7) on the extended phase space $(\boldsymbol {x},u)\rightarrow (\boldsymbol {x},u,w,{\rm \pi} )$. We discretize this action in a manner similar to that in (6.10), namely

(6.22)\begin{equation} S_{{\rm gc}1}=\sum_{k}hL_{d}(\boldsymbol{x}_{k},u_{k},w_{k},{\rm \pi}_{k},\boldsymbol{x}_{k+1},u_{k+1},w_{k+1},{\rm \pi}{}_{k+1}),\end{equation}

with

(6.23)\begin{gather} hL_{d}=A_{2}^{{{\dagger}}}(\boldsymbol{x}_{k+1},u_{k+1})(x^2_{k+1}-x^2_{k})+A_{z}^{{{\dagger}}}(\boldsymbol{x}_{k+1},u_{k+1})(x^3_{k+1}-x^3_{k})\nonumber\\ +{\rm \pi}_{k+1}(w_{k+1}-w_{k})-hK(\boldsymbol{x}_{k+1},u_{k+1},{\rm \pi}_{k+1}), \end{gather}

where

(6.24)\begin{equation} K(\boldsymbol{x}_{k+1},u_{k+1},{\rm \pi}_{k+1})=\frac{H_{{\rm gc}}(\boldsymbol{x}_{k+1},u_{k+1})+{\rm \pi}_{k+1}}{\rho(\boldsymbol{x}_{k+1},u_{k+1})}. \end{equation}

The DEL equations are

(6.25)\begin{gather} A_{2,1}^{{{\dagger}}}(\boldsymbol{x}_{k},u_{k})(x^2_{k}-x^2_{k-1})+A_{3,1}^{{{\dagger}}}(\boldsymbol{x}_{k},u_{k})(x^3_{k}-x^3_{k-1}), \end{gather}

(6.26)\begin{gather} -hK_{1}(\boldsymbol{x}_{k},u_{k},{\rm \pi}_{k})=0, \end{gather}

(6.27)\begin{gather} A_{2,2}^{{{\dagger}}}(\boldsymbol{x}_{k},u_{k})(x^2_{k}-x^2_{k-1})+A_{3,2}^{{{\dagger}}}(\boldsymbol{x}_{k},u_{k})(x^3_{k}-x^3_{k-1})\nonumber\\ +A_{2}^{{{\dagger}}}(\boldsymbol{x}_{k},u_{k})-A_{2}^{{{\dagger}}}(\boldsymbol{x}_{k+1},u_{k+1})-hK_{2}(\boldsymbol{x}_{k},u_{k},{\rm \pi}_{k})=0, \end{gather}

(6.28)\begin{gather} A_{2,3}^{{{\dagger}}}(\boldsymbol{x}_{k},u_{k})(x^2_{k}-x^2_{k-1})+A_{3,3}^{{{\dagger}}}(\boldsymbol{x}_{k},u_{k})(x^3_{k}-x^3_{k-1})\nonumber\\ +A_{3}^{{{\dagger}}}(\boldsymbol{x}_{k},u_{k})-A_{3}^{{{\dagger}}}(\boldsymbol{x}_{k+1},u_{k+1})-hK_{3}(\boldsymbol{x}_{k},u_{k},{\rm \pi}_{k})=0, \end{gather}

(6.29)\begin{gather} b_{2,k}(x^2_{k}-x^2_{k-1})+b_{3,k}(x^3_{k}-x^3_{k-1})-hK_{u}(\boldsymbol{x}_{k},u_{k},{\rm \pi}_{k})=0,\end{gather}

(6.30)\begin{gather} {\rm \pi}_{k}-{\rm \pi}_{k+1}=0, \end{gather}

(6.31)\begin{gather} w_{k}-w_{k-1}=h\frac{\partial}{\partial p_{w,k}}K(\boldsymbol{x}_{k},u_{k},{\rm \pi}_{k})\nonumber\\ =\frac{h}{\rho(\boldsymbol{x}_{k},u_{k})}. \end{gather}

The assumed time independence of the fields leads to the simple form in (6.30). Similar to the uniform time-step case of Ellison et al. (Reference Ellison, Finn, Burby, Kraus, Qin and Tang2018), the discrete Hessian has rank six, consistent with a first-order system in $(x,y,z,u,w,{\rm \pi} )$. We start by taking (6.26) and (6.29), written as

(6.32)\begin{equation} \left[\begin{array}{@{}cc@{}} A_{2,1}^{{{\dagger}}}(\boldsymbol{x}_{k},u_{k}) & A_{3,1}^{{{\dagger}}}(\boldsymbol{x}_{k},u_{k})\\ b_{2,k} & b_{3,k} \end{array}\right]\left[\begin{array}{@{}c@{}} x^2_{k}-x^2_{k-1}\\ x^3_{k}-x^3_{k-1} \end{array}\right] =\left[\begin{array}{@{}c@{}} hK_{1}(\boldsymbol{x}_{k},u_{k},w_{k},{\rm \pi}_{k})\\ hK_{u}(\boldsymbol{x}_{k},u_{k},w_{k},{\rm \pi}_{k}) \end{array}\right]. \end{equation}

Solving for $x^2_{k}-x^2_{k-1}$ and $x^3_{k}-x^3_{k-1}$, which are written in terms of quantities with index $k$, we substitute these into (6.27) and (6.28), incrementing $k\rightarrow k+1$ in (6.27), (6.28) and (6.31). The resulting equations involve time steps labelled with only $k$ and $k+1$. That is, consistent with the discrete Hessian condition, the scheme is a single-step scheme, a DVI, and parasitic modes cannot occur. As for the magnetic field line case discussed in the last subsection, we have obtained similar numerical results to uniform time-stepping results of § 5.3, and defer adaptive integration to a future publication.

6.4. Extensions for higher accuracy

From the formulation in the last two sections, it is clear from (6.9) and (6.21) that the modification to prescribe non-uniform time stepping leads to an addition to the phase-space Lagrangian of a term ${\rm \pi} dw/\textrm {d}\zeta$ or ${\rm \pi} \,\textrm {d}z/\textrm {d}\zeta$ and a modification to the Hamiltonian $H\to (H+{\rm \pi} )/\rho$, and these terms can be discretized in exactly the same manner as in the uniform time-step case. This means that the modifications in this section can be applied to any discretization of the phase-space Lagrangian that leads to a DVI scheme. Therefore, it should be straightforward to construct a non-uniform time-step scheme for either of the second-order-accurate DVI methods of § 3.

It is also clear that such discretizations can be applied to any time-step density $\rho$, so that it should be straightforward to use an optimum density $\rho$ based on an error estimator, to minimize the integrated error over an orbit for the scheme at hand, as done in Richardson & Finn (Reference Richardson and Finn2011) and Finn (Reference Finn2015). Therefore, it is possible to combine the formulations of this paper to give an adaptive second-order-accuratevariational integrator. We leave further details to a future publication.