Frédéric Legoll

Tony Lelièvre

Upanshu Sharma

An adaptive parareal algorithm: application to the simulation of molecular dynamics trajectories

The aim of this article is to design parareal algorithms in the context of thermostated molecular dynamics. In its original setup, the fine and coarse propagators used in the parareal algorithm solve the same dynamics with different time-steps, with the goal of achieving accuracy in the limit of small time-step of the integrators involved. This is typically not useful in molecular dynamics, where one is interested in extremely long trajectories and where the time-step of the fine propagator is in practice chosen as large as possible, that is close to the limit of stability of the numerical scheme. In this article, we consider a version of the parareal algorithm which is better suited to molecular dynamics simulations, and wherein the propagators involved use the same time-step while employing different potential energy landscapes to drive the dynamics.

Although the parareal algorithm always converges, it suffers from various limitations in this context: intermediate blow-up of the trajectory (which makes it impossible to postprocess) may be observed; in certain cases the trajectory encounters undefined values before converging (the way the algorithm handles them might depend on the computer architecture); more critically, the computational gain of the algorithm in terms of wall-clock time (compared to a standard sequential integration) converges to one when the time-horizon increases. We highlight these issues with numerical experiments and provide some elements of theoretical analysis. We then present a modified version of the parareal algorithm wherein the algorithm adaptively divides the entire time-horizon into smaller time-slabs where the aforementioned issues are circumvented. Using numerical experiments on toy examples, we show that the adaptive algorithm overcomes the various limitations of the standard parareal algorithm, thereby allowing for significantly improved gains.

Introduction

A central objective in molecular dynamics (MD) is the computation of ensemble averages and dynamical quantities, both of which involve averages over long trajectories of stochastic dynamics. Computations using such dynamics are often difficult in practice since typical quantities of interest evolve on time-scales which are orders of magnitude larger than those accessible by classical numerical methods. Consequently, the development of numerical algorithms which speed up computations is a challenging perspective in MD. More specifically, these issues arise due to extremely high-dimensional non-convex potential energy landscapes, which trap the system in a certain region for extremely long times, before letting it explore other parts of the state space (a phenomenon called metastability). In practice, this issue may be mitigated by using approximate potentials, which are cheaper to compute (and thus make it less challenging to reach the requested long time horizons) while retaining the essential features of the original potential. However this strategy suffers from the obvious drawback that it only captures the original dynamics approximately.

In this article we present a numerical method that aims at efficiently simulating the original dynamics over long trajectories. To this end, we consider an adjusted version of the parareal algorithm, using an approximate potential as a predictor and the original potential as a corrector. We recall that the parareal algorithm has been originally proposed in [START_REF] Lions | A "parareal" in time discretization of PDEs[END_REF] to efficiently solve initial-value problems using parallelin-time computations. More precisely, the algorithm exploits the presence of multiple processors to reduce 1 the wall-clock time needed to obtain a solution on some time interval. It is based on a decomposition of the time interval into subintervals and makes use of a predictor-corrector strategy, in which the calculation of the corrections is performed concurrently, in parallel, on the different processors that are available. The algorithm has been successfully applied to many problems, see e.g. [START_REF] Farhat | Time-decomposed parallel time-integrators: theory and feasibility studies for fluid, structure, and fluid-structure applications[END_REF][START_REF] Fischer | A parareal in time semi-implicit approximation of the Navier-Stokes equations[END_REF][START_REF] Garrido | A convergent algorithm for time parallelization applied to reservoir simulation[END_REF]. Several variants of the algorithm have been proposed for specific applications: we refer e.g. to [START_REF] Blouza | Parallel in time algorithms with reduction methods for solving chemical kinetics[END_REF][START_REF] Maday | Parareal in time algorithm for kinetic systems based on model reduction[END_REF], [START_REF] Legoll | A micro-macro parareal algorithm: application to singularly perturbed ordinary differential equations[END_REF] and [START_REF] Engblom | Parallel in time simulation of multiscale stochastic chemical kinetics[END_REF] for multiscale-in-time problems, to [START_REF] Dai | Symmetric parareal algorithms for Hamiltonian systems[END_REF][START_REF] Dai | Stable parareal in time method for first-and second-order hyperbolic systems[END_REF] for Hamiltonian ODEs or PDEs, to [START_REF] Sharma | Parallelization in time of (stochastic) ordinary differential equations[END_REF][START_REF] Legoll | Parareal computation of stochastic differential equations with time-scale separation: a numerical convergence study[END_REF] for stochastic differential equations (with a specific emphasis, in [START_REF] Legoll | Parareal computation of stochastic differential equations with time-scale separation: a numerical convergence study[END_REF], to high-dimensional SDEs stemming from MD and for which qualitative features, such as positivity and norm conservation of the associated probability distribution, are to be preserved), and also to [START_REF] Farhat | Time-decomposed parallel time-integrators: theory and feasibility studies for fluid, structure, and fluid-structure applications[END_REF][START_REF] Garrido | Convergent iterative schemes for time parallelization[END_REF].

Let us now briefly introduce the classical version of the parareal algorithm (as first presented in [START_REF] Lions | A "parareal" in time discretization of PDEs[END_REF] and next reinterpreted in [START_REF] Bal | A parareal time discretization for nonlinear PDE's with application to the pricing of an American put[END_REF]). To fix the ideas, we consider the simple evolution problem given by

dx dt = f (x), x(0) = x 0 , x(t) ∈ R d , t ∈ [0, T ].
(1.1)

Assume that we have at hand two propagators for this system: F ∆t is a fine, expensive propagator which accurately approximates the exact flow of (1.1), and C ∆t is a coarse, less expensive propagator which is also less accurate. In the classical approach, F ∆t and C ∆t are often integrators of the same dynamics (1.1) over the time range ∆t with a given discretization scheme, using a small time-step for F ∆t and a large time-step for C ∆t . More precisely, F ∆t advances the system over the time ∆t by performing several time-steps of small length δt F , whereas C ∆t advances the system over the time ∆t by performing a few (or even only one) time-steps of length δt C > δt F . The parareal algorithm iteratively builds a sequence of N -tuples {x k n } 1≤n≤N (with N = T /∆t) such that, at every parareal iteration k ≥ 0, x k n is an approximation of x(n∆t). For k = 0, the initial approximation is obtained using the coarse propagator C ∆t in a sequential manner:

x k=0 n+1 = C ∆t (x k=0 n ), x k=0 n=0 = x 0 .
The subsequent parareal iterations are given by

x k+1 n+1 = C ∆t (x k+1 n ) + F ∆t (x k n ) -C ∆t (x k n ), x k+1 n=0 = x 0 .
(1.2)

The solution to (1.2) can be computed efficiently in parallel using the following procedure. Once the solution at a parareal iteration k has been computed, i.e. once {x k n } 0≤n≤N is known, the corrections F ∆t (x k n )-C ∆t (x k n ) are first computed in parallel over each subinterval [n∆t, (n + 1)∆t], 0 ≤ n ≤ N -1. These corrections are then propagated sequentially, by adding C ∆t (x k+1 n ) to the stored correction F ∆t (x k n ) -C ∆t (x k n ), which yields the solution at parareal iteration k + 1. Note that the fine propagator is never used sequentially in the algorithm. A crucial property of the parareal algorithm is that, for sufficiently large k, it always converges to the reference solution (given by the fine propagator), since it can be checked (see e.g. [BBM + 02, BM02]) that ∀k ≥ n, x k n = F n ∆t (x 0 ).

(1.3) This result is however not very useful to understand the computational gain one can expect from (1.2), that we now discuss. The computational gain of the parareal algorithm arises from the key observation that, in (1.2), the fine propagator can be computed in parallel over different processors (since it only depends on the previous iteration). Suppose that the cost of a single evaluation of F ∆t is considerably larger than the cost of propagating the system according to C ∆t over the complete time range [0, T ]. Assuming that the cost of one call to the fine-scale propagator F ∆t is c 0 , the cost of K iterations of the parareal algorithm is K c 0 . This cost is to be compared to the cost of computing the reference solution using the fine-scale propagator sequentially, which is equal to N c 0 . Assuming that the algorithm has converged after K parareal iterations, the computational gain is then N/K. This gain is always larger than one since the number of parareal iterations required to obtain convergence is smaller than or equal to N , in view of (1.3). In practice, for many applications, convergence is actually reached for a number of iterations K considerably smaller than N , thus the practical efficiency of the algorithm. A noteworthy exception is the case of hyperbolic problems, where a larger number K of iterations is often needed to reach convergence, as observed e.g. in [START_REF] Dai | Symmetric parareal algorithms for Hamiltonian systems[END_REF][START_REF] Dai | Stable parareal in time method for first-and second-order hyperbolic systems[END_REF].

In this article, we are interested in accurately and efficiently computing long trajectories in the context of MD. Specifically, we focus on the so-called Langevin dynamics (see (2.1) below), which is a popular dynamics to simulate molecular systems at fixed temperature. In applications, it is observed that the limiting timestep is not determined by the required accuracy of the dynamics, but by its stability: simulations are thus performed with a time-step at the limit of stability of the numerical method. Therefore, the classical approach to parareal algorithms as described above, with different time-steps for the fine and coarse propagators, is not viable. Indeed, a typical choice would consist in choosing the time-step of the fine propagator close to the limit of stability (because this is how a reference computation would be performed), and thus a coarse propagator based on the same numerical scheme for a larger time-step would be unstable. Instead, it is more appropriate in MD to use propagators with the same time-step but different potentials: the fine propagator solves the Langevin dynamics with the original potential V f , while the coarse propagator uses an approximate potential V c . Using different evolutions for the propagators is atypical for the parareal literature, with the exception of [BBM + 02] (see also [START_REF] Dai | Symmetric parareal algorithms for Hamiltonian systems[END_REF]Section 7]), where the authors employ propagators with different dynamics and different time-steps in the context of short-time computations in ab-initio MD. More precisely, in [BBM + 02, Sect. 5], the fine and the coarse propagators both perform abinitio MD simulations, using pseudo-potentials (for the computation of the electronic configuration) of very different complexity (thus the difference in cost between C ∆t and F ∆t , which both perform a unique time-step to advance the system over the time range ∆t). In contrast to our work here, the dynamics considered there is the Hamiltonian dynamics (i.e. the Newton equations of motion, integrated using a Verlet algorithm of time-step ∆t) rather than the Langevin dynamics, and it is integrated until a rather short time horizon: the final time corresponds to T = 30∆t. On this short time range, the parareal algorithm performs nicely, providing accurate trajectories and a very good energy preservation after only 5 parareal iterations.

The article is organised as follows. In Section 2, we present our MD setting, including two examples that are used throughout the article. Subsequently, in Section 3, we discuss various features of the parareal algorithm in the context of MD. We furthermore detail two limitations when dealing with long trajectories:

• The parareal trajectory, as a function of k, typically goes through an intermediate regime wherein the error is extremely large before converging.

• The gain N/K of the parareal algorithm converges to one as the time-horizon over which the algorithm is used gets larger. Both of these limitations are discussed with numerical examples, and obviously affect the utility of the parareal algorithm in MD. In Section 4, we present an adaptive version of the parareal algorithm, which overcomes these issues by adaptively choosing the size of the time-slab over which the trajectory is computed. On this smaller time-slab (say [0, T 1 ] ⊂ [0, T ], to fix the ideas), the parareal trajectory quickly stabilises.

Once the computation has converged on [0, T 1 ], the algorithm moves to the subsequent part of the time domain (here [ T 1 , T ]) and again adaptively chooses the size of the time-slab (say [ T 1 , T 2 ] ⊂ [ T 1 , T ]) so that convergence is quickly reached over that time range. The procedure is repeated until we reach the targeted time horizon T . The numerical results reported in Section 4 show that the gain dramatically improves compared to the standard version of parareal, which consists in running the algorithm at once on the whole time range [0, T ]. In Section 5 we present conclusions and additional discussions. We next collect in Section 6 some theoretical analysis on a toy example, which sheds light on the limitations of the standard parareal algorithm. From that theoretical viewpoint, studying the Langevin dynamics is challenging, and we therefore prefer to consider a simpler problem, namely the case of (1.1) in dimension d = 1, with forces driven by harmonic potentials. Even though simple, this case is useful for a two-fold reason. First, we observe on that case the same limitations of the classical parareal algorithm as those pointed out above and discussed in Section 3 for more realistic problems. Second, it allows to understand the interplay (on the qualitative behaviour of the algorithm along the parareal iterations k) between the time horizon T and the difference between the coarse and fine propagators. The analysis helps to quantify (on admittedly a toy-problem case) the expected fact that the algorithm can be used on longer time windows whenever the coarse propagator is closer (in a sense made precise) to the fine propagator. Stated otherwise, for a given couple of coarse and fine propagators, there is a maximum time horizon beyond which the parareal algorithm does not behave satisfactorily. Our main result in that section is Proposition 6.5.

The molecular dynamics model problem

In this article we consider the Langevin dynamics

dq t = p t dt, dp t = -∇V (q t ) dt -γ p t dt + 2 γ β -1 dW t ,
(2.1) with initial condition (q 0 , p 0 ). Here (q t , p t ) ∈ R 2md are the positions and momenta at time t of m particles in dimension d, V : R md → R is the potential energy of the system, γ > 0 is the friction coefficient, β is the inverse temperature and W t is a standard Brownian motion in R md . We have set the mass of each particle to unity for simplicity. Throughout this article, we use the following numerical scheme (see [JS12, Eq.(4.2)]) to integrate the Langevin dynamics:

p n+ 1 4 = p n - ∆t 2 ∇V (q n ), q n+1 = q n + ∆t p n+ 1 4 , p n+ 1 2 = p n+ 1 4 - ∆t 2 ∇V (q n+1 ), p n+1 = e -γ∆t p n+ 1 2 + (1 -e -2γ∆t ) β -1 G n , (2.2)
where G n is a standard Gaussian random variable. This scheme uses a Verlet integration for the Hamiltonian part of (2.1) (which corresponds to the first term in the right hand side of both lines in (2.1)) and an explicit solution for the Ornstein-Uhlenbeck part of (2.1) (which corresponds to the last two terms in the right-hand side of the second line of (2.1)).

As stated in the introduction, in the MD context, it is more appropriate to use propagators with the same time-step but different potentials V , and thereby the propagators solve different dynamics. Henceforth, the fine (resp. coarse) propagator F ∆t (resp. C ∆t ) corresponds to (2.2) with the potential V f (resp. V c ). Note that F ∆t and C ∆t amount to advancing the system for a time range ∆t by performing a single time-step of the integrator (2.2). Of course, for the parareal algorithm to converge, it is crucial that the coarse and fine dynamics are driven by the same noise. We thus use Gaussian increments G n in (2.2) which, at a given time-step n, are the same for both propagators F ∆t and C ∆t and throughout all parareal iterations. Remark 2.1 (Overdamped Langevin dynamics). While here we focus on the Langevin dynamics since it is preferred in practice, other choices are also possible. For instance, one can use the overdamped Langevin dynamics

dq t = -∇V (q t ) dt + 2 β -1 dW t , (2.3)
which arises as the large-friction limit of (2.1) (see [LRS10, Section 2.2.4]). It turns out that the difficulties we encountered when applying the parareal algorithm described in this article to the Langevin dynamics also appear when working with its overdamped limit (2.3).

Two examples

We now present two pairs of fine and coarse potentials that are used throughout this article.

One-dimensional harmonic potentials (Har-1d)

In the first example, we consider V f and V c to be two different harmonic potentials in dimension d = 1, i.e., for some ω > 0 and for any q ∈ R,

V f (q) = q 2 2 , V c (q) = ω q 2 2 .
(2.4) We consider this toy example since it has sufficient complexity to exhibit, in the parareal context, the same issues as those appearing with more complex choices while also being amenable to theoretical analysis (see Section 6 below).

2.1.2 Two-dimensional Lennard-Jones 7-atom cluster (LJ7-2d)

For a more realistic model, we consider a 7-atom (m = 7) cluster in dimension d = 2 with pairwise Lennard-Jones interaction. The energy landscape has four local minima which correspond to four stable configurations (see Figure 1). The fine potential V f : R 14 → R is

V f (q) = 1 2 i,j∈{1,...,7} i =j φ f (|q i -q j |), φ f (r) = r -12 -2 r -6 ,
where | • | is the 2-dimensional Euclidean norm. The coarse potential V c is the harmonic approximation of V f around the initial well, i.e. the well in which the initial condition lies (see Appendix A for explicit formulae).

In practice, we choose the initial position to be the global minimizer of the potential energy, that is the deepest well (see Figure 1a). The choice of initial momenta is discussed in the following Section 2.2.1.

Implementation issues

We now briefly discuss some implementation issues.

Initial conditions

We first discuss how initial conditions are chosen. The invariant measure for the Langevin dynamics is the Boltzmann-Gibbs measure with density

µ(q, p) = Z -1 exp -β |p| 2 2 + V f (q) ,
where | • | is the Euclidean norm in R md and where Z is the normalisation constant which ensures that µ is a probability density. In the Har-1d case, since V f (q) = q 2 /2, we simply draw the initial condition q 0 and p 0 independently according to N (0, β -1 ). Throughout this article, for Har-1d, we work with β = 3, and the corresponding initial condition is (q 0 , p 0 ) = (-0.386593, 0.48453).

For the LJ7-2d, we cannot use the same strategy as above to choose an initial condition, since it is in practice very challenging to draw q 0 according to the correct probability measure. Instead we follow a burnin procedure which we now describe. We start the particles with random initial momenta (drawn from the invariant measure N (0, β -1 Id) for the momenta, where Id is the identity matrix of size md × md) and initial positions at the bottom of the deepest well corresponding to Figure 1a (six atoms located π/3 radians apart on a unit circle and one in the centre), and we next run the Langevin dynamics. We integrate this dynamics until the time N burn-in ∆t with N burn-in chosen such that the average kinetic energy (averaged over time, dimension and number of atoms)

1 14N burn-in N burn-in n=1 i∈{1,...,7} 1 2 p i (n∆t)
2 is sufficiently close to its target value 1/(2β). This configuration, obtained at time N burn-in ∆t, is then chosen as the initial condition for the parareal procedure.

Translation-invariance for LJ7-2d

We impose translation invariance, since it ensures that translating the entire system does not contribute to the errors. This is implemented by using zero-mean noise in both directions. More precisely, for each particle i ∈ {1, . . . , 7} and at each time-step, we draw a Gaussian random variable G i x for the x-coordinate, and G i y for the y-coordinate (recall (2.2)). The idea is then to replace G i x by G i x -

1 7 7 j=1
G j x , and likewise G i y by

G i y - 1 7 7 j=1 G j y . Using that V f satisfies 7 j=1 ∂V f ∂q j x = 7 j=1
∂V f ∂q j y = 0, we then deduce from (2.2) that, if

7 j=1 p j n = 0, then 7 j=1
p j n+1 = 0. The center of mass of the system hence does not move. We note that this scheme actually corresponds to integrating the Langevin dynamics projected on the manifold (q, p) ∈ R 14 , 7 j=1 p j = 0 . We underline that the conclusions drawn on that case also hold in the absence of translation invariance, that is when we directly implement (2.2) without modification of the Gaussian increments. Since it is commonly employed in MD, we have decided here to implement this translation-invariance feature.

Parareal within a single well of the potential energy

In all the experiments on LJ7-2d, we work with parameters (time horizon T , friction coefficient γ, inverse temperature β, . . . ) such that the fine trajectory remains within the deepest well in which it initially lies. This requirement is imposed keeping in mind that using a harmonic approximation around the initial well is not a reasonable choice when the trajectory is close to a saddle point or leaves the well. The question of choosing a good coarse potential is a challenging one and often requires an understanding of the system being studied. In this article we are interested in the general features of the parareal algorithm, and consequently the study of better coarse potentials on the entire state space is left to future work (see Section 5.2 below for additional details).

Choice of the friction coefficient

One possible way to choose the friction coefficient γ in the Langevin equation (2.1) is by using the synchronisation rate. Considering two initial conditions which are close one to each other, we write the two dynamics (2.1) starting from these two initial conditions. Assuming that we can linearize the potential, we take the difference between both dynamics (the Brownian motions thus cancelling each other) and we look for γ such that the difference (both in position and momentum) converges to zero with the largest possible rate. This is the case when γ is set at the so-called synchronisation rate, which is equal to twice the square-root of the smallest eigenvalue of the Hessian of the potential at the bottom of the well.

In the Har-1d case, where V f (q) = q 2 /2, the synchronisation rate is thus equal to 2. In the LJ7-2d case, we ignore the first three eigenvalues of the Hessian, which correspond to translation and rotation in dimension d = 2, and are very close to zero. The first non-zero eigenvalue for this setup is 69.1, which gives a synchronisation rate of 16.62.

We will not systematically use exactly these values, but we will consider values of γ of this order of magnitude for each system. We will observe below that the efficiency of the parareal algorithm significantly depends on γ.

Numerical results using the classical parareal algorithm

We now introduce some quantities useful to study the parareal algorithm. For k ∈ {1, . . . , N }, the relative error (for the entire trajectory) between two consecutive parareal trajectories is defined as

E(k, N, γ, β) = N n=1 |q k n -q k-1 n | N n=1 |q k-1 n | . (3.1)
Since the parareal algorithm always converges, for a fixed parameter δ conv > 0 (convergence criterion), we stop the algorithm at the first parareal iteration k for which E(k, N, γ, β) < δ conv . For a given δ conv , the gain of the parareal algorithm is then defined as

gain(N, γ, β) = N k . (3.2)
Note that this definition of gain assumes in particular that the cost of the coarse propagator is negligible in comparison with the cost of the fine propagator. This definition is also centered around the notion of wall-clock time, in contrast to the notion of total CPU cost. As soon as several parareal iterations are used, the total CPU cost of the parareal algorithm is larger than the cost of a single, sequential computation of the trajectory using the fine propagator. But the wall-clock time may be smaller, and this is precisely quantified by (3.2). For the sake of simplicity we choose

δ conv = 10 -5 (3.3)
for all the numerical experiments in this article. The relative error (for the entire trajectory) between the fine trajectory (which, we recall, is denoted (q(t), p(t)) for 0 ≤ t ≤ T = N ∆t) and the parareal trajectory is

E f (k, N, γ, β) = N n=1 |q(n∆t) -q k n | N n=1 |q(n∆t)| . (3.4)
Since the parareal algorithm always converges, we have

E(N + 1, N, γ, β) = 0 = E f (N, N, γ, β).
Remark 3.1 (Error at terminal time). Note that the relative error E defined by (3.1) is quite strong since it requires convergence over the entire trajectory. Another (weaker) measure of error is

|q k N -q k-1 N | |q k-1 N |
, which keeps track of the relative error between two consecutive parareal trajectories only at the terminal time-step.

On the examples we have considered, it turns out that this error scales the same way as E. Consequently, in this article, we use the relative errors for the entire trajectory (3.1) or (3.4) as a measure of error. We comment in Section 5.1 below on other measures of error, well adapted to the MD context.

The standard parareal algorithm exhibits two crucial features when computing long trajectories (i.e., when T = N ∆t is large) which we now discuss.

Intermediate explosion of error for large times

For small N , the parareal trajectory converges to the reference trajectory fairly quickly when k increases. This is illustrated in Figure 2, where we plot the relative error E f between the parareal and the fine trajectory as a function of k for various values of γ. For k of the order of 10, we have a converged trajectory over a time horizon N ∆t with N of the order of 100 (thus a gain of the order of 10). When γ increases, this convergence improves as well. This is explained by the fact that γ is the friction parameter in the Langevin dynamics (2.1), and therefore increasing γ leads to the system becoming increasingly dissipative and thus less movement of the particles. In addition, for very large values of γ, the physical forces term in the second line of (2.1) is dominated by the next two terms, and thus the coarse and the fine propagators integrate dynamics which are close one to each other. We now comment upon the non-monotone profile of E f for γ = 15 that we observe in Figure 2b. As we shall now see, this is a typical feature of the parareal algorithm. When T = N ∆t is large, the parareal trajectory has an intermediate non-monotone error profile with respect to the fine trajectory: the error on [0, T ] typically increases to very large values when k increases, before eventually decreasing and converging to zero. This is obvious in Figure 3, where we have run the same experiments as in Figure 2 (with essentially the same values of γ) but for a much larger value of N . We observe that the error E f reaches extremely large values for intermediate values of k. In fact, in some cases, the parareal trajectory has undefined (NaN) values for certain parareal iterations, before converging when k becomes closer to N . For instance, the curve for γ = 30 in Figure 3b corresponds to such a situation: it is not complete and reaches NaN values at k ≈ 240. This phenomenon is explained by the fact that the parareal trajectory goes far away from the fine trajectory (and in fact takes extremely large values) before eventually converging. This is illustrated in Figure 4 for Har-1d (where we have set N = 3000 as in Figure 3a), where, for an intermediate value of k, the parareal trajectory has values that are orders of magnitude larger than the fine, reference trajectory. We work here with γ = 3, which corresponds to the purple curve in Figure 3a. We observe there that the error is maximal for k ≈ 40. On the top right of Figure 4, we show the parareal trajectory for k = 40: for n ≥ 1400, it takes values that are extremely large compared to the reference values, thus the very large error. When we perform additional parareal iterations and consider k = 80, the error is much smaller, as can be seen in Figure 3a, although still large. The parareal trajectory for k = 80 (see bottom left of Figure 4) now blows up only for n ≥ 2700, and is a least of the correct order of magnitude (if not accurate) for n ≤ 2700. We can expect from Figure 3a that the error is small for k = 150, a fact confirmed by the plot on the bottom right of Figure 4: the parareal trajectory for k = 150 can essentially be superimposed to the reference trajectory shown on the top left of Figure 4. Remark 3.2. In Figures 2, 3 and4, we have chosen to work with ω = 0.1 for the coarse potential in the Har-1d case (which is then a bad approximation of the fine potential) to mimick the fact that, for realistic examples, it may be difficult to design coarse potentials V c which are at the same time less expensive to compute than V f and sufficiently close to V f . On purpose, we have chosen here a bad coarse potential.

The fact that the parareal trajectory goes extremely far away for intermediate values of k is not the sole reason for slow convergence of the iterations. Actually, for large T , the parareal algorithm converges slowly even when stable. This is illustrated in Figure 5, where we use the following "trust-region based" fine and coarse potentials:

V trust f (q) =    q 2 2 if q ∈ [-3, 3] 0 otherwise , V trust c (q) =    ω q 2 2 if q ∈ [-3, 3] 0 otherwise . (3.5)
We consider the same numerical parameters as in Figure 4. We see on the top left plot of that figure that the fine trajectory always satisfies |q(t)| ≤ 3. Changing the fine potential from V f to V trust f hence does not affect the fine, reference trajectory. By using these potentials V trust f and V trust c , we make sure that the parareal trajectory experiences zero force when leaving the trust-region. As a result, the parareal trajectory (and consequently the error) stays bounded (see Figure 5). However, the parareal trajectory using these modified potentials converges to the reference trajectory as slowly as when using the original propagators. This substantiates the claim that the parareal algorithm, even when stable, converges slowly for large time horizon T > 0. This is consistent with the theoretical analysis for harmonic potentials presented in Section 6 (see in particular Figures 12a and12b, in which the error remains bounded but may converge extremely slowly to zero).

The intermediate explosion discussed above leads to various practical issues. In MD, one is typically interested in certain features of the trajectory that need to be calculated online, such as the current basin of attraction. However, these quantities cannot be computed (or are meaningless) for such parareal trajectories due to the unphysically large values. Furthermore, the NaN values pose obvious issues since they are architecture dependent.

In Figure 6, we plot E (the relative error between two consecutive iterations) and E f (the relative error with respect to the reference trajectory) as a function of k, on a simulation performed with again the same numerical parameters as in Figure 3. Recall that the parareal algorithm only uses E, and not E f , as a n , p 80 n ) in log-scale for k = 80 (the gaps in the graph correspond to negative values). Bottom right: plot of (q 150 n , p 150 n ) for k = 150. While the trajectory is stable at larger k, it blows up for intermediate k. stopping criterion. However, as clearly seen on the figure, the explosion of the parareal trajectory and therefore E f can be correlated to an increase in E beyond a certain threshold. This observation plays a key role in the adaptive parareal algorithm introduced in Section 4.

Gain converges to one for large N

We now focus our investigations on the gain provided by the algorithm, which is defined by (3.2). We numerically observe that the gain converges to one when N (or equivalently T ) becomes larger, and consequently the parareal algorithm does not offer any advantage, for long time simulations, over using the fine integrator in a sequential manner. This is illustrated in Figure 7, where we plot the gain for both examples as a function of N for different values of γ. Note that, while increasing γ improves the gain at fixed N (recall that we have already observed in Figure 2 that the error is smaller when considering larger values of γ), the gain nevertheless converges to one when N → ∞, whatever the value of γ. An expected feature of the parareal algorithm is that the gain decreases as the coarse potential becomes a worse approximation of the fine integrator. This is illustrated in Figure 8, where we plot the gain for Har-1d case, as a function of ω, for fixed N , γ and β. As expected, the gain decreases as ω takes values further away from the reference value ω ref = 1. A similar conclusion also follows from the theoretical analysis presented in Section 6. Comparing the left and the right plot in Figure 8, we also observe that the algorithm behaves differently whether ω > ω ref or ω < ω ref , a fact again consistent with the analysis of Section 6, where the sign of y defined by (6.4) (or equivalently the sign of ω -ω ref ) matters. In Table 1 (given at the end of Section 4 below), we collect in blue the gains for LJ7-2d for various γ and β. Note that, for smaller β, the gain is rather insensitive to the value of β, whereas it considerably increases for β ≥ 18. This is due to the fact that for larger β, the effect of the noise becomes less pronounced in the Langevin dynamics and the system tends to stay close to the bottom of the well, where the harmonic approximation of the potential is an excellent approximation. In contrast, for smaller β, the noise has a considerable effect on the trajectory, which explores a much larger part of the state space.

Adaptive algorithm

Let us summarise the two important observations from Section 3:

• The relative error between the fine and the parareal trajectory decreases when T decreases. We have observed this behaviour in all the above numerical experiments. This behaviour is also consistent with the theoretical analysis presented in Section 6.

• If the relative error between consecutive parareal trajectories remains larger than some threshold (say larger than 1), then the gain of the parareal algorithm is very limited (see in particular the discussion at the end of Section 3.1, where values of E larger than 1 are correlated with blow-up of the error).

Based on these observations, we now introduce an adaptive parareal algorithm.

Starting from the time-slab [0, T ], the idea is to run the parareal algorithm until the relative error between consecutive trajectories is either larger than an explosion threshold δ expl > 0 (we denote k cur the parareal iteration at which this is attained), or smaller than the convergence threshold δ conv (which, we recall, has been fixed at δ conv = 10 -5 , see (3.3)). In the latter case, the method has converged on the time range [0, T ]. In the former case, for the parareal iteration k cur , we find the first time iteration N final ≤ N for which the relative error exceeds δ expl , and change the slab to [0, (N final -1)∆t]. We then proceed with the parareal iterations on the smaller time-slab [0, (N final -1)∆t] (note that we do not need to perform again the first k cur iterations, we can just proceed from the knowledge of {x kcur n } 0≤n≤N final -1 ), until the relative error (now of course on [0, (N final -1)∆t]) between consecutive trajectories is either larger than δ expl , or smaller than δ conv . In the latter case, the method has converged on the time range [0, (N final -1)∆t]. In the former case, we again reduce the size of the time-slab. We hence continue to find smaller and smaller slabs until the parareal algorithm converges, and we denote [0, T 1 ] = [0, m 1 ∆t] the slab which has eventually been found.

Once we have a converged parareal trajectory on [0, m 1 ∆t], we proceed with the next part of the time range and define the new time-slab as [ m 1 ∆t, T ]. We proceed with the parareal iterations (starting from the converged configuration at time m 1 ∆t), and possibly reduce the size of the time-slab to [ m 1 ∆t, m 2 ∆t] in order to reach convergence.

This procedure is repeated until we have reached the time horizon T .

The pseudocode for this algorithm is presented in Algorithm 1, where we use the following notation for the relative error over a given time interval: for N init , N final ∈ N with N init ≤ N final and two trajectories {x prev n } Ninit≤n≤N final and {x cur n } Ninit≤n≤N final , the relative error is

E(N init , N final ; x prev , x cur ) = N final n=Ninit |x cur n -x prev n | N final n=Ninit |x prev n | .
The relative error between two consecutive parareal trajectories is computed in line 11 of the algorithm. We actually do not wait to reach the current end point of the time-slab to compute it, but rather compute it for each new time point in [N init ∆t, N final ∆t], and exit the computation (to redefine a new endpoint) whenever this relative error is larger than δ expl . Note that it is of course advantageous to compute independently the numerator and the denominator of this relative error (which are easy to increment for each new value of

x cur n+1 , with N init + 1 ≤ n + 1 ≤ N final ).
The while instruction on line 5 corresponds to the parareal iterations: as long as δ > δ conv (no convergence yet) and δ < δ expl (no explosion), we proceed with the iterations. Note that we always exit this while loop, since we know that, for sufficiently large k, the parareal algorithm converges (recall (1.3)): if the trajectory does not blow up (i.e. if δ remains smaller than δ expl ), then the algorithm reaches some iteration where δ < δ conv . The else condition on line 20 corresponds to the case δ < δ conv (since we have exited the while loop, we can only have δ > δ expl or δ < δ conv ).

The key parameter in the adaptive algorithm is the threshold of explosion, which we have denoted by δ expl . More precisely, the algorithm adaptively chooses the slab size for which the relative error is smaller than δ expl . This means that for smaller δ expl , these slabs are of smaller size and we need a larger number of slabs to reach the final time horizon T . In contrast, if δ expl were chosen very large, the relative error would always be smaller than δ expl and a unique slab, namely [0, T ], would be used. This is illustrated in Figure 9 for the Har-1d case for two large values of N . The influence of δ expl on the gain is discussed below. We simply point out for now that, in this particular case and for N = 5000 (see Figure 10a and also, for larger values of N , Figure 11a), the optimal choice for best gain is δ expl = 5. Figure 9 shows that the algorithm builds 4 slabs (resp. 20 slabs) to converge on [0, N ∆t] with N = 5000 (resp. N = 30000), and the size of these slabs varies from one to the next. The number of CPUs required in practice to obtain a perfect parallelization of the fine propagators thus corresponds to the largest slab-size. For instance, the case N = 30000 and δ expl = 1 (resp. δ expl = 5) requires close to 440 (resp. 2700) CPUs. While this may be considered as a large number (and while some of these CPUs are not used when the slab size is shorter), it is considerably less compared to the classical parareal algorithm which would require N = 30000 CPUs. Furthermore the number of CPUs required is fairly robust when N increases (compare the left and the right plots of Figure 9, e.g. for δ conv = 1 and for δ conv = 3). We also observe that the fluctuation of the size of the slabs seems to depend on δ expl , and seems to be larger when δ expl is larger.

We define the cost (as above, in terms of wall clock time and neglecting the cost of the coarse propagator compared to the fine propagator) for the adaptive parareal algorithm as follows. Consider a tentative timeslab [ m j ∆t, T ], which has possibly been eventually shortened to [ m j ∆t, m j+1 ∆t] in order for the parareal iterations to reach convergence on that interval. Since the fine propagations are performed in parallel, the cost of the algorithm over the time-slab [ m j ∆t, m j+1 ∆t] is the total number C j of parareal iterations that have been performed for the trajectory to converge on [ m j ∆t, m j+1 ∆t] (or equivalently the number of times the fine integrator has been called to propagate the system from m j ∆t to (1 + m j )∆t; of course, on the time-slab [ m j ∆t, m j+1 ∆t], the fine propagator has to be called to propagate the system further in time than (1 + m j )∆t, but this is performed in parallel and thus should not be accounted for in terms of wall-clock time). We then define the cost of the adaptive variant as the sum of these C j over all the time-slabs that have been needed to reach the final time horizon T : if the time interval [0, T ] has been decomposed by the adaptive algorithm as [0, T ] = ∪ 0≤j≤J-1 [ m j ∆t, m j+1 ∆t], then the cost is defined as

0≤j≤J-1 C j . The reason
for this definition is that the computation from one time-slab to the next is performed in a sequential manner (the trajectory should be converged on [ m j ∆t, m j+1 ∆t] before we move on to [ m j+1 ∆t, m j+2 ∆t]), and thus the wall clock times should add up. The computation of this cost is performed in line 8 of the algorithm. Following the definition (3.2) of the gain for the classical parareal algorithm, we define the gain of the adaptive parareal algorithm (compared to a sequential computation) as the ratio of N divided by the cost defined above. Note that, if the adaptive variant turns out to use a unique slab to reach the final time horizon, then both definitions of gain, (3.2) and the present one, coincide. In Figure 10 For instance, when N = 5000 and δ expl = 1, the total time interval [0, T ] is eventually divided into 31 slabs, each of which of size of the order of 160 time-steps. The 12-th slab is the longest and its length is of the order of 350 time-steps. The 31st slab, which is the last one needed to reach the final horizon, is the shortest with a length of the order of 50 time-steps.

When N = 5000 and δ expl = 5, the total time interval [0, T ] is eventually divided into 4 slabs, of length close to 1300, 1500, 1200 and 1000 time-steps.

for the LJ7-2d case), the gain seems to be fairly robust and independent of N . Comparing the left and the right plots of Figure 10, we note that the gain increases when γ increases. As expected, for the largest values of δ expl , the adaptive parareal algorithm is identical to the classical parareal algorithm. It thus leads to the same gain, which converges to 1 when N increases, as already discussed in Section 3.2. For large values of N (which is our regime of interest here), the adaptive algorithm always outperforms the classical version. For a good choice of δ expl , we can reach a gain of the order of 30 for the Har-1d case, and of 7 for the LJ7-2d case. In the LJ7-2d case, and in contrast to the Har-1d case, we also observe that the adaptive variant outperforms the classical variant even for small values of N .

In Figure 11, we plot the gain as a function of δ expl , for fixed N , γ and β. Note that the gain is fairly stable for long-times: we have already observed this in Figure 10, and we confirm this observation by comparing the two curves in the Har-1d case in Figure 11a. The large plateau in the LJ7-2d case is due to the fact that the relative error E has extremely high spikes (recall Figure 6c): whether we set δ expl = 0.1 or δ expl = 100, adaptivity of the time-slabs is triggered at the same moment, corresponding to the high spikes. The comparison between Figures 11b and11c allows to investigate the influence of γ. We observe that, when γ increases, the plateau increases in size as well, which is in agreement with the fact that the algorithm is more stable when γ increases. In all the plots of Figure 11, the gain is maximal for a certain optimal choice of δ expl , which we expect to be system-dependent. When δ expl becomes larger than the optimal choice, the adaptive algorithm allows for larger slabs and thereby tends to behave like the classical parareal algorithm, which explains the diminishing gain. In contrast, when δ expl is much smaller than the optimal choice, the time-slabs become very short, and there is almost no parallelisation over different subintervals anymore, thus again a diminishing gain.

We collect in Table 1 the gains for LJ7-2d for various γ and β, for the classical parareal algorithm (in blue) and for the adaptive variant (in black). We observe that, when the classical variant does not perform satisfactorily (and provides a gain close to 1), then the adaptive variant significantly outperforms it. We have not observed gains smaller than 4 with the adaptive variant. When the classical variant already performs in an efficient manner, then the adaptive variant provides similar (if not larger) gains. For smaller values of 1: This table collects the gain(N, γ, β) for LJ7-2d with ∆t = 0.01 and N = 1000 for various γ and β. The values in blue correspond to the classical parareal algorithm and those in black to the adaptive parareal algorithm with δ expl = 1. Note that, for each different value of β, a different initial condition has been chosen, using the burn-in procedure described in Section 2.2.1.

β, the adaptive variant performances are slightly increasing when β increases. They dramatically improve for larger values of β. In contrast to the classical variant, the evolution of the performances of the adaptive variant as a function of β is smoother. In addition, when γ increases, the gain of the adaptive variant steadily increases, in contrast to the classical variant for which the gain increase may be very limited.

Discussions

In this article, we have considered a version of the parareal algorithm using integrators based on different dynamics and we have outlined its limitations, which include possible intermediate explosion and poor gain for longer trajectories. We have then introduced a modified version of the parareal algorithm which adaptively chooses the size of the time-slabs. The adaptive variant does not show the above mentioned shortcomings. We have observed that it systematically outperforms the classical version of the algorithm, and that it always provides significant gains on the problems we have studied.

Some comments on alternative ideas

We comment here on some related ideas from the literature about the parareal algorithm in the MD context and some associated numerical tests that we have performed.

Stabilising the parareal algorithm with projections. The authors of [START_REF] Dai | Symmetric parareal algorithms for Hamiltonian systems[END_REF] consider Hamiltonian systems (which have traditionally been challenging test cases for the parareal algorithm) and introduce symmetrised variants of the parareal algorithm (we also refer to [START_REF] Dai | Stable parareal in time method for first-and second-order hyperbolic systems[END_REF] for a study of the wave equation). The additional idea in [START_REF] Dai | Symmetric parareal algorithms for Hamiltonian systems[END_REF] is to perform a projection within the parareal iteration, such that the energy (i.e. the Hamiltonian) of the parareal trajectory matches with that of the fine trajectory. Such a projection is of course motivated by the fact that, in the Hamiltonian setting, the energy is preserved along the exact dynamics. In principle, similar ideas can be applied to the Langevin dynamics, wherein we could ensure that the parareal trajectory matches the value of either the Hamiltonian (physical energy) or of some Lyapunov function (see e.g. [Tal02, Eq. (2.1)]) along the fine trajectory. Recall in particular that such Lyapunov functions (which are built by considering the Langevin dynamics without noise) are very useful functions to explain convergence to equilibrium of the Langevin equation. Furthermore, in contrast to the physical energy, they depend on γ, and may thus appear to be a more adapted choice to a given Langevin dynamics than the Hamiltonian energy. It turns out that these matching algorithms lead to extremely slow convergence (if not lack of convergence at all). The reason is that the Langevin dynamics does not have a purely dissipative structure (which would be relevant for Lyapunov matching), neither of course a Hamiltonian structure. γ-adaptive parareal algorithm. Recently, the authors of [START_REF] Maday | An adaptive parareal algorithm[END_REF] consider a variant of the classical parareal algorithm, where the time-step of the fine integrator is adaptively decreased along the parareal iterations. The heuristics behing this strategy is that it is useless, for the first few parareal iterations, to accurately compute the fine propagator, since these first iterations only serve as an "initial guess" for the latter ones. This strategy is shown to lead to better gains. While decreasing the step-size is not a feasible strategy in the MD context, a similar idea could be implemented by choosing a fine integrator which starts at a larger value of γ (since parareal is more stable at larger γ) and then adaptively reduces γ to the original value along the parareal iterations. While this idea offers positive results when the original dynamics starts with large γ, it does not even converge for intermediate values of γ such as the ones discussed in this article. A typical result is that the first parareal iterations are somewhat stable. Then, when γ decreases below a certain threshold (irrespective of how slowly γ decreases along the parareal iterations), instabilities appear and prevent the algorithm from converging.

Trust-region based parareal algorithm. Another possible strategy, inspired by trust-region optimisation, is to introduce a modified coarse propagator which uses the coarse potential when inside the trust-region and the fine potential otherwise (it is thus a different -although related -idea than (3.5)). A motivation is that instabilities appear when the trajectory leaves the domain where it should be, and that reverting to the fine potential in that case may help fix the instabilities. It turns out that this strategy is costly in terms of wall-clock time, since the expensive fine propagator happens to be evaluated in a sequential manner, which defeats the purpose of the parareal algorithm in the first place. In fact, even when a zero potential (instead of V f ) is used for the coarse-propagator outside the trust-region, the gain behaves as in the classical parareal case (see Equation (3.5), Figure 5 and the corresponding discussion).

Trajectorial versus statistical properties. In the present context of MD, we recall that statistical properties of the system (as opposed to individual trajectories, which is the focus of this article in view of the definition (3.4) of the error) are of interest. They are typically less demanding in terms of accuracy for the numerical scheme: examples abound for which the numerical trajectory is plainly wrong but its statistical properties are very accurate. However, as we have observed, the individual trajectories provided by the classical parareal algorithm are far away from the fine trajectories at long times. They thus cannot be expected to yield correct averages. From the perspective of statistical quantities, the classical variant thus presents the same limitations as for trajectory computations.

Open questions

We now summarise some open problems and describe some perspectives of this work.

Why does the parareal trajectory explode before converging? We have shown that, in several MD cases, the parareal trajectory blows-up before eventually converging to the fine trajectory. It is yet unclear to us why exactly such an explosion occurs. Note that this feature is not observed in the (classical) setting when the fine and coarse propagators solve the same parabolic (i.e. dissipative) dynamics with different time-steps. We have observed it for parabolic problems when the two propagators are based on different dynamics (see Section 6 below, in particular Figure 12). This issue has also been observed when solving hyperbolic problems (see e.g. [DLBLM13, Figs. 4 and5]).

Global coarse potential for LJ7-2d. In this article, we have used the harmonic approximation of the fine potential around the deepest well as the coarse potential, and we have considered initial conditions in that deepest well. However, such an approximation is a bad choice if one is interested in conformational changes which involve a trajectory visiting multiple wells. A natural next step of this work is thus to explore global coarse potentials and ascertain the gain achieved by the (classical and adaptive) parareal algorithms in this setting. We intend to explore this question in future works on more realistic physical systems, using various force fields.

Focus on dynamical quantities. Quantities of interest in MD include thermodynamic quantities and dynamical quantities. For the computation of the former, several methods have been proposed in the literature (see e.g. [START_REF] Lelièvre | Free Energy Computations[END_REF]). Some of them rely on the long-time integration of the Langevin dynamics (and then a time-parallelisation tool such as the adaptive strategy introduced in this article may be very useful). But many other strategies, based on different ideas, have been proposed. In contrast, dynamical quantities (such as velocity autocorrelation plots, diffusion coefficients of defects -vacancies, interstitials, . . . -in a material, exit times from a given well, to name but a few) are defined on the basis of the dynamics. In this context, it is of paramount interest to be able to speed-up the computation of long time trajectories, for instance using the adaptive strategy introduced here. Investigating the efficiency of parareal algorithms for the computation of such dynamical quantities is a very interesting question, and we intend to explore it on more realistic systems in forthcoming works.

Theoretical analysis on a toy problem

In the above sections, we have studied parareal algorithms for the Langevin dynamics. These are not straightforwardly amenable to numerical analysis, because the variable in a Langevin dynamics is necessarily multi-dimensional. In this section, we analyse the error between the fine (reference) and the parareal trajectories, for a simple ODE in dimension d = 1 driven by harmonic potentials (we are going to actually consider the same potentials V f and V c as in the Har-1d case, see (2.4)). As we shall see, this simpler setting offers explicit expressions of the propagators and of the solution x k n for any n and k, which allows us to theoretically discuss various aspects of the classical parareal algorithm that we have outlined in Section 3 on the basis of numerical simulations. Our main result is Proposition 6.5, which describes the behaviour of the relative error as a function of n and k, depending on the quality of the coarse propagator (namely, how close the coarse and the fine propagators are).

Let us briefly mention that several theoretical analysis of parareal algorithms have already been proposed in the literature. We refer e.g. to [START_REF] Maday | A parareal in time procedure for the control of partial differential equations[END_REF][START_REF] Maday | The parareal in time iterative solver: a further direction to parallel implementation[END_REF], to [START_REF]On the convergence and the stability of the parareal algorithm to solve partial differential equations[END_REF][START_REF] Staff | Stability of the parareal algorithm[END_REF] for stability investigation, to [START_REF] Gander | Nonlinear convergence analysis for the parareal algorithm[END_REF] for a numerical analysis in a nonlinear context and to [START_REF] Gander | Analysis for parareal algorithms applied to Hamiltonian differential equations[END_REF] for an analysis in the specific Hamiltonian context. To the best of our knowledge, our result below is the first one to shed light on the interplay on the error between the quality of the coarse propagator (with respect to the fine one) and the length of the time interval on which the parareal algorithm is used. Remark 6.1. Our analysis below applies to the setting of this article, namely when the coarse and fine propagators integrate different dynamics (governed by different potentials). But it also applies to the classical case when the coarse and fine propagators integrate the same dynamics with different time-steps and/or different numerical schemes. All what matters (see e.g. our main result, Proposition 6.5) is the value of y defined in (6.4) below.

Consider the reference ODE

dx dt = -V f (x) = -x, x(0) = x 0 , (6.1)
for x : [0, T ] → R, where T > 0 is the finite-time horizon over which we study the dynamics and where V f is the reference potential defined by (2.4). Without loss of generality, we assume that the initial condition x 0 is positive: x 0 > 0. As before, we divide the interval [0, T ] into subintervals of length ∆t with N = T /∆t. The fine propagator F ∆t corresponding to (6.1) is linear and given by F ∆t (x) = F x, where F = exp(-∆t).

For ω > 0, we define the coarse propagator as the propagator (over the time length ∆t) of the coarse model dx t dt = -V c (x t ) = -ω x t (where V c is the coarse potential also defined by (2.4)). We therefore have

C ∆t (x) = Cx, where C = exp(-ω∆t).
Denoting by x k n the solution at the k-th parareal iteration and n-th time-step, the parareal algorithm reads

x k+1 n+1 = C ∆t (x k+1 n ) + F ∆t (x k n ) -C ∆t (x k n ) = Cx k+1 n + (F -C)x k n , x k+1 0 = x 0 (6.2)
and the 0-th parareal iteration is given by

x 0 n+1 = C ∆t (x 0 n ) = Cx 0 n , x 0 0 = x 0 . (6.3)
In what follows, we fix n, and study the relative error at the terminal point between the parareal and the fine trajectory

R n,k = x(n∆t) -x k n x(n∆t) = x n n -x k n x n n
where k ∈ {0, . . . , n -1}, where we have used the parareal algorithm property that x n n = x(n∆t). Note that R n,k is mathematically more tractable than the relative error for the entire trajectory that we have used to present the numerical experiments in the above sections (recall Remark 3.1). A critical parameter in this study is

y = F C -1, (6.4)
which quantifies how much the coarse and the fine propagators are far apart. Using the explicit expressions of F and C, we have y = exp (ω -1)∆t -1, which satisfies

y > 0 if ω > 1 and y ∈ (-1, 0) if ω ∈ (0, 1).
Of course, if the coarse and fine propagators are identical, namely if ω = ω ref = 1, then y = 0.

The following result provides an explicit expression for x k n in terms of y, and is useful to prove the results below. Lemma 6.2. For any integer n ≥ 1, define (0) n = C n and, for any p ∈ {1, . . . , n}, define (p) n = C n y p n p .

Then, for any n ≥ 1 and any k ∈ {1, . . . , n}, we have

x k n = k p=0 (p) n x 0 . (6.5) Furthermore, we have x n n = F n x 0 = x(n∆t).
Proof. Using (6.2) and (6.3), the tuple z n+1 := (x 0 n+1 , x 1 n+1 , . . . , x k n+1 ) T ∈ R k+1 , which contains the (n + 1) th time-iteration for the first k parareal iterations, evolves according to the matrix evolution

z n+1 = Lz n , with L =      C 0 . . . 0 F -C C . . . . . . 0 F -C C     
, where L ∈ R (k+1)×(k+1) . Using this relation and denoting z 0 = (x 0 , . . . , x 0 ) T the initial condition, the iterative formula for z n is

z n = L n z 0 .
Since L is a lower triangular matrix, we find

L n =        C n 0 0 . . . 0 (1) n C n 0 . . . 0 (2) n (1) n C n . . . . . . . . . . . . (k) n . . . (2) n (1) n C n       
, where, for any p ∈ {1, . . . , k},

(p) n = i1,...,ip=0 i1+...+ip≤n-p C i1 (F -C)C i2 (F -C) . . . C ip (F -C) C n-p- p α=1 iα .
Using that C and F -C are scalar numbers (a fact that would not hold for a true Langevin equation) and therefore commute, we get

(p) n = i1,...,ip=0 i1+...+ip≤n-p C n-p (F -C) p = C n-p (F -C) p n! p!(n -p)! = C n F C -1 p p .
Using the definition (6.4) of y, we obtain the expression of (p) n claimed in the lemma statement and thus (6.5).

In the particular case when k = n, we recover the fine propagator, since

x 0 n p=0 (p) n = x 0 C n n p=0 F C -1 p 1 n-p n p = C n F C -1 + 1 n x 0 = F n x 0 .
This concludes the proof of Lemma 6.2.

Using Lemma 6.2, the relative error R n,k can be written as

R n,k = 1 F n x 0 n p=0 (p) n x 0 - k p=0 (p) n x 0 = 1 F n n p=k+1 C n y p n p = n! (1 + y) n n p=k+1 y p Γ(p + 1)Γ(n -p + 1) , (6.6)
where Γ is the gamma function satisfying Γ(n + 1) = n! for any integer n. Since we are interested in the analysis of R n,k as a function of k for fixed n, the term inside the summation in the right-hand side of (6.6) is the key term that dictates the behaviour of R n,k . The following lemma (the proof of which is postponed until Appendix B) provides a complete characterisation of this term.

Lemma 6.3. Consider some integer n > 1 and some real number y ∈ (-1, 0) ∪ (0, ∞). For p ∈ {1, . . . , n}, define

η(p) = |y| p Γ(p + 1)Γ(n -p + 1)
. (6.7)

We introduce the digamma function defined by ψ = Γ Γ . Then y satisfies one (and only one) of the following three cases (hereafter, log denotes the natural logarithm function):

1. If y satisfies log |y| < ψ(2) -ψ(n), (6.8)
then η(p) is a strictly decreasing sequence in p ∈ {1, . . . , n}.

2. If y satisfies log |y| > ψ(n + 1) -ψ(1), (6.9) then η(p) is a strictly increasing sequence in p ∈ {1, . . . , n}. As pointed out in the proof below, the function ψ is increasing on (0, ∞). Assuming that n > 2, the condition (6.8) (resp. (6.9)) implies that y ∈ (-1, 0) ∪ (0, 1) (resp. y > 1). Remark 6.4 below discusses in more details the two conditions (6.8) and (6.9), and shows that the values of y such that either (6.8) or (6.9) is satisfied become somehow "less and less likely" when n increases. In the case (6.10) (which is hence the most "frequent" case when n is large), we see that (6.12) means that p essentially scales (when n 1) as

If

p ≈ n|y| 1 + |y| . (6.13)
Remark 6.4. We discuss here the conditions (6.8) and (6.9), in terms of y (and hence of the parameter ω of the coarse model) and n. In view of (B.2) below, we note that ψ(n) is equivalent to log n when n → ∞.

More precisely, there exists some c > 0 such that ∀n ∈ N , -c + log n ≤ ψ(n) ≤ c + log n.

Condition (6.8) hence yields not only that |y| < 1, but also that |y| ≤ c/n for some universal c > 0. The parameter y is thus very close to 0, and hence the coarse and the fine models are close one to each other. Using y = exp (ω -1)∆t -1, it follows that (6.8) is equivalent to

ω ∈ 1 + 1 ∆t log 1 -exp ψ(2) -ψ(n) , 1 ∪ 1, 1 + 1 ∆t log 1 + exp ψ(2) -ψ(n) .
For n = 1000 and ∆t = 0.05, which are typical values used in this article, this rewrites (after rounding) as ω ∈ (0.97, 1) ∪ (1, 1.03), which thus corresponds to ω being very close to the reference value ω ref = 1. If n increases, this domain shrinks even further.

Similarly, Condition (6.9) not only yields that y > 1, but also that y ≥ c n for some universal c > 0. The parameter y is thus very large, and the coarse and the fine models are very far apart each other. Using again y = exp (ω -1)∆t -1, it follows that (6.9) for y ∈ (1, ∞) holds if and only if

ω > 1 + 1 ∆t log 1 + exp ψ(n + 1) -ψ(1) .
As before, for n = 1000 and ∆t = 0.05, this condition rewrites as ω > 150.72, namely ω taking very large values compared to the reference value ω ref = 1. If n increases, this threshold value increases further.

For the sake of completeness, and to provide some rough orders of magnitude, we note that the value y = 1 corresponds to F = 2C, and thus ω = 1 + log 2 ∆t . For ∆t = 0.05, this corresponds to ω = 14.86.

Using Lemma 6.3, the following result, which is the main result of this section, describes the behaviour of the relative error R n,k as a function of k for fixed n. Its proof is postponed until Appendix C. Proposition 6.5. Fix n ∈ N with n > 1. If y ∈ (0, ∞), then the sequence {0, . . . , n -1} k → R n,k is strictly decreasing and we have the bound

0 < R n,k ≤ 1 - 1 (1 + y) n .
(6.14) Furthermore, we have three (and only three) distinct possibilities:

1. If y satisfies (6.8), then k → R n,k is convex for k ∈ {1, . . . , n -1}.

2. If y satisfies (6.9), then k → R n,k is concave for k ∈ {1, . . . , n -1}.

3. If there exists p ∈ [1, n] such that (6.10) is satisfied, then k → R n,k is concave for k ∈ {1, . . . , p } and convex for k ∈ { p , . . . , n -1}, where we recall that p satisfies the bounds (6.12).

If y ∈ (-1, 0), then the series R n,k = n! (1 + y) n n p=k+1 (-1) p η(p) is an alternating series and we have two (and only two) distinct possibilities:

1. If y satisfies (6.8), then R n,2 ≤ 0 and R n,1+2 ≥ 0 for any n and . The sequences → |R n,2 | and → |R n,1+2 | are both decreasing and we have the bound n! q!(n -q)! |y| q (1 + y) n .

|R n,k | ≤ max(|R n,0 |, |R n,1 |) = max 1 (1 + y) n -1, 1 - 1 (1 + y) n - ny (1 + y) n (6.
(6.17) Some estimate of the right-hand side of (6.17), in terms of y and n, is given in (6.18)-(6.19) below.

When y > 0, the estimate (6.14) shows that, when n is fixed (i.e. the time-horizon is fixed) and y gets closer to 0 (i.e. the coarse and the fine models get closer one to each other), then the bound on the relative error becomes smaller (and converges to 0 when y → 0). Conversely, for a given y > 0, the estimate (6.14) shows that, when n increases, the relative error remains bounded but may be not small. This is confirmed by numerical experiments, as illustrated below (see Figures 12a and12b and the corresponding discussion).

When y < 0 and for a fixed value of n, Condition (6.8) is satisfied for sufficiently small |y|. The estimate (6.15) again shows that the bound on the relative error converges to 0 when y → 0.

When we consider a fixed value of y < 0 and when n increases, then Condition (6.8) cannot be satisfied for sufficiently large n. Condition (6.10) is thus satisfied. The error (in absolute value) is maximal for k equal to p or p ± 1, as shown by (6.16). This is confirmed by numerical experiments (see Figure 12c below). The estimate (6.17) provides a lower bound on this error. Recasting (6.13) in the form p ≈ α n for α = |y|/(1 + |y|) and using Stirling formula, we see that the right-hand side of (6.17) is of the order of 1 2 min

q∈{ p , p } n! q!(n -q)! |y| q (1 + y) n ≈ 1 2 √ 2πn (n/e) n √ 2πp (p/e) p 2π(n -p) ((n -p)/e) n-p |y| p (1 + y) n ≈ 1 2 n n 2πα(1 -α) n (α n) α n ((1 -α) n) (1-α) n |y| α n (1 + y) n = 1 2 1 2πα(1 -α) n α α n (1 -α) (1-α) n |y| α n (1 + y) n = C α √ n (D α,y ) n , (6.18)
for some constant C α (only depending on α, which is itself a function of |y|) and with

D α,y = 1 α α (1 -α) (1-α) |y| α 1 + y = 1 α α (1 -α) (1-α) α α (1 -α) α 1 1 -α/(1 -α) = 1 1 -2α , (6.19)
where we have used that |y| = α/(1 -α) to obtain the second equality. Recalling that y ∈ (-1, 0), we have α ∈ (0, 1/2), and thus D α,y > 1. We thus infer from (6.18) that the maximum value of the error goes to ∞ when n → ∞, exponentially fast with respect to n. The equation (6.19) shows that the rate is increasing when |y| ∈ (0, 1) increases. All these facts are confirmed by numerical experiments (see Figure 12d below).

Figure 12 illustrates the findings of Proposition 6.5. Except stated otherwise, we work with the parameters n = 1000 and ∆t = 0.05. On the top left (Figure 12a), we plot k → R n,k as a function of k for various ω. All the chosen values of ω are larger than 1, hence y > 0 in all these cases. The value ω = 1.02 falls in the case when (6.8) holds (recall the numerical values from Remark 6.4), and hence R n,k is convex with respect to k (and in practice very close to 0 for most values of k). The value ω = 160 falls in the case when (6.9) holds (recall again Remark 6.4), and hence k → R n,k is concave (and in practice very close to 1 for most values of k). The remaining choices of ω correspond to intermediate values of y, falling in the case when (6.10) holds, and where R n,k is concave up to some p and then convex. We have shown with (6.13) that p essentially scales as n|y|/(1 + |y|). For instance, when ω = 14.86, we have y = 1 and we thus expect p to be of the order of n/2. This is indeed what we observe, with an inflection point of the function k → R n,k close to k = 500. When ω becomes larger (resp. smaller) than 14.86, y increases (resp. decreases) and thus n|y|/(1 + |y|) increases (resp. decreases) as well. We thus expect the inflection point to move to larger (resp. smaller) values of k, which is indeed what we observe in Figure 12a.

On the top right (Figure 12b), we again plot k → R n,k as a function of k, for the same ω as in Figure 12a, but for a value of n = 2000 twice larger. Comparing Figures 12a and12b, we check that the inflection point lies at values twice larger when n is twice as large. We have also performed computations with n = 500 and n = 1500 (results not shown), which are also in perfect agreement with (6.13).

On the bottom left (Figure 12c), we plot k → R n,k for ω = 0.1 and n = 1000. We are in the case when y ∈ (-1, 0) and Condition (6.10) holds for some p. We compute that y ≈ -0.044 and thus, according to (6.12), we have 41 ≤ p ≤ 43. We have stated in Proposition 6.5 that the maximum oscillation of R n,k occurs for k ≈ p, and that |R n,k | is maximal also for k ≈ p (see (6.16)), which is indeed what we observe in Figure 12c.

On the bottom right (Figure 12d), we consider several values of ω all smaller than 1 (hence corresponding to y ∈ (-1, 0)) and we consider values of n sufficiently large such that Condition (6.10) holds for all the values of ω that we consider (this occurs for any n ≥ 700 for ω = 0.9, and for any n ≥ 75 for ω = 0.1). We then plot R n, p as a function of n, where p is approximated by (6.13). The figure shows that the maximum error dramatically increases as a function of n: as expected from (6.18), it grows exponentially with respect to n, with a rate which, as expected from (6.19), decreases when ω gets closer to ω ref = 1.

Overall, the best possible case is when the coarse and fine models and the time window are such that Condition (6.8) holds. If y > 0, this corresponds to a convex positive decreasing function R n,k (and thus an error which converges very fast to 0). If y < 0, the function R n,k has decreasing oscillations and we essentially have that |R n,k | is decreasing. In both cases, the error decreases when the parareal iterations proceed. As pointed out in Remark 6.4, Condition (6.8) only holds for y close to 0, i.e. ω close to ω ref = 1 (and thus coarse and fine models close to each other), or for small enough values of n. In the intermediate case when (6.10) holds, and if y > 0, there is an inflection point in the behaviour of R n,k at k ≈ p. The error R n,k is always bounded, but it may take many parareal iterations to converge to zero (this is also true in the case (6.9)). If y < 0, then neither R n,k nor |R n,k | is decreasing, and the error can reach extremely large values.

A Explicit formulae for LJ7-2d

The fine potential V f for LJ7-2d is the sum of Lennard-Jones contributions:

V f (q) = 1 2 i,j∈{1,...,7} i =j φ f (|q i -q j |), φ f (r) = r -12 -2 r -6 ,
We successively study these cases. Since f is increasing, the first case corresponds to log |y| < f (1) = ψ(2) -ψ(n), and we recover (6.8). In that case, we have η (p) < 0 for any p ∈ [1, n], and thus η is strictly decreasing on [1, n].

Since f is increasing, the second case corresponds to log |y| > f (n) = ψ(n + 1) -ψ(1), and we recover (6.9). In that case, we have η (p) > 0 for any p ∈ [1, n], and thus η is strictly increasing on [1, n].

Again using the fact that f is increasing, the third case corresponds to the fact that there exists some p ∈ [1, n] such that log |y| = f (p) = ψ(p + 1) -ψ(n -p + 1) and we recover (6.10). Using (B.1), we next obtain (6.11).

We now establish explicit bounds on p. Using the fact that

ψ(p + 1) = -γ EM + p i=1 1 i for any positive integer p, (B.2)
where γ EM is the Euler-Mascheroni constant, we have

f ( p ) =                    - n-p i= p +1 1 i if 2 p < n 0 if 2 p = n p i=n-p +1 1 i if 2 p > n and likewise f ( p ) =                    - n-p i= p +1 1 i if 2 p < n 0 if 2 p = n p i=n-p +1 1 i if 2 p > n . (B.3)
We now write some easy and useful estimates. For any integer k such that k + 1 < n -k + 1, we have the following upper and lower bounds:

log n -k + 1 k + 1 = n-k+1 k+1 dx x = n-k i=k+1 i+1 i dx x ≤ n-k i=k+1 1 i , n-k i=k+1 1 i = n-k-1 j=k 1 j + 1 ≤ n-k-1 j=k j+1 j dx x = n-k k dx x = log n -k k . (B.4)
Proceeding similarly, we find, for any integer and hence (6.12).

k such that k + 1 > n -k + 1, that log k + 1 n -k + 1 = k+1 n-k+1 dx x = k i=n-k+1 i+1 i dx x ≤ k i=n-k+1 1 i , k i=n-k+1 1 i = k-1 j=n-k 1 j + 1 ≤ k-1 j=n-k j+1 j dx x = k n-k dx x = log k n -k . (B.
We next turn to the case (II), when n < 2 p . We then have log |y| ≥ f ( p ) > 0, thus |y| > 1 and hence y > 1 (since y ∈ (-1, 0) ∪ (0, ∞)). Inserting (B.8) and (B.9) into (B.10) and proceeding as above, we obtain the bounds

(n -1)|y| -1 1 + |y| ≤ p ≤ (n + 1)|y| -1 1 + |y| , n|y| 1 + |y| ≤ p ≤ (n + 2)|y| 1 + |y| ,
and therefore (n -1)|y| -1 1 + |y| ≤ p ≤ (n + 2)|y| 1 + |y| , which again implies (6.12).

We eventually turn to the case (III), when 2 p ≤ n ≤ 2 p . If p is an integer, then 2p = n. In addition, in view of (B.3), we have f (p) = 0, and thus (B.10) implies that |y| = 1, and hence y = 1 (since y ∈ (-1, 0) ∪ (0, ∞)). The bound (6.12) again holds.

If p is not an integer, then either (i) n = 2 p or (ii) n = 1 + 2 p or (iii) n = 2 + 2 p = 2 p . In the case (i) (i.e. when n = 2 p ), we have, in view of (B.3), that f ( p ) = 0 and f ( p ) = 1/ p + 1/( p -1). . This shows that (6.12) again holds.

In a similar fashion, in the case (iii) (i.e. when n = 2 p ), we have, in view of (B.3), that f ( p ) = 0 and f ( p ) = -1/( p + 2) -1/( p + 1). In view of (B.10), this implies that 0 ≥ log |y| ≥ -

1 p + 2 - 1 p + 1 ≥ -log p + 2 p = -log n + 2 n -2
and hence 1 ≥ |y| ≥ n -2 n + 2 .

The left and right-hand sides of (6.12) then satisfy . This shows that (6.12) again holds.

We eventually consider the case (ii), when n = 1 + 2 p = 2 p -1. In view of (B.3), we have that f ( p ) = -1/( p + 1) and f ( p ) = 1/ p . In view of (B.10), this implies that

-log n + 1 n -1 = -log p + 1 p ≤ - 1 p + 1 ≤ log |y| ≤ 1 p ≤ log p p -1 = log n + 1 n -1 and hence n -1 n + 1 ≤ |y| ≤ n + 1 n -1 .
and each term η(j -1) -η(j) (for odd values of j) is positive. If n is even, then R n,1 = n! (1 + y) n η(2) -η(3) + . . . + η(n -2) -η(n -1) + η(n) and again, each term η(j -1) -η(j) (for odd values of j) is positive, while η(n) is also positive. We thus have that R n,1 > 0 for any n ≥ 2. With similar arguments, we can show that R n,k ≥ 0 for any odd k, and that the sequence → R n,1+2 (namely, the relative error after an odd number of parareal iterations) is decreasing.

Since Case B2. We eventually consider the case when y ∈ (-1, 0) and (6.10) is satisfied. Using Lemma 6.3, we know that s → η(s) is increasing (resp. decreasing) for s ≤ p (resp. s ≥ p). In view of (C.1), we get that For any q ∈ {1, . . . , n}, we have max k∈{0,...,n} |R n,k | ≥ |R n,q | ≥ |R n,q -R n,q-1 | -|R n,q-1 | ≥ |R n,q -R n,q-1 | -max |R n,q -R n,q-1 |. Since this holds for any q ∈ {1, . . . , n}, we deduce, using our previous result on the largest oscillation, that max k∈{0,...,n}

|R n,k -R n,k-1 | = n! (1 + y) n η(k).
|R n,k | ≥ 1 2 min q∈{ p , p }
|R n,q -R n,q-1 | = 1 2 min q∈{ p , p } n! q!(n -q)! |y| q (1 + y) n , which is (6.17). We are now left with showing (6.16).

Case B2.1. In order to study |R n,k |, we first assume that n and p are even. We claim that the sequence → R n,2 is decreasing when 0 ≤ ≤ p /2 and increasing when p /2 ≤ ≤ n/2. Indeed, using (C.1), we have R n,2 -R n,2+2 = n! (1 + y) n η(2 + 2 ) -η(1 + 2 ) , which is positive when 0 ≤ ≤ p /2 -1 and negative when ≥ p /2 (since then 1 + 2 ≥ 1 + p > p). In addition, for the final value = n/2, we of course have R n,2 = 0.

Likewise, we claim that the sequence → R n,1+2 is increasing when 0 ≤ ≤ p /2 -1 and decreasing when p /2 ≤ ≤ n/2. Indeed, using again (C.1), we have R n,1+2 -R n,3+2 = n! (1 + y) n η(2 + 2 ) -η(3 + 2 ) , which is negative when 0 ≤ ≤ p /2-2 and positive when ≥ p /2. We also note that we cannot compare the values of R n,1+2 when = p /2 -1 and when = p /2, since their difference is R n, p -1 -R n, p +1 = n! (1 + y) n η( p ) -η(1 + p ) , the sign of which is unknown. In addition, for the final value = n/2-1, we have R n,1+2 = n! (-1) n η(n)

(1 + y) n > 0 since we have assumed n to be even.

We now identify max k∈{0,...,n}

|R n,k |, supposing first that R n,0 ≤ 0. Since → R n,2 is first decreasing and next increasing and since it vanishes for = n/2, we have R n,2 ≤ 0 for all 0 ≤ ≤ n/2, and therefore max ∈{0,...,n/2}

|R n,2 | = R n, p . If R n,1 ≥ 0, then using that → R n,1+2 is first increasing and next decreasing and that it is positive for = n/2 -1, we have R n,1+2 ≥ 0 for all 0 ≤ ≤ n/2 -1. This implies that max ∈{0,...,n/2-1}

|R n,1+2 | = max R n, p -1 , R n, p +1 , and we thus obtain (6.16). If R n,1 < 0, then using again that → R n,1+2 is first increasing and next decreasing and that it is positive for = n/2 -1, we have that max ∈{0,...,n/2-1}

|R n,1+2 | = max R n,1 , R n, p -1 , R n, p +1 . We also observe that R n,0 -R n,1 = -n! (1 + y) n η(1) < 0, and thus 0 > R n,1 > R n,0 . This implies that R n,1 < R n,0 ≤ max ∈{0,...,n/2}

|R n,2 |, and we thus again obtain (6.16).

We next identify max |R n,1+2 |. We thus deduce (6.16).

Case B2.2. We now revisit the above arguments studying |R n,k |, under the assumption that n is even and that p is odd. We then have that the sequence → R n,2 is decreasing when 0 ≤ ≤ ( p -1)/2 and increasing when ( p + 1)/2 ≤ ≤ n/2. We also note that we cannot compare the values of R n,2 when = ( p -1)/2 and when = ( p + 1)/2, and that, for the final value = n/2, we have R n,2 = 0. Similarly, the sequence → R n,1+2 is increasing when 0 ≤ ≤ ( p -1)/2 and decreasing when ( p -1)/2 ≤ ≤ n/2. For the final value = n/2 -1, we again have R n,1+2 > 0.

Using the same arguments as above for identifying max k∈{0,...,n}

|R n,k | again leads to (6.16).

Case B2.3. We eventually consider the case when n is odd:

• if p is even, then the sequence → R n,2 is decreasing when 0 ≤ ≤ p /2, increasing when p /2 ≤ ≤ (n -1)/2, and is negative for the final value = (n -1)/2. The sequence → R n,1+2 is increasing when 0 ≤ ≤ p /2 -1, decreasing when p /2 ≤ ≤ (n -1)/2, and vanishes for the final value = (n -1)/2.

Figure 1 :

 1 Figure 1: Local minima of the potential energy for the Lennard-Jones cluster with seven atoms in 2D. Also indicated are the associated potential energy levels.

Figure 2 :

 2 Figure 2: Plot of E f (k, N, γ, β) as a function of k for various γ. Left: Har-1d with ∆t = 0.05, ω = 0.1, β = 3 and N = 200. Right: LJ7-2d with ∆t = 0.01, β = 10 and N = 100.

Figure 3 :

 3 Figure 3: Plot of E f (k, N, γ, β) as a function of k for various γ. Left: Har-1d with ∆t = 0.05, ω = 0.1, β = 3 and N = 3000. Right: LJ7-2d with ∆t = 0.01, β = 10 and N = 500. The gap in the right figure for γ = 30 corresponds to undefined (NaN) values for E f .

Figure 4 :

 4 Figure 4: Plot of fine and parareal trajectories as a function of the time index for Har-1d with ∆t = 0.05, N = 3000, β = γ = 3 and ω = 0.1. Top left: plot of the fine trajectory. Top right: plot of (q 40n , p 40 n ) in log-scale for k = 40 (the gaps in the graph correspond to negative values). Bottom left: plot of (q 80 n , p 80 n ) in log-scale for k = 80 (the gaps in the graph correspond to negative values). Bottom right: plot of (q 150 n , p 150 n ) for k = 150. While the trajectory is stable at larger k, it blows up for intermediate k.

Figure 5 :

 5 Figure 5: Plot of E f (k, 3000, 3, 3) as a function of k for Har-1d with ω = 0.1, ∆t = 0.05, N = 3000 and β = γ = 3 with different propagators. The classical parareal curve corresponds to the potentials defined by (2.4). The trust-region curve corresponds to the potentials defined by (3.5).

Figure 6 :

 6 Figure6: Left and Center: Plot of E(k, N, γ, β) and E f (k, N, γ, β) for Har-1d for the first 50 (resp. 500) parareal iterations (the left plot is a zoom of the center plot for the first parareal iterations; it clearly shows the maximum value attained by the error E, a value which is to be compared with the threshold δ expl introduced in Section 4 below). Right: Plot of E(k, N, γ, β) and E f (k, N, γ, β) for LJ7-2d for the first 350 iterations. Note that an increase in E(k, N, γ, β) beyond a threshold indicates explosion of parareal trajectory.

Figure 7 :

 7 Figure 7: Plot of gain(N, γ, β) as a function of N for various γ. Left: Har-1d with ∆t = 0.05 and β = 3. Right: LJ7-2d with ∆t = 0.01 and β = 10.

Figure 8 :

 8 Figure 8: Plot of gain(5000, 3, 3) for Har-1d (with β = γ = 3, ∆t = 0.05 and N = 5000, hence T = N ∆t = 250) as a function of ω. Left: ω < 1. Right: ω > 1. Note the different scale on the vertical axis of the two curves.

Algorithm 1 :}:= x cur n ; 7 Compute9if δ > δ expl then 19 N 21 N

 171921 Adaptive parareal algorithmNumerical parameters: N , δ conv , δ expl 1 Set N init := 0, N final := N and δ := (δ conv + δ expl )/2; 2 Compute {x cur n } Ninit≤n≤N final : x cur 0 := x 0 and sequential update x cur n+1 := C ∆t (x cur n ); 3 Set Cost := 0;4 while N init < N do 5 while δ ∈ [δ conv , δ expl ] do 6 Define {x prev n Ninit≤n≤N final as x prev n {J n } Ninit≤n≤N final (in parallel): J n := F ∆t (x prev n ) -C ∆t (x prev n for n ← N init to N final -1 by 1 do 10 x cur n+1 := C ∆t (x cur n ) + J n ; 11 Compute the relative error δ = E(N init , n + 1; x prev , x cur ); 12 if δ > δ expl then 13 N expl := n + 1; 14 break; // exit the for loop if condition satisfied; we also exit the while loop since δ is too large final := N expl -1; 20 else init := N final ; 22 Compute {x cur n } Ninit≤n≤N : x cur Ninit := x cur N final and sequential update x cur n+1 := C ∆t (x cur n ); 23 N final := N ; 24 end 25 Set δ := (δ conv + δ expl )/2; 26 end Output of the algorithm: {x cur n } 0≤n≤N and Cost

Figure 9 :

 9 Figure9: We use the Adaptive parareal algorithm for the Har-1d case with ω = 0.1 and β = γ = 3. For two different values of the time horizon (left: T = 250; right: T = 1500), and different choices of δ expl , we plot the list of the sizes of the time-slabs found by the algorithm. For instance, when N = 5000 and δ expl = 1, the total time interval [0, T ] is eventually divided into 31 slabs, each of which of size of the order of 160 time-steps. The 12-th slab is the longest and its length is of the order of 350 time-steps. The 31st slab, which is the last one needed to reach the final horizon, is the shortest with a length of the order of 50 time-steps. When N = 5000 and δ expl = 5, the total time interval [0, T ] is eventually divided into 4 slabs, of length close to 1300, 1500, 1200 and 1000 time-steps.

NNFigure 10 :Figure 11 :

 1011 Figure 10: Plot of gain(N, γ, β) as a function of N for various γ for the classical and the adaptive parareal algorithm. Top: Har-1d with ∆t = 0.05, β = 3 and γ = 3 (left) or γ = 5 (right). Bottom: LJ7-2d with ∆t = 0.01, β = 10 and γ = 16 (left) or γ = 30 (right).

Figure 12 :

 12 Figure12: We work with ∆t = 0.05. Top left: Plot of k → R n,k as a function of k for various ω, for n = 1000 (all these ω are larger than 1, hence y > 0 in all cases; the three possible cases (6.8), (6.9) and (6.10) are represented). Top right: Plot of k → R n,k as a function of k for various ω, for n = 2000. Bottom left: Plot of k → R n,k for ω = 0.1 and n = 1000 (y < 0 and (6.10) holds). Bottom right: Plot of |R n, p | as a function of n for various ω (in these cases, y < 0).

  → |R n,2 | and → |R n,1+2 | are both decreasing, we have max k∈{0,...,n-1}|R n,k | = max(|R n,0 |, |R n,1 |), y) n (1 + y) n -1 = 1 (1 + y) n -1 and 0 < R n,1 = 1 (1 + y) n n p=2 n! y p p!(n -p)! = 1 (1 + y) n (1 + y) n -1 -ny = 1 -1 (1 + y) n -ny (1 + y) n ,thus the bound (6.15).

  The largest oscillation for R n,k is thus given by max k∈{1,...,n}|R n,k -R n,k-1 | = n! (1 + y) n max k∈{1,...,n} η(k) = n! (1 + y) n max η( p ), η( p ) = max R n, p -R n, p -1 , R n, p -R n, p -1 .

  k∈{0,...,n} |R n,k |, and therefore max k∈{0,...,n} |R n,k | ≥ 1 2

  k∈{0,...,n} |R n,k | under the assumption that R n,0 > 0. Using the same arguments as above, we have max ∈{0,...,n/2} |R n,2 | = max R n,0 , R n, p . Since R n,1 > R n,0 , we have R n,1 > 0, and thus max ∈{0,...,n/2-1} |R n,1+2 | = max R n, p -1 , R n, p +1 . In addition, we have 0 < R n,0 < R n,1 ≤ max ∈{0,...,n/2-1}

  there exists a real number p ∈ [1, n] such that

					log |y| = ψ(p + 1) -ψ(n -p + 1),				(6.10)
	then
				η(p)	is strictly increasing for p ∈ {1, . . . , p }, is strictly for p ∈ { p , . . . , n},				(6.11)
	where • (resp. • ) is the floor (resp. ceiling) function. In this case, p satisfies the following bound:
	min	n|y| -2 1 + |y|	,	(n -1)|y| -1 1 + |y|	≤ p ≤ max	(n + 1)|y| + 1 1 + |y|	,	(n + 2)|y| 1 + |y|	.	(6.12)

  15)for any k ∈ {0, . . . , n}. Furthermore, the oscillations of R n,k , namely the sequence k → |R n,k+1 -R n,k |, are decreasing when k increases.2. If there exists p ∈ [1, n] such that (6.10) is satisfied, then the largest oscillation of R n,k when k varies

	(namely max k∈{1,...,n}	|R n,k -R n,k-1 |) is either R n, p -R n, p -1 or R n, p -R n, p -1 , where we recall
	that p satisfies (6.12). In addition, we have the bound
	∀k ∈ {0, . . . , n},	|R n,k | ≤ max R n, p -1 , R n, p , R n, p +1 ,	(6.16)
	and the lower bound	max k∈{0,...,n}	|R n,k | ≥	1 2	min q∈{ p , p }

  Since f is increasing, it follows from (6.10), that we write in the form log |y| = f (p), that

	and likewise, if 2 p > n, then f ( p ) > 0 and it satisfies
			log	p + 1 n -p + 1	≤ f ( p ) ≤ log	p n -p	.	(B.9)
						f ( p ) ≤ log |y| ≤ f ( p ).	(B.10)
	We now distinguish three (mutually exclusive) cases:
	(I) either 2 p ≤ 2 p < n,
	(II) or n < 2 p ≤ 2 p ,
	(III) or 2 p ≤ n ≤ 2 p ,
	that we successively study.
	In the case (I), we have 2 p < n, hence log |y| ≤ f ( p ) < 0 and therefore |y| < 1. Inserting (B.6)
	and (B.7) into (B.10), we find
			log	p n -p	≤ log |y| ≤ log	p + 1 n -p + 1	,
	and thus				p ≤	n|y| 1 + |y|	,	(n + 1)|y| -1 1 + |y|	≤ p .
	Using p ≥ p -1, we infer that
	n|y| -2 1 + |y|	≤ p ≤	n|y| 1 + |y|	,		(n + 1)|y| -1 1 + |y|	≤ p ≤	(n + 1)|y| + 1 1 + |y|	,
	and hence, since p and p are integer numbers, we get
	n|y| -2 1 + |y|	≤ p ≤	n|y| 1 + |y|		,		(n + 1)|y| -1 1 + |y|	≤ p ≤	(n + 1)|y| + 1 1 + |y|	.
	Using that p ≤ p ≤ p , we eventually get the bound
					n|y| -2 1 + |y|		≤ p ≤	(n + 1)|y| + 1 1 + |y|	,	5)
	Using (B.3) and (B.4), we see that, if 2 p < n, then f ( p ) < 0 and it satisfies
			log	p n -p	≤ f ( p ) ≤ log	p + 1 n -p + 1	.	(B.6)
	Similarly, if 2 p < n, then f ( p ) < 0 and it satisfies
			log	p n -p	≤ f ( p ) ≤ log	p + 1 n -p + 1	.	(B.7)
	Alternatively, if 2 p > n, we infer from (B.3) and (B.5) that f ( p ) > 0 and that it satisfies
			log	p + 1 n -p + 1	≤ f ( p ) ≤ log	p n -p	,	(B.8)

  In view of (B.10), this implies that

			0 ≤ log |y| ≤	1 p	+	1 p -1	≤ log	p p -2	= log	n + 2 n -2
	and hence						1 ≤ |y| ≤	n + 2 n -2	.
	The left and right-hand sides of (6.12) then satisfy
		min	n|y| -2 1 + |y|	,	(n -1)|y| -1 1 + |y|	=	(n -1)|y| -1 1 + |y|	≤	n 2	= p ≤ p
	and	max	(n + 1)|y| + 1 1 + |y|	,	(n + 2)|y| 1 + |y|	=	(n + 2)|y| 1 + |y|	≥	n + 2 2	= p + 1 ≥ p,
	where we have used that |y| ≥ 1 (for the first equalities) and the fact that |y| →	(n -1)|y| -1 1 + |y|	=	n 1 + 1/|y|	-1
	and |y| →	(n + 2)|y| 1 + |y|	are both increasing functions (for the first inequalities)

Acknowledgements

The work of FL and US was supported in part by the ANR project CINE-PARA (ANR-15-CE23-0019-06). US also acknowledges support from the Alexander von Humboldt foundation and from DFG under Germany's Excellence Strategy-MATH+: The Berlin Mathematics Research Center (EXC-2046/1)-project ID:390685689 (subproject EF4-4). This work was also partially funded by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 810367; project EMC2) and from the European High-Performance Computing Joint Undertaking (JU) under grant agreement No 955701 (project Time-X). The JU receives support from the European Union's Horizon 2020 research and innovation programme and Belgium, France, Germany, Switzerland. The authors would like to thank Claude Le Bris, Yvon Maday, Danny Perez and Arthur Voter for stimulating interactions, and Olga Gorynina and Gabriel Stoltz for their remarks on a draft version of this article.

where | • | is the d-dimensional Euclidean norm (here, d = 2). The corresponding force on the system is F f = -∇V f , and the force on the i-th particle (here onwards we do not explicitly state that i, j ∈ {1, . . . , 7}) is

The coarse potential is a harmonic approximation of V f around the well where the system initially lies. To calculate this coarse force, we need the Hessian of the fine potential V f , which is given (for any α, β ∈ {x, y} and any i, j ∈ {1, . . . , 7}) by

where δ αβ = 1 if α = β and zero otherwise, and

Using q for the position of the well corresponding to the initial condition of the system, the (harmonic) coarse potential is given by

and the corresponding force on the j-th particle in the β-th coordinate is

B Proof of Lemma 6.3

Differentiating the function [1, n] p → η(p) and using the digamma function ψ(p) = Γ (p)/Γ(p), we find

Using that ψ > 0 (a fact which directly stems from the series expansion of the polygamma function of order one ψ 1 = ψ ), we observe that the mapping

We can thus distinguish three (mutually exclusive) cases:

• or inf

Since the functions that map |y| to the left and right-hand sides of (6.12) are both increasing, we obtain

which again implies (6.12). This concludes the proof of Lemma 6.3.

C Proof of Proposition 6.5

The proof falls in two cases.

Case A. We first consider the case when y ∈ (0, ∞). In view of (6.6), the sequence k → R n,k is a decreasing sequence and therefore, for any k ∈ {0, . . . , n -1}, we have

which yields (6.14). In view of (6.6) and (6.7), we have

η(p). The convexity properties of R n,k directly follow from the monotonicity properties of p → η(p) established in Lemma 6.3.

Case B. We next consider the case when y ∈ (-1, 0). We have pointed out below Lemma 6.3 that Condition (6.9) cannot be satisfied in that case, and we thus only have to consider (6.8) or (6.10). In view of (6.6) and (6.7), we have

which is an alternating series.

Case B1. If y ∈ (-1, 0) satisfies (6.8), we know from Lemma 6.3 that η is a decreasing function. In view of (C.1), we obtain that the oscillations of R n,k are decreasing when k increases. We claim that R n,0 < 0. Indeed, if n is even, then we write

and each term -η(j -1) + η(j) (for even values of j) is negative. If n is odd, then

and again, each term -η(j -1) + η(j) (for even values of j) is negative, while -η(n) is also negative. We hence have that R n,0 < 0 for any n ≥ 1. With similar arguments, we can show that R n,k ≤ 0 for any even k. Furthermore, the sequence → R n,2 (namely, the relative error after an even number of parareal iterations) is increasing.

We next claim that R n,1 > 0. Indeed, if n is odd, then we write R n,1 = n! (1 + y) n η(2) -η(3) + . . . + η(n -1) -η(n)