Error analysis of the transport properties of Metropolized schemes

Description of the model

Standard results (see for instance the references in [13, Section 2.2]) show that (1) admits the Boltzmann-Gibbs measure

as its unique invariant probability measure (note that $Z$ is finite since the position space $\mathcal{M}$ is compact and $V$ is smooth hence bounded). In fact, (1) is ergodic with respect to this measure, where ergodicity is understood both as (i) the long-time (almost-sure) convergence of averages along trajectories

for any initial condition $q_{0}\in\mathcal{M}$ and all observables $f\in L^{1}(\mu)$; or as (ii) the convergence of the law $\psi(t,q)\,dq$ of the process (1), happening here at an exponential rate, for instance in total variation: Denoting with some abuse of notation the measure $\psi(t,q)\,dq$ by $\psi(t)$, there exist $C,\lambda>0$ such that

where the total variation distance between two measures $\nu_{1},\nu_{2}$ is defined as

the suprema being taken over all measurable sets of $\mathcal{M}$ for the first one, and over all bounded, measurable functions for the second one.

For further purposes, we introduce the generator of (1), namely the operator

This operator (defined with domain $D(\mathcal{L})=H^{2}(\mu)$) is self-adjoint on the Hilbert space $L^{2}(\mu)$ endowed with the scalar product

The operator $-\mathcal{L}$ moreover has a positive spectral gap (see for instance [13, Section 2] and references therein). Indeed, a simple computation shows that

The Poincaré inequality $\|\varphi\|_{L^{2}(\mu)}\leqslant C_{\mathcal{M},V}\|\nabla\varphi\|_{L^{2}(\mu)}$, valid for any function belonging to

which shows that the spectral gap is larger or equal to $C_{\mathcal{M},V}^{-1}$. In particular, the resolvent $\mathcal{L}^{-1}$ is a well-defined operator on $\widetilde{L}^{2}(\mu)$, and the following estimate holds:

Here and in the sequel, for a given Banach space $X$, we denote by $\mathcal{B}(X)$ the Banach space of bounded operators on $X$, endowed with the norm

2. Self-diffusion

where the expectation is over all realizations of the continuous dynamics (1), starting from initial conditions distributed according to the Boltzmann-Gibbs measure (2). The following result shows that the diffusion tensor (8) is well defined, and naturally arises in a diffusive time-rescaling of the dynamics (1).

Consider for $\varepsilon>0$ the diffusively rescaled process $Q_{t}^{\varepsilon}=\varepsilon Q_{t/\varepsilon^{2}}$. Then, as $\varepsilon\to 0$, the process $Q_{t}^{\varepsilon}$ starting from a given initial condition $Q_{0}$ weakly converges on finite time intervals to an effective Brownian motion starting from $Q_{0}$ and with covariance matrix $\mathscr{D}$ given by (8). Moreover, $\mathscr{D}$ is a real, positive definite $dN\times dN$ matrix, satisfying

in the sense of symmetric matrices, and which can alternatively be expressed as

where the expectation is over all realizations of the continuous dynamics (1), starting from initial conditions distributed according to the Boltzmann-Gibbs measure (2).

The proof of this statement is standard, and follows from arguments presented in [1, Chapter 3] for instance. We nonetheless provide a short proof in Section 5.1 since the proofs of the discrete counterparts of Theorem 1 rely on an appropriate extension of the argument used in the continuous case (see Section 5.5).

A straightforward consequence of Theorem 1 is that the self-diffusion constant $\mathcal{D}$, defined as the average mean-square displacement of the individual particles, is well defined and has two equivalent expressions:

The expression (10) is called the Einstein formula. The second expression (11) involves an integrated autocorrelation function. In accordance with the standard physics and chemistry nomenclature, we call (11) the Green-Kubo formula for the self-diffusion in the sequel.

3. Numerical estimation of the self-diffusion

In order to compute approximations of formulas such as (10) or (11), the first task is to numerically integrate realizations of the continuous dynamics (1). The most straightforward way would be to resort to a Euler-Maruyama scheme: given a time-step $\Delta t>0$ and denoting by $q^{n}$ an approximation of $q_{n\Delta t}$, this scheme reads

where $(G^{n})_{n\geqslant 0}$ is a sequence of independent and identically distributed (i.i.d.) $dN$-dimensional standard Gaussian random variables. However, this simple scheme has been shown to fail to be ergodic when the dynamics is considered on unbounded spaces and the potential energy function is not globally Lipschitz . In simulations of Brownian dynamics for ionic solutions, potential energy functions with Coulomb-type singularities are used and it has been observed that the energy blows up along trajectories of (12) (see Section 3.2).

is the probability transition of the Markov chain (12). When the proposition is accepted, we project $\widetilde{q}^{n+1}$ into the periodic simulation cell $\mathcal{M}$. If the proposal is rejected, the previous configuration is counted twice: $q^{n+1}=q^{n}$ (It is very important to count rejected configuration as many times as needed to ensure that the Boltzmann-Gibbs measure $\mu$ is invariant). In conclusion,

where $U^{n}$ are i.i.d. uniform random variables in $$, and$\mathbf{1}_{U^{n}\leqslant\alpha}$is an indicator function whose value is 1 when$U^{n}\leqslant\alpha$ and 0 otherwise. The average rejection rate is

In order to avoid confusion, we call the scheme (13) “Metropolized Euler-Maruyama” in the sequel, and denote by $P_{\Delta t}$ its evolution operator:

Of course, the fact that some configurations are rejected destroys the trajectorial accuracy of the dynamics, see for precise statements. The resulting strong errors and, as a consequence, errors on finite time correlation functions have been quantified in , with prefactors which unfortunately depend on time. The estimates provided by these works therefore do not provide error estimates on diffusion coefficients, obtained either as infinite time integrals of correlation functions as in the Green-Kubo formula (11) or as the infinite time average mean square displacement as in Einstein’s formula (10).

The next section quantifies the errors in the approximation of (10) and (11) when the Metropolized Euler-Maruyama scheme is used. Although the formulas (10) and (11) are equivalent for continuous dynamics, they lead to different numerical methods. Let us already emphasize that the errors on the diffusion coefficients are in fact determined by the expansion of the evolution operator $P_{\Delta t}$. This expansion is the same as the one used to establish weak error estimates. From a technical viewpoint, the techniques used in the proofs of our main results are therefore quite different from the techniques of , which are based on strong error estimates obtained with Gronwall’s lemma.

A priori error estimates on the self-diffusion

As discussed in , error bounds on transport properties in fact depend on approximation properties of the evolution operator (similar to the one used to prove weak error estimates), rather than strong error estimates. A key building block in this framework is the following expansion of the evolution operator, obtained by a slight extension of [5, Lemma 4.7] and [4, Lemma 5.5].

There exist an operator $A$ and ${\Delta t}^{*}>0$ such that, for any $0<{\Delta t}\leqslant{\Delta t}^{*}$ and any smooth function $\psi$,

with a remainder $r_{\psi,{\Delta t}}$ uniformly bounded for $0<{\Delta t}\leqslant{\Delta t}^{*}$. Moreover,

where the expectation is over all possible realizations of the standard $dN$-dimensional Gaussian random variable $G$.

The precise expression of the operator $A$ is unimportant. It is however given in the proof of this result, see Section 5.2. Note that the numerical scheme can be proved to be weakly first order accurate by relying on standard techniques , in view of the equality

Another important result which we will repeatedly use in the analysis below is the following uniform-in-${\Delta t}$ geometric ergodicity of the Metropolized Euler-Maruyama scheme, easily obtained by adapting the results of to the case of compact position spaces (for completeness, we nonetheless provide elements of proof in Section 5.3). To state the result, we introduce the following functional space

As a consequence, there exists $K>0$ such that

We first give error estimates on (11) by appropriately adapting the results from . The result is stated for two smooth observables $\psi,\varphi$ with average 0 with respect to $\mu$. Define to this end

Error estimates for (11) are obtained by setting $\psi=\varphi=\partial_{q_{i,\alpha}}V$, with $1\leqslant i\leqslant N$ and $1\leqslant\alpha\leqslant d$. Note indeed that a simple integration by parts shows that $\partial_{q_{i,\alpha}}V$ has average 0 with respect to $\mu$, so $\partial_{q_{i,\alpha}}V\in\widetilde{C}^{\infty}(\mathcal{M})$.

Consider two observables $\psi,\varphi\in\widetilde{C}^{\infty}(\mathcal{M})$, and define the modified observable

where the operator $A$ is defined in (15). Then, there exists ${\Delta t}^{*}>0$ such that, for any $0<{\Delta t}\leqslant{\Delta t}^{*}$,

with $r_{\psi,\varphi,{\Delta t}}$ uniformly bounded (with respect to ${\Delta t}$), and where the expectation on the left hand side of the above equation is with respect to initial conditions $q_{0}\sim\mu$ and over all realizations of the dynamics (1), while the expectation on the right hand side is with respect to initial conditions $q^{0}\sim\mu$ and over all realizations of the Metropolized Euler-Maruyama scheme (13).

As a corollary, we obtain first order error bounds on the computation of the self-diffusion through (11):

where $\widetilde{\mathcal{D}}^{\rm GK}_{\Delta t}$ is uniformly bounded for ${\Delta t}$ sufficiently small, and where the numerically computed self-diffusion reads

The expression of the correction term is obtained by replacing the modified observable by its expression:

The appearance of subleading fractional correction term in (20) (here, of order ${\Delta t}^{3/2}$) is typical of Metropolis algorithms, and is usually not encountered for standard, un-Metropolized discretizations of SDEs (compare with the results of ).

2. Error bounds on the Einstein formula

In this section, we investigate the discretization error made when using the Metropolized Euler-Maruyama scheme to approximate the self-diffusion using (10). In accordance with the definition (7), we introduce a discrete additive functional allowing to keep track of the diffuse behavior of the Markov chain: Starting from $Q^{0}=q^{0}$,

where, as in Theorem 2, the expectation on the right hand side is with respect to initial conditions $Q^{0}=q^{0}\sim\mu$ and for all realizations of the Metropolized Euler-Maruyama scheme.

There exists ${\Delta t}^{*}>0$ such that, for any $0<{\Delta t}\leqslant{\Delta t}^{*}$,

The proof of this result can be read in Section 5.5. In fact, a slight extension of our technique of proof would allow to show that the diffusively rescaled process generated by the Metropolized Euler-Maruyama scheme, namely $\varepsilon Q^{\lfloor t/({\Delta t}\,\varepsilon^{2})\rfloor}$ (where $\lfloor x\rfloor$ denote the unique integer such that $\lfloor x\rfloor\leqslant x<\lfloor x\rfloor+1$), weakly converges on finite time intervals to a Brownian motion with covariance matrix $\mathscr{D}_{\Delta t}$.

An immediate corollary of Theorem 3 is the following a priori error estimate on the self-diffusion:

where $\widetilde{\mathcal{D}}^{\rm Einstein}_{\Delta t}$ is uniformly bounded for ${\Delta t}$ sufficiently small. Some more work would allow to prove that the subleading correction term is of order ${\Delta t}^{3/2}$, as in the Green-Kubo case (see Remark 4).

Numerical illustration

The aim of this section is to illustrate the errors bounds (20) and (24). We perform long computations so that the statistical errors are negligible. Let us mention that numerical simulations illustrating timestep errors for velocity autocorrelation functions were already presented in .

The self-diffusion coefficient $\mathcal{D}_{\Delta t}^{\rm Einstein}$ for the Einstein approach is approximated by fitting the unnormalized self-diffusion

by a linear function $\mathcal{D}^{{\rm Einstein},M}_{\Delta t}\,n{\Delta t}$, the slope being the estimation of the self-diffusion for the time step under consideration. This is indeed confirmed by Figure 1 (Left), which presents the evolution of $D_{n}^{M}$ as a function of the physical time $n{\Delta t}$. The results produced in Figures 1 and 3 have been obtained with $M=10^{7}$ replicas and $n_{\rm Einstein}=3\times 10^{5}$ steps. Let us also note that, in accordance with (17), the rejection rate scales as ${\Delta t}^{3/2}$.

A numerical approximation of (21) requires both a discretization using finitely many replicas, but also a truncation of the integration in time with an upper bound $\tau$. We consider the following numerical estimation of the self-diffusion coefficient obtained with the Green-Kubo formula:

The correlation functions shown in Figure 2 suggest that the autocorrelation of $V^{\prime}$ is exponentially decreasing. It can be considered as negligible for times larger than 0.2. The numerical results reported in Figure 2 and 3 have been obtained with $M=2\times 10^{8}$ replicas and a time cut-off $\tau=0.3$.

The results presented in Figure 3 indeed confirm that, for small time steps ${\Delta t}$, the error in the self-diffusion is of order ${\Delta t}$ for both methods. For larger time steps, nonlinearities appear. Note also that the errors on the coefficients computed with the Green-Kubo formula are smaller in this simple case. In any case, in accordance with the statements of Theorem 1, the self-diffusion is between 0 and 1 when ${\Delta t}\to 0$.

2. The more realistic case of solvated ions

We additionally consider a large ion, modeled as a fixed particle at position $q_{\rm ion}$ (the center of the simulation box), whose interaction with the solvent particles is described by an attractive Yukawa potential (screened Coulomb interaction) plus some repulsive potential preventing the small particles from coming too close to the ion. More precisely, for a solvent particle at position $q$, the interaction reads $v_{\rm ion}(|q-q_{\rm ion}|)$, with

The potentials $v$ and $v_{\rm ion}$ are plotted in Figure 4.

Expectations are approximated using $M$ trajectories of the system. We integrate trajectories one after the other, using the Metropolized Euler scheme (13), with initial conditions for the $(m+1)$th trajectory obtained by taking the last configuration of the $m$th trajectory. The self-diffusion coefficient calculated with the Green-Kubo formula is obtained as in the previous section, using the estimator (26) (upon introducing the correct normalization factor $1/(3N)$ in the autocorrelation and replacing $V^{\prime}$ by $\nabla V$). The values obtained by the Einstein formula are computed by dividing the unnormalized diffusion $D_{n}^{M}$ defined in (25) by the final time of the simulation: for a simulation time $\tau$,

The results presented in Figure 5 show that the unnormalized mean squared displacement indeed grows linearly in time, as expected. Note also that, in accordance with (17), the rejection rate scales as ${\Delta t}^{3/2}$. The results of Figure 5 and (7) have been obtained with $M=10^{5}$ trajectories, with integrations performed up to $\tau=20$.

The results presented in Figure 6 suggest that the decay of the force autocorrelation cannot be represented by a single exponential function. The force autocorrelation can be considered to be small in relative value for times of the order of $0.1$. The numerical results reported in Figure 6 and 7 were obtained by averaging $M=10^{6}$ trajectories with an integration time $\tau=0.3$.

Error estimates for the diffusion coefficients are gathered in Figure 7. For small time steps $\Delta t$, the error in the self-diffusion is linear in $\Delta t$, while nonlinearities appear for larger time steps. The continuous lines are linear fits obtained over the values corresponding to the 10 smallest time steps. As in the simple example discussed in the previous section, estimates obtained with Green-Kubo’s formula seem more reliable than those obtained with Einstein’s formula. Note also that, in accordance with Theorem 1, the self-diffusion is between 0 and 1 in all cases for sufficiently small time-steps.

Possible work tracks to reduce the error on the estimation of the self-diffusion

Both the Green-Kubo and the Einstein approaches lead to discretization errors of order ${\Delta t}$, as proved theoretically and verified numerically. A natural question is how to reduce this error. In the chemistry literature, it was proposed in to renormalize the time in Einstein’s method by replacing the simulation time $n{\Delta t}$ appearing in the denominator of (23) by some effective time $\theta_{\Delta t}n{\Delta t}$, where $\theta_{\Delta t}$ is the average acceptance rate of the Metropolis algorithm for a given time step ${\Delta t}$. However, since the average acceptance rate is of order $1-C{\Delta t}^{3/2}$ for small time steps (see (17)), such a correction cannot possibly cancel out the ${\Delta t}$-error in the diffusion coefficient.

A more promising work track, which we started working on at the end of our stay at CIRM, is to modify the proposed move (instead of simply considering the Euler-Maruyama scheme (12)), and possibly the invariant measure as well, in order to increase the weak and strong orders of the associated Metropolized scheme. Some steps in this direction have already been pursued in . In fact, the proofs of Theorem 2 and 3 show that, in order to gain accuracy on the computation of transport coefficients, it is sufficient to find a numerical scheme such that

A key element to obtain such equalities is to decrease the rejection rate for small time steps.

Proof of the results

In all the proofs, the constants $C>0$ and the critical time steps ${\Delta t}^{*}$ may change from line to line. Upon changing $V$ into $\beta V$, we may also assume that $\beta=1$.

The idea is to rewrite the part of the additive functional involving $-\nabla V(q_{t})$ as an approximate martingale. To this end, we introduce the solutions $\Phi_{0}=(\Phi_{0,1},..,\Phi_{0,dN})$ of the following Poisson equations

In view of the resolvent estimate (6), the functions $\Phi_{0,j}$ are well defined elements of $\widetilde{L}^{2}(\mu)$ since

In addition, by elliptic regularity, the functions $\Phi_{0,j}$ are smooth. By Ito’s lemma, we therefore obtain

so that the integrated displacement from the origin can be rewritten as

The ergodic properties of the diffusion process $q_{t}$ allow to prove that

Since $-\mathcal{L}$ has a positive spectral gap on $L^{2}(\mu)$ (see (5)), we can write the following operator equality on $\widetilde{L}^{2}(\mu)$ :

Therefore, for general functions $\psi,\varphi$ with vanishing average with respect to $\mu$,

Combining this result with (32) leads to the expression (9) of the diffusion matrix $\mathscr{D}$. To prove the convergence of the processes, two arguments should be made precise (see for an elementary account):

the convergence of the finite-dimensional laws, which can be obtained very simply here by considering the exponential martingales

which are such that the conditional expectations converge to those of a Brownian motion as $\varepsilon\to 0$:

$\mathcal{F}_{s/\varepsilon^{2}}$ denoting the filtration of events until the time $s/\varepsilon^{2}$. Finite-dimensional laws are then obtained by a simple induction, as made precise in for instance.

the tightness of the process, proved using Prohorov’s criterion (see for instance ):

This criterion is satisfied in view of the tightness of the martingale $\mathscr{M}^{F,\varepsilon}_{t}$, itself easily obtained using Doob’s inequality (see ).

2. Proof of Lemma 1

We first determine the magnitude of the acceptance probability in the Metropolis algorithm, which reads

Using the following expansions (with integral remainders)

while there exists a constant $C>0$ such that $\left|\widetilde{\xi}_{\Delta t}(q,G)\right|\leqslant C(1+|G|^{6})$ for any $0\leqslant{\Delta t}\leqslant 1$. We next use the inequality

obtained by distinguishing the cases $x\leqslant 0$ and $x\geqslant 0$. This shows that

with $\xi_{+}(q,G)=\max(0,\xi(q,G))$ and where $\left|\widehat{\xi}_{\Delta t}(q,G)\right|\leqslant C(1+|G|^{12})$. The estimate on the average acceptance rate (17) is obtained by taking the expectation over all possible realizations of $G^{n}$, with the definition

To obtain the action of the operator $A$, we start from the expression of the Metropolis transition operator (see for instance [13, Section 2.1.2])

We now write $q^{\prime}=q-{\Delta t}\,\nabla V(q)+\sqrt{2{\Delta t}}\,g$, so that, using (37) to estimate the rejection rate,

where we have used for the first integral a Taylor expansion at fourth order similar to (34) to obtain (see the computations in [12, Section 4.9])

and a Taylor expansion at first order for the term involving the rejection rate to obtain

The remainder $r_{\psi,{\Delta t}}$ is uniformly bounded in $L^{\infty}(\mathcal{M})$ for ${\Delta t}$ sufficiently small. The conclusion follows by setting $A=(\mathcal{L}^{2}+D_{1}+D_{2})/2$.

Finally, the invariance of $\mu$ by $P_{\Delta t}$ implies

This equality, together with the expansion (15) proves (16).

3. Proof of Lemma 2

The proof below is a simplification of the argument presented in , made possible since we work on a compact state space. The idea of the proof is to compare the Metropolis dynamics to the continuous dynamics, using the standard, un-Metropolized Euler-Maruyama scheme as an intermediate. Alternatively, it would be possible to directly compare the Metropolized and un-Metropolized schemes, by proving as in [12, Section 4.2] that the standard, un-Metropolized Euler-Maruyama scheme is geometrically ergodic since it satisfies a minorization condition.

Since we work on a compact state space, it is enough to show by Harris’ theorem (see the presentation in ) that there exist $\alpha\in(0,1)$ such that

uniformly in $0<{\Delta t}\leqslant{\Delta t}^{*}$ and $(q,q^{\prime})\in\mathcal{M}^{2}$. We now introduce the transition kernel $Q_{t}$ of the continuous dynamics (1), defined as

and consider it at time $t=1$. The transition kernel is well defined and regular since the generator is elliptic. By the triangle inequality,

From [4, Lemma 2.7], since we work on a compact space, we know that there exists $\varepsilon>0$ such that

To control the second term in (39), we introduce the transition kernel of the standard, un-Metropolized Euler-Maruyama scheme (12), denoted by $\widetilde{P}_{\Delta t}$, and write

By [4, Lemma 4.2], the transition kernel of the standard, un-Metropolized dynamics is uniformly close to the transition kernel of the continuous dynamics when the state space is compact: There exists $C_{\rm EM}>0$ and ${\Delta t}^{*}>0$ such that, for any $0<{\Delta t}\leqslant{\Delta t}^{*}$,

It therefore remains to control the distance between the transition rate of the Metropolized and un-Metropolized dynamics. It is at this stage that the argument of [4, Lemma 4.6] can be simplified. As in the proof of this lemma, we use a coupling argument, and consider two chains $q^{i}$ and $\widetilde{q}^{i}$, corresponding respectively to the Metropolized dynamics and the standard un-Metropolized one, starting from the same initial condition $q^{0}\in\mathcal{M}$. The probability that $q^{n}\neq\widetilde{q}^{n}$ is bounded from above by the probability that $q^{i}\neq\widetilde{q}^{i}$ for some $1\leqslant i\leqslant n$, i.e. at least one rejection occurred along the discrete trajectory. Since the probability to reject the move from $q^{i}$ to $q^{i+1}$ is $1-R_{\Delta t}(q^{i},q^{i+1})$, it holds

The combination of all previous estimates finally gives (38) provided ${\Delta t}$ is sufficiently small.

from which the bound (19) immediately follows.

4. Proof of Theorem 2

We follow the strategy of [12, Section 4.8]. The proof starts by noticing that the integrated correlation function can be written using $\mathcal{L}^{-1}$, as made precise in (33). The strategy of the proof is to write an approximation of $\mathcal{L}^{-1}$ using the discrete evolution operator $P_{\Delta t}$.

In view of (40) and (15), and since $\widetilde{C}^{\infty}(\mathcal{M})$ is stable under $\mathcal{L}^{-1}$, it holds

where the sum is convergent in view of (18). In conclusion,

5. Proof of Theorem 3

We start by highlighting the martingale part of the increments $\delta_{\Delta t}(q^{n},G^{n},U^{n})=Q^{n+1}-Q^{n}$, similarly to the continuous case (compare (31)) :

where the expectation on the right-hand side is over the Gaussian random variable $G$ and the uniform variable $U$ (the configuration $q^{n}$ being fixed). It will be useful to decompose the increment as

In view of Lemma 3 (which shows in particular that $\delta^{\rm reject}_{\Delta t}(q,G,U)$ can be thought of as being of order ${\Delta t}^{2}$ by setting $p=1$ in (43)), the first term on the right-hand side of (41) is equal to $\sqrt{2{\Delta t}}\,G^{n}$ at dominant order in ${\Delta t}$. This term therefore corresponds to the term $\sqrt{2}W_{t}$ in the decomposition (31). The second term in the right-hand side of (41) is handled by introducing an appropriate Poisson equation (see (44) below), which is the discrete analogue of (30).

The average increment has the following expansion in powers of ${\Delta t}$:

where $r_{\delta,{\Delta t}}$ is uniformly bounded for ${\Delta t}$ sufficiently small. In addition, there exists a constant $K>0$ such that $\left|\delta^{\rm reject}_{\Delta t}(q,G,U)\right|\leqslant K\left(1+\sqrt{{\Delta t}}|G|\right)$ and, for any $p>0$,

There exists, for any ${\Delta t}>0$, a unique function $\Phi_{\Delta t}=(\Phi_{{\Delta t},1},\dots,\Phi_{{\Delta t},dN})\in\left(\widetilde{L}^{\infty}(\mathcal{M})\right)^{dN}$ such that

Moreover, recalling the definition (30) of $\Phi_{0}$,

where $\Psi_{\Delta t}$ is uniformly bounded for ${\Delta t}$ sufficiently small, and $\widetilde{\Phi}^{1}$ is the unique solution of the Poisson equation

Since $\Phi_{\Delta t}$ is uniformly bounded as ${\Delta t}\to 0$, we obtain

We now expand $M^{0}_{\Delta t}$ in powers of ${\Delta t}$. By Lemma 4, it is possible to replace the function $\Phi_{\Delta t}$ in the second term on the right-hand side of (46) by $\Phi_{0}+{\Delta t}\widetilde{\Phi}^{1}$, up to a remainder of order ${\Delta t}^{3/2}$. We next use the following lemma to compute the cross correlation between the various functions of ${\Delta t}$ appearing in $M^{0}_{\Delta t}\otimes M^{0}_{\Delta t}$.

For any smooth functions $f,g$, growing at most polynomially,

The conclusion then follows by applying this result with the functions $f,g$ replaced by $x\mapsto F^{T}x$, $x\mapsto F^{T}\Phi_{0}(q+x)$ and $x\mapsto F^{T}\widetilde{\Phi}^{1}(q+x)$ for a given test direction $F$. Indeed,

with $|r_{F,{\Delta t}}|/|F|^{2}$ uniformly bounded as ${\Delta t}\to 0$. Manipulations similar to the ones leading to (32) finally give the claimed result.

In order to characterize the leading order term in the error and prove that the subleading order term indeed is of order ${\Delta t}^{3/2}$, it would be necessary to compute correlation terms involving components of the remainder term $\Psi_{\Delta t}$. This is not possible as such because the regularity of $\Psi_{\Delta t}$ is not established, and obtaining regularity result from (44) is difficult. An expansion of $\Phi_{\Delta t}$ up to ${\Delta t}^{2}$ terms (as $\Phi_{0}+{\Delta t}\widetilde{\Phi}^{1}+{\Delta t}^{3/2}\widetilde{\Phi}^{3/2}+{\Delta t}^{2}\widetilde{\Psi}_{\Delta t}$) is therefore needed in order not have to treat correlations involving the remainder $\Psi_{\Delta t}$. This, in turn, would require an expansion of $P_{\Delta t}$ up to remainders of order ${\Delta t}^{3}$, instead of ${\Delta t}^{5/2}$ as in (15). Although this does not pose any problem in principle, we chose not to follow this path in order to keep the arguments as simple as possible.

We conclude this section with the proofs of the technical results quoted above.

The bound (43) is a straightforward consequence of the equality

while $\left|\delta^{\rm reject}_{\Delta t}(q,G,U)\right|\leqslant\left|\sqrt{2{\Delta t}}\,G-{\Delta t}\,\nabla V(q)\right|$ immediately gives $\left|\delta^{\rm reject}_{\Delta t}(q,G,U)\right|\leqslant K\left(1+\sqrt{{\Delta t}}|G|\right)$. ∎

We introduce the normalized average increment, defined as the following periodic function :

Lemma 2 shows that $\Phi_{\Delta t}$ is well defined provided the periodic function $\overline{\delta}_{\Delta t}$ has a vanishing average with respect to $\mu$. To prove this statement, we start from

Since $q$ and $q^{\prime}$ play symmetric roles, $I_{0}$ vanishes. For $n\geqslant 0$, we obtain, by first exchanging the names of the dummy variables $q$ and $q^{\prime}$ and then using $T_{\Delta t}(q+nL,q^{\prime})=T_{\Delta t}(q,q^{\prime}-nL)$ as well as the invariance of $\mu$ by translations of the periodic cell,

Therefore, $I_{n}=-I_{-n}$, which allows to conclude that

To obtain an expansion of $\Phi_{\Delta t}$ in terms of fractional powers of ${\Delta t}$, we rely on (42), which implies that

This gives the result upon defining $\Psi_{\Delta t}={\Delta t}^{-3/2}\left(\Phi_{\Delta t}-\Phi_{0}-{\Delta t}\,\widetilde{\Phi}^{1}\right)$. ∎

A Taylor expansion with integral remainder gives

A simple computation using (42) and (43) shows that

where $r_{f,{\Delta t}}$ is uniformly bounded for ${\Delta t}$ sufficiently small. Therefore,

where the remainder $r_{f,g,{\Delta t}}$ is uniformly bounded. ∎

References