Liouville Formalism: a Tool for Building Symplectic and Reversible Integrators

In the previous paragraphs we have seen that it is highly beneficial for an integrator to be symplectic. We may now wonder if there exists a general way for obtaining symplectic and possibly, reversible integrators from ``first principles''. To this end, we start by noting that for any property which depends on time implicitly through $p,q \equiv {\bf x}$ we have

$${dA(p,q) \over dt} = \sum_{q,p} \left ( \dot q {\partial A \over \partial q} + \dot p ... ...rtial q} - {\partial H \over \partial q} {\partial A \over \partial p} \right )$$

$$= iL A$$ (2.13)

where the sum is extended to all $n$ degrees of freedom in the system. $L$ is the Liouvillean operator defined by

$$iL \equiv \sum_{q,p} \left ( \dot q {\partial \over \partial q} + \dot p {\... ...artial q} - {\partial H \over \partial q} {\partial \over \partial p} \right ).$$ (2.14)

Eq. (2.13) can be integrated to yield

$$A(t) = e^{iL t} A(0).$$ (2.15)

If $A$ is the state vector itself we can use Eq. (2.15) to integrate Hamilton's equations:

$$\left [ \begin{array}{r} q(t) \\ p(t) \end{array}\right ] = e^{iLt} \left [ \begin{array}{r} q(0) \\ p(0) \end{array}\right ].$$ (2.16)

The above equation is a formal solution of Hamilton's equations of motion. The exponential operator $e^{iLt}$ times the state vector defines the t-flow of the Hamiltonian system which brings the system phase space point from the initial state $q_{0},p_{0}$ to the state $p(t),q(t)$ at a later time $t$. We already know that this transformation obeys Eq. (2.8). We may also note that the adjoint of the exponential operator corresponds to the inverse, that is $e^{iLt}$ is unitary. This implies that the trajectory is exactly time reversible. In order to build our integrator, we now define the discrete time propagator $e^{iL\Delta t}$ as

$$e^{iLt} = \left [ e^{iLt/n} \right ]^{n} ; ~~~~~~~~~~~ \Delta t = t/n$$ (2.17)

$$e^{iL\Delta t} = e^{\sum_{q,p} \left ( \dot q {\partial \over \partial q} + \dot p {\partial \over \partial p} \right )\Delta t }.$$ (2.18)

In principle, to evaluate the action of $e^{iL\Delta t}$ on the state vector $p,q$ one should know the derivatives of all orders of the potential $V$. This can be easily seen by Taylor expanding the discrete time propagator $e^{il\Delta t}$ and noting that the operator $\dot q \partial / \partial q$ does not commute with $- \partial V /\partial q (\partial /\partial p)$ when the coordinates and momenta refer to same degree of freedom. We seek therefore approximate expressions of the discrete time propagator that retain both the symplectic and the reversibility property. For any two linear operators $A,B$ the Trotter formula [74] holds:

$$e^{(A+B) t} = \lim_{n \rightarrow \infty} (e^{ A t/n} e^{B t/n} )^{n}$$ (2.19)

We recognize that the propagator Eq. (2.18) has the same structure as the left hand side of Eq. (2.19); hence, using Eq. (2.19), we may write for $\Delta t$ sufficiently small

$$e^{iL\Delta t} = e^{\left ( \dot q {\partial \over \partial q} + \dot p {\partial ... ...a t } e^{\dot p {\partial \over \partial p} \Delta t } + O (\Delta t^{2}) ~~~~.$$ (2.20)

Where, for simplicity of discussion, we have omitted the sum over $q$ and $p$ in the exponential. Eq. (2.20) is exact in the limit that $\Delta t \rightarrow 0$ and is first order for finite step size. Using Eq. (2.8) it is easy to show that the t-flow defined in Eq. (2.20) is symplectic, being the product of two successive symplectic transformations. Unfortunately, the propagator Eq. (2.20) is not unitary and therefore the corresponding algorithm is not time reversible. Again the non unitarity is due to the fact that the two factorized exponential operators are non commuting. We can overcome this problem by halving the time step and using the approximant:

$$e^{(A+B)t} \simeq e^{At/2}e^{Bt/2}e^{Bt/2}e^{At/2} = e^{At/2}e^{Bt}e^{At/2} .$$ (2.21)

The resulting propagator is clearly unitary, therefore time reversible, and is also correct to the second order [75]. Thus, requiring that the product of the exponential operator be unitary, automatically leads to more accurate approximations of the true discrete time propagator [76,75]. Applying the same argument to the propagator (2.18) we have

$$e^{iL\Delta t} = e^{\left ( \dot q {\partial \over \partial q} + \dot p {\partial ... ...lta t } e^{\dot p {\partial \over \partial p} \Delta t/2 } + O (\Delta t^{3}) .$$ (2.22)

The action of an exponential operator $e^{(a \partial / \partial x )}$ on a generic function $f(x)$ trivially corresponds to the Taylor expansion of $f(x)$ around the point $x$ at the point $x+a$, that is

$$e^{a \partial / \partial x} f(x) = f(x+a) .$$ (2.23)

Using Eq. (2.23), the time reversible and symplectic integration algorithm can now be derived by acting with our Hermitian operator Eq. (2.22) onto the state vector at $t=0$ to produce updated coordinate and momenta at a later time $\Delta t$. The resulting algorithm is completely equivalent to the well known velocity Verlet:

$$p(\Delta t /2) = p(0) + F(0) \Delta t /2$$

$$q(\Delta t) = q(0) + \left ( {p(\Delta t /2) \over m } \right )\Delta t$$

$$p(\Delta t) = p(\Delta t/2) + F(\Delta t) \Delta t /2 .$$ (2.24)

We first notice that each of the three transformations obeys the symplectic condition Eq. (2.8) and has a Jacobian determinant equal to one. The product of the three transformation is also symplectic and, thus, phase volume preserving. Finally, since the discrete time propagator (2.22) is unitary, the algorithm is time reversible.

One may wonder what it is obtained if the operators $\dot q \partial / \partial q$ and $- \partial V /\partial q (\partial /\partial p)$ are exchanged in the definition of the discrete time propagator (2.22). If we do so, the new integrator is

$$q(\Delta t /2) = q(0) + { p(0) \over m} \Delta t /2$$

$$p(\Delta t) = p(0) + F[q (\Delta t /2)] \Delta t$$

$$q(\Delta t) = q(\Delta t/2) + { p(\Delta t) \over m} \Delta t /2 .$$ (2.25)

This algorithm has been proved to be equivalent to the so-called Leap-frog algorithm [77]. Tuckerman et al. [20] called this algorithm position Verlet which is certainly a more appropriate name in the light of the exchanged role of positions and velocities with respect to the velocity Verlet. Also, Eq. (2.21) clearly shows that the position Verlet is essentially identical to the Velocity Verlet. A shift of a time origin by $\Delta t/2$ of either Eq. (2.25) or Eq. (2.24) would actually make both integrator perfectly equivalent. However, as pointed out in Ref. [21], half time steps are not formally defined, being the right hand side of Eq. (2.21) an approximation of the discrete time propagator for the full step $\Delta t$. Velocity Verlet and Position Verlet, therefore, do not generate numerically identical trajectories although of course the trajectories are similar.

We conclude this section by saying that is indeed noticeable that using the same Liouville formalism different long-time known schemes can be derived. The Liouville approach represent therefore a unifying treatment for understanding the properties and relationships between stepwise integrators.