Simulating Stochastic Dynamics Using Large Time Steps

We present a novel approach to investigate the long-time stochastic dynamics of multi-dimensional classical systems, in contact with a heat-bath. When the potential energy landscape is rugged, the kinetics displays a decoupling of short and long time scales and both Molecular Dynamics (MD) or Monte Carlo (MC) simulations are generally inefficient. Using a field theoretic approach, we perform analytically the average over the short-time stochastic fluctuations. This way, we obtain an effective theory, which generates the same long-time dynamics of the original theory, but has a lower time resolution power. Such an approach is used to develop an improved version of the MC algorithm, which is particularly suitable to investigate the dynamics of rare conformational transitions. In the specific case of molecular systems at room temperature, we show that elementary integration time steps used to simulate the effective theory can be chosen a factor ~100 larger than those used in the original theory. Our results are illustrated and tested on a simple system, characterized by a rugged energy landscape.


I. INTRODUCTION
The investigation of a vast class of physical phenomena requires the understanding of the long-time dynamics of classical systems, in contact with a heat-bath. Examples include critical dynamics, molecular aggregation, protein folding, to name a few.
The most natural strategy to describe these processes is to integrate numerically the equations of motion, i.e. to perform MD simulations. Unfortunately, when the number of degrees of freedom is very large, or in the presence of large free energy barriers, MD approaches become extremely costly 1 , or even impracticable. The problem arises because the time scale associated with the system's local conformational changes can be many orders of magnitude smaller that the time scales of the dynamics one is interesting in studying. As a result, most of the computational time is invested in simulating uninteresting thermal oscillations.
This situation is exemplified in Fig.1, where we show the stochastic motion of a point particle, interacting with a 2-dimensional external potential. The solid line was obtained by means of a MD simulation and illustrates how, at short time scales, the dynamics of this system is dominated by fast modes associated to thermal diffusion. However, when the evolution of the system is described using much lower time resolution power, the effect of such short-time thermal fluctuations tends to average out and to become unimportant. This is evident from the comparison between the solid line and the dashed line, which was obtained by averaging over blocks of consecutive frames in the original MD trajectory. At long times, the dynamics of system is mostly sensitive to the structure of the external energy landscape, which was chosen to be spherically symmetric.
Clearly, an important question to ask is whether it is possible to develop theoretical/computational frameworks which yield directly the correct long-dynamics, but avoid investing computational time in simulating the short-time thermal oscillations. Significant progress in this direction has been recently made developing approaches based on Markov State Models 2,3,4 . A potential difficulty of such approaches resides in correct identification of the metastable states. In addition, for each different system, one needs to perform a large set of independent MD simulations in order to accurate calculation of the rate coefficients.
In this work, we present an alternative approach to simulate the dynamics over long times. We develop a rigorous effective theory which (i) yields by construction the correct long-dynamics and (ii) does not require to identify metastable states, nor to evaluate the transition matrix by MD. To our goal, we use a field theory approach, based on Renormalization Group (RG) ideas and on the notion of effective field theory 5 . Such a powerful tools have been already successfully applied to describe the low-energy dynamics of a vast variety of of quantum and statistical systems characterized by a separation of scales -see e.g. 6,7 -. To the best of our knowledge, this method has never been applied to develop an effective theory to efficiently simulate the long-time stochastic dynamics of a system in contact with a heat bath.
The main idea of our approach is to exploit the decoupling of time scales in the system in order to define a perturbative series, in which the expansion parameter is the ratio of short-over large-time scales. In such a perturbative framework, the average over the short-time fluctuations can be computed analytically, to any desired level of accuracy. The average over the fast thermal oscillations gives rise to new terms in the stochastic path integral, which represent corrections both to the interaction and to the diffusion coefficient. Such new terms implicitly take into account of the dynamics of the fast degrees of freedom, which have been integrated out from the system.
Once a finite number of such effective terms corresponding to a given accuracy have been calculated analytically, it is possible to simulate the dynamics of the system using much larger time steps. By construction, in the regime of decoupling of fast and slow modes, one is guaranteed that the effective long time theory generates the same probability distributions of the underlying, more fundamental stochastic theory. It is important to emphasize the fact that the present approach is not equivalent to simply including higher-order corrections in the Trotter expansion 8 . Indeed, the assumption of decoupling of time scales leads to further simplifications with respect to such an approach.
The paper is organized as follows. In section II, we review the path integral formulation of the Langevin dynamics and we outline the formal connection between stochastic dynamics and evolution of a quantum particle in imaginary time. Such a connection is used in section III to identify and isolate the dynamics of the fast degrees of freedom. In sections IV and V we present our perturbative scheme which allows to integrate out the fast modes and derive the effective interactions and diffusion coefficients. In section VI we discuss how the effective theory for the dynamics at long time scales can be simulated using the diffusion MC algorithm, which is briefly reviewed in appendix B. Section VI is devoted to simple examples, which illustrate how this method works in practice. In section VIII we discuss the applicability of the present approach to simulate the Langevin dynamics of molecular systems. Results and conclusions are summarized in section IX.

II. LANGEVIN DYNAMICS
We consider a system defined by a stochastic d-dimensional variable x obeying the Langevin Eq.n: where U (x) is a potential energy function, m is the mass, γ is the friction coefficient and ξ(t) is a δ−correlated Gaussian noise. In many molecular systems of interest, the acceleration term mẍ is damped at time scales of the order 10 −13 s, which much smaller than the time scale associated to local conformational changes. If such a term is dropped one obtains the so-called over-damped or velocity Langevin Eq.: where η(t) is a rescaled delta-correlated Gaussian noise, satisfying the fluctuation-dissipation relationship: The over-damped Langevin Eq. defines a Markovian process. The probability distribution P (x, t) generated by such a stochastic differential equation obeys the Fokker-Planck Eq.
By performing the substitution the Fokker-Planck Eq.(4) can be recast in the form of a Schrödinger Equation in imaginary time: where the effective "Quantum Hamiltonian" operator readŝ and V ef f (x) is called the effective potential and reads Hence, the problem of studying the diffusion of a classical particle can be mapped into the problem of determining the quantum-mechanical propagation in imaginary time of a virtual system, defined by the effective quantum Hamiltonian (7), interacting with the effective potential V ef f (x).
Let G(x f , t f |x i ) be the Green's function of the Fokker-Planck operator, subject to initial condition The interpretation of such a Green's function is the probability for the system to be in x at t, conditioned to start from x i at the initial time. Formally, such a conditional probability can be related to the "quantum" propagator of the effective Hamiltonian (7): Hence, it is immediate to derive a path integral representation of the Green's function G(x, t|x i ): where S ef f [x] is the effective "action", The pre-factor e −β/2(U (x)−U (xi) in Eq. (11) can be transformed away, noticing that dU (x) dτ =ẋU (x). One than obtains a path integral in which the statistical weight contains the Onsager-Machlup functional Eq. (12) provides an expression for the conditional probability in terms of the microscopic stochastic dynamics governing the system. It represents the starting point of the Dominant Reaction Pathway approach 9,10,11,12 , which deals with the problem of finding the most probable transition pathways between the given configurations x i and x f , which are visited at the initial and final time On the other hand, in this work we are interested in the corresponding initial value problem, i.e. we are want to develop an effective theory which yields directly the long-time evolution of the probability density P (x, t), solution of Eq. (4), starting from a given initial probability density P (x, t = 0) = ρ 0 (x). The probability density P (x, t), the Green's function G(x f , t|x i , t i ) and the initial distribution ρ 0 (x) are related by the Eq.
Hence, for positive time intervals, the conditional probability G(x, t|x i ) can be considered as the propagator associated to the stochastic Fokker-Planck Eq. (4).

III. SEPARATION OF FAST AND SLOW MODES
Without loss of generality, let us focus on the stochastic path integral (12), with periodic boundary conditions: We observe that the inverse temperature 1 β plays the role of , in the analogy with the quantum mechanical formalism. Hence, the loop expansion of the path integral (16) generates an expansion in powers of 1 β . Let us introduce the Fourier conjugate: where ω n are the Matsubara frequencies: In numerical simulations, the integration of the over-damped Langevin Eq. is performed by choosing a finite elementary time step ∆t. In frequency space, this implies the existence of an ultra-violet cut-off Ω, which is inversely proportional to ∆t: Such a relationship becomes a strict equality in the case of periodic boundary conditions, as in Eq. (16). In general, when the boundary conditions are not periodic, it represents just an order-of-magnitude estimate of the largest Fourier frequencies, which are associated to a given choice of the integration time step ∆t. Let us now introduce a parameter 0 < b < 1 and split the frequency interval (0, Ω) as (0, b Ω) ∪ (b Ω, Ω). Then the Fourier decomposition of a path contributing to (16)) can be split as: where x < (t) and x > (t) will be referred to as the slow -and fast-modes respectively: The main purpose of this work is to develop a perturbation series to systematically integrate out from the path integral the modes with frequencies ω n > bΩ. To this end, we begin by re-writing the "kinetic" term which appears in the effective action (13) of the path integral (12) as a sum of the kinetic energy of slow and fast modes: where S b denotes the shell of hard modes S b = (bΩ, Ω). Let us now consider the potential term and expand around the slow modes x < (t) The complete path integral (16) can be split in the following way: in this expression, the action functional S ef f is evaluated on the slow-modes only and depends on the original effective potential V ef f (which we also shall refer to as the "tree-level" effective potential). S > [x < (τ )] is a correction term action which accounts for the dynamics of the fast modes which are integrated out: where the S int is an effective interaction term. In such an Eq., the integration over the hard modes is performed in Fourier space, Eq. (26) is formally exact. In the next section, we evaluate the effective action S > [x < (τ )] perturbatively. The effective interaction which includes the correction coming from S > [x < ] will be referred to as the renormalized effective interaction.

IV. RENORMALIZED EFFECTIVE INTERACTION
In the previous section, we have seen that the the integration over the fast modes generates an additional term in the effective action for the slow modes: where In this section we formally perform such an integration. We begin by re-writing e −β S>[x<(τ )] as where the notation · 0 denotes the expectation value evaluated in the free theory To evaluate the matrix element e −βSint , we represent the e −βSint "operator" by its power series: Next, we expand the interaction action S int [x > + x < ] around the slow modes 1 : Notice that each term in the perturbative expansion (34) generates a new vertex, with an increasing power of the x > (τ ) field. The couplings to the fast modes depend implicitly on the time τ , through the slow modes x < (τ ). By Wick theorem, each term in the series (33) can be related to a Feynman graph with vertexes given by (35) and propagators given by -see appendix A -: The expansion (33) can be re-organized as the exponent of the sum performed over only connected diagrams: Hence, the path integral (26) for the slow modes can be given the following exact diagrammatic representation Below we give a classification of all the connected diagrams that may give a contribution to the expansion above.
Firstly note that all diagrams that involve an odd numbers of fast field vanish thanks to the Wick theorem. We are thus left with the following sets of, a priori nonvanishing, diagrams: • 1PR (one-particle-reducible) diagrams, namely diagrams that can be topologically separated into two distinct subdiagrams by cutting one internal fast mode line (propagator): they have the topology of a dumbbell. The simplest examples of dumbbell diagrams are depicted in the upper part of Fig. 2.
The main assumption of this work is the existence of a gap between slow modes and fast modes. Under such assumption all the 1PR diagrams give vanishing contributions. From the physical point of view, this can be understood as a consequence of energy conservation: in order for the total energy flowing through a vertex with a single hard mode to be conserved, at least one of the external modes has to be hard. On the other hand, our working assumption implies that all the modes in the external legs of diagrams are soft. This result can be rigorously proven for all 1PR. As an example, we explicitly compute the upper left diagram of Fig. 2. We have We note that the effective potentials depend smoothly on time, through the periodic functions This allows to perform the time integrals, which simply yield t 2 δ ωn+νn,0 δ ωn−νm,0 . Due to such delta-functions, only hard ν-modes survives, which are projected in a term These modes thus yield negligible contributions under the physical assumption of large separation of frequency scales. On the other hand, if one does not assume a separation of time scale, this diagram gives finite contribution and has to be accounted for. Note that this term has the same structure as the first correction which appears when one performs higher-order Trotter expansion 8 .
It is not difficult to check that such result holds for all 1PR diagrams, so that we can reduce our effective action to the sum of 1PI (one-particle-irreducible) diagrams, i.e. diagrams that cannot be disconnected by cutting a single internal line. They can be classified in two main groups: • 1PI "daisy" diagrams, namely diagrams with a single vertex. Such diagrams only involve equal-time hard propagators and only give rise to contributions to the renormalized effective action which are local in time: they have the topology of a daisy, hence the name. Examples of daisy diagrams are depicted in the middle part of Fig. 2. It is not difficult to compute a generic daisy diagram with K petals (propagators). It is due to the vertex with 2K hard fields and reads where ∆ ≡ δ ij ∂ i ∂ j is the Laplacian operator, and the numerical factor in front is a combinatorial factor. The sum, i.e. the equal-time propagator, can be easily performed by taking the continuum limit ω → t 2π dω that simply yields 2 so that we finally obtain where we have reinstated the diffusion coefficient D = 1/βγ. Hence, one can even formally resum all the daisy diagrams into the compact expression • 1PI non-daisy diagrams: all other non-local diagrams. The simplest examples of such diagrams are depicted in the lower part of Fig. 2. These diagrams generate contributions to the renormalized effective action that are non-local in time and give rise to infinite series of local diagrams. For example, the evaluation of the lower left diagram of Fig. 2 yields a contribution of the form: where the 2 in front is a combinatorial factor. After Fourier transforming the potentials (see discussion below eq. (39)), the integrals over times yield t 2 δ ωm+ωn+νn,0 δ ωm+ωn−νm,0 . Hence, 2 Here for later use we consider a generic even power 2p. It is easy to check that the error one makes in considering the continuum limit is of order Now we again make use of the assumption that slow modes and fast modes of physical processes under study are separated by a large gap. Under such assumption we can safely expand the second fraction in the latter expression in power series of slow modes ν n and rewrite (48) as higher-time-derivative expansion. Let us reintroduce the integral over time as 1 = 1 t t 0 dτ νm e i(νn+νm)τ so that powers of ν n can be traded with time-derivative of the potential (note that odd powers vanish upon symmetric sum; in fact they would give rise to total time-derivative terms that are zero upon integration thanks to periodicity in time.) We are thus left with Sums over hard frequencies can be performed in the continuum limit with the help of formula (44) and timederivatives can be partially integrated in order to rewrite the latter in a more symmetric form The infinite higher-derivative expansion is the legacy of non-locality in time: such an expansion can diagrammatically represented as an infinite sum of local (daisy-like) diagrams, as depicted in Fig. 3.
It is intuitive to expect that, in the presence of decoupling low and high frequency modes, the higher-derivative terms should be suppressed. In the next section, we shall generalize this statement and present a quantitative method to systematically organize all contributions to the effective action in terms of a perturbative series.

V. SLOW-MODE PERTURBATION THEORY
The diagrammatic representation of the path integral given by Eq. (38) is formally exact, but rather useless. In fact, it is obviously impossible to evaluate and re-sum exactly all the infinitely many Feynmann graphs appearing in the exponent. On the other hand, in this section we show that it is possible to compute the renormalized effective action S ef f [x < (t)]] to an arbitrary level of precision, by calculating only a finite number of Feynmann graphs. This way, the low-frequency effective theory retains predictive power.
The idea is to exploit the decoupling of the short-time dynamics from the long-time dynamics to organize the sum over all possible graphs as a perturbative expansion. We shall refer to such a systematic evaluation of the renormalized low-frequency effective action as to the slow-mode perturbation theory.
The first step in the construction of our perturbation series is to identify all the dimensionless combinations of the physical quantities which appear in the Feynmann graphs contributing to (38), evaluated in stationary phase approximation. Let us first define the quantities where k is the typical wave vector on the spatial Fourier transform of V ef f (x) and ω < is the typical frequency in temporal Fourier transform of V ef f (x < (τ )). Using these combinations, we can thus construct the following dimensionless combinations: We are interested in describing the dynamics of physical systems for which each of these parameters can be considered small. In order to illustrate the physical interpretation of the condition α 1 1, we observe that the probability for the system to remain in the same configuration x, during an elementary time interval dt is Hence, the combination βV ef f represents 3 the typical time scale associated to local conformational changes, and the condition α 1 1 expresses the condition that the time spent on average by the system in each configuration is large compared to the elementary short-time scale, dt ∼ 1 bΩ . The condition α 2 1 implies that the effective potential varies over length scales which are large, compared with the mean distance covered by Brownian motion in an elementary time interval dt. Finally, the condition α 3 1 implies that the typical slow mode frequencies are small compared to the ultra-violet cut-off, which is of the order of the inverse of the elementary time interval dt.
It is easy to see that any local diagram in the expansion of the renormalized effective action comes about with integer powers of these coefficients, when compared to the tree level effective action. In particular, any diagram composed by r vertices of M 1 , . . . , M r hard fields will involve M = r i=1 M i spatial derivatives and M 2 propagators each of which yields a power of 1 bΩ . Finally, each additional vertex yields a power of βV and each time derivative yields a power of ω. So, the above diagram, at the lowest level in time derivatives will be of order α r−1 1 α M 2 2 with respect to the tree level expression. Higher time derivative terms will add powers α 3 . It is thus natural to define a degree of slowness L for a local diagram, given by where N v is the number of vertices, N τ the number of time derivatives and N x the number of spatial derivatives. The definition in Eq. (54) is normalized in such a way that L(tree) = 0. It is easy to check that the degree of slowness L corresponds to the power of 1 bΩ of the local diagrams. Note also that for daisy diagrams and for all other diagrams where N τ = 0, the degree L is nothing but the number of loops. One can thus easily write down and compute the finite set of local diagrams that renormalize the effective action up to a fixed (yet arbitrary) level of precision L max . Let us consider a few simple examples.
• L ≤ 1 corresponds to a single daisy diagram with L = 1, N v = 1 and N x = 2, represented in the left panel of Fig.4. The expression of this diagram is given by Eq. (45) and gives a correction to the effective action of the form • L ≤ 2 corresponds to two further diagrams, one daisy diagrams with either N x = 4 and the two-vertex local diagram with N x = 2 and no time-derivatives. This latter however is 1PR and gives no contribution. We are thus left with the corrections • L ≤ 3 corresponds to two further diagrams, one daisy diagrams with N x = 6 and the two-vertex local diagram with N x = 4 and no time-derivatives. This latter can be simply read off from eq. (50). Hence that involves in the last term the trace of the square of Hessian of the tree level potential (∂ i ∂ j V (x < )) 2 = Tr H 2 V . In order to see the first time derivatives appearing into the renormalized effective action we need to consider L ≤ 5 where, along with several other corrections, we have the correction coming from the second term in (50) that yields − 3βD 2 5π that can be also recast as a correction of the kinetic action Some comments on the results obtained in this section are in order. First of all, we emphasize that the effective interactions have been derived under the assumption that the modes which are relevant for the long-time dynamics vary over time scales much longer than that of the fast modes, which enter in the loop diagrams. This is the crucial assumption of all renormalization group approaches. Our results confirm the intuitive picture that if one adopts a low "time-resolution power", then the effective interactions generated by the ultra-violet modes can be regarded as instantaneous. This is in fact general property of renormalization group theory, which is preserved to any order in the perturbative expansion. Finally, we note that the correction terms generated by the integration over the fast modes is suppressed, in the small temperature limit.

VI. RENORMALIZATION GROUP IMPROVED MONTE CARLO
The usefulness of the renormalization procedure resides in the fact that it gives rise to an effective theory, in which the largest frequency scale is lowered form Ω to bΩ. Equivalently, the shortest time scale is increased form ∆t to 1 b ∆t. By construction, in the regime of decoupling of fast and slow modes, the probability density generated by the new slow-mode effective theory must be the same as that of the original (i.e. tree-level) theory. In this section, we show how it is possible to use the slow-mode effective theory to develop improved MC algorithms for the time evolution of the probability density P (x, t), in which the elementary time step used to propagate the configurations is increased by a factor 1/b.
The starting point of the MC approach 13 is to write the probability of observing the system in configuration x at time t in terms of the Green's function of the Fokker-Planck Eq. G(x, t|x i , t i ): where ρ 0 (y) is the density of states at the initial time.
One then uses Trotter's formula to write the transition probability as a sequence of intermediate elementary propagation steps: If a sufficiently large number of intermediate steps N is adopted, then the time steps ∆t = t k+1 − t k can be considered infinitesimal and the (unnormalized) transition probability G(y k+1 , t k+1 |y k , t k ) can be calculated analytically G(y + dy, t + ∆t|y, t) = const. × e "Completing the square" in the first exponent, one finds Now and recalling the definition of the effective potential (8) in the second exponent, this Green's function can be written as: In the MC algorithm, one starts from a set of initial system's configurations, sampled according to he distribution ρ 0 (x i ). Such an ensemble is evolved in time, according to the following procedure. Each configuration is propagated for an elementary time interval ∆t, by sampling from the Gaussian in Eq. (64). Such a configuration is then re-weighted according to the factor The iteration of such a procedure for many consecutive elementary propagations gives rise to a set of diffusive trajectories, called walkers. In the so-called diffusion MC algorithm, the term W is used to replicate or annihilate the walkers. The ensemble of configurations obtained according to this procedure is distributed according to the probability density (60). For the MC algorithm to be efficient, the fluctuations in the statistical weight of the walkers -or, equivalently, in the number of walkers-should remain small, throughout the entire time-evolution. This condition is verified if the factor W(y) is always of order one. Note however that this term tends to enhance (suppress) the weight of configurations in the vicinity of the local minima (maxima) of U (y), where the Laplacian is positive (negative). Hence, if the energy landscape varies very rapidly in space, then the fluctuations in the statistical weights -or in the number of walkerswill in general be large, unless the elementary time step ∆t is chosen very small. This feature represents a limiting factor of MC simulations, which makes the sampling of the probability density at large times very computationally expensive.
Clearly, the elementary propagation time step ∆t is the shortest time scale in the simulation. On the other hand, in the slow-mode effective theory one integrates out the dynamics in the time scale range (∆t, 1/b∆t). Hence, we expect that by taking into account of the corrections associated to the renormalized effective interaction it is possible to perform MC simulations in which the integration time step ∆t is chosen a factor 1/b larger.
In practice, the effective slow-mode theory introduces a correction in the re-weighting -or branching-. To order L = 1 one has Notice that this expression contains a factor of the inverse frequency cut-off 1/Ω in the exponent. Such a term is proportional to the elementary time step ∆t. The corresponding proportionality factor reads 2π only for periodic path integral. For a generic initial value MC one can write in general where the constant κ is to be determined from simulations. Hence, we obtain The unknown constant κ can be determined by matching the results obtained by running a short simulation in the tree-level theory -i.e. using an integration step ∆t and the tree-level weighting term (66)-with those obtained in the effective theory -i.e. using an integration step 1/b ∆t and the renormalized weighting term (69)-. In the regime of decoupling of fast and slow modes, once the matching has been done, the two algorithms must generate the same evolution for the probability density at any later times.
In the next session, we shall provide an example which illustrates how this procedure works in practice and show that the fundamental and the effective theory do indeed generate the same long-time stochastic dynamics.

VII. AN ILLUSTRATIVE EXAMPLE
In order to illustrate how the renormalization of the effective interaction works in a simple example, let us consider the dynamics of a point particle, diffusing in a rugged asymmetric harmonic oscillator: with h 1 = 2, h 2 = 1, h 3 = 1, w = 4. The viscosity coefficient is set to γ = 5 and inverse temperature to β = 5. Note that this potential has been chosen in such a way that the average value of x at thermal equilibrium is non-vanishing. The diffusion Monte Carlo algorithm used in our numerical simulations is presented in the appendix B. The factor Ω, which appears in the L = 1, L = 2 improvement terms was determined from the time interval ∆t using Eq. (68). The proportionality constant κ in Eq. (69) was determined once and for all, by matching the result of x(t) of the unimproved (i.e. L = 0) simulations after 10 integration time steps with ∆t = 0.01, with those of the RG-improved (i.e. L = 1, L = 2) MC simulations after a single elementary time step, with ∆t = 0.1. We found κ = 0.35, with no appreciable difference between the L = 1 and L = 2 estimates.
Let us now discuss the results of our simulations. We begin by analyzing the effects of accounting for the factor W(x) defined in Eq. (66), in numerical MC simulations. Fig. 5 shows the average position, once the system has attained thermal equilibrium, obtained by diffusion MC simulations with and without branching the walkers according to W(x). We recall that neglecting such a term is equivalent to simulating the dynamics in the Ito calculus, while the branching is expected to improve the time discretization to order ∆t 2 . Indeed, our results show that, when one chooses small discretization steps, the two approaches are consistent with each other and yield the exact equilibrium average -which was computed directly from the Boltzmann distribution-. On the other hand, at large discretization steps, accounting for the factor W significantly improves the result. The same discussion can be trivially repeated in simulations in which the factor W(x) is interpreted as a re-weighting term, while the number of walkers is held constant.
We now discuss the use of our effective theory to simulate the stochastic dynamics, using large time steps. Fig. 6 shows the time evolution of the average particle position at time t, computed using a small discretization time step -∆t = 0.01-and a large discretization step -∆t = 0.1-. The two curves obtained in the original -i.e. treelevel-theory are compared with the results of the effective theory at order L = 1 and L = 2, which were obtained using an integration time step which was one order of magnitude larger, ∆ t = 0.1.
The time evolution of the observable x(t) , obtained in the tree-level theory using large integration time steps (squares) is inconsistent with the same quantity obtained using small time steps ∆t = 0.01 (circles). This is expected, because for ∆t = 0.1 the numerical simulations of the tree-level theory start to be affected by significant discretization errors -see The results of simulations with large discretization time steps are significantly improved if one uses the effective theory, already at order L = 1 (diamonds). At order L = 2 the dynamics of the tree-level theory simulated at ∆t = 0.01 is indistinguishable from the dynamics of the effective theory simulated with ∆t = 0.1 (triangles). These results show that the hard-mode dynamics in the short time range from 0.01 to 0.1 has been correctly taken into account by means of the renormalized effective interaction. As a consequence, the use of the effective theory allows to obtain very accurate predictions, using larger time steps.

VIII. LONG TIME DYNAMICS OF MOLECULAR SYSTEMS
The improvement of the MC algorithm based on our effective theory is expected to be most efficient when the gap between the slow and the fast modes is very large. In fact, in this regime, the slow-mode perturbation theory remains reliable even when one integrates out a large frequency shell, i.e. when b 1. Hence, in this case, by RG-improvement it is possible to simulate the time evolution using elementary time steps ∆t which are significantly larger than the original elementary time step ∆t, which would be used in the usual (unimproved) MC algorithm.
A natural application of the RG-improved MC is the investigation of the long-time dynamics of macromolecules, for which standard MD or MC algorithms can be extremely computationally expensive. Hence, it is interesting to address the question of what is the typical range of reliability of the slow-mode perturbation theory for a typical molecular interaction, at room temperature. To this end, let us consider the over-damped diffusion at temperature 300 K of two molecules of mass m = 30 amu, interacting through a Van-Der-Waals potential: where = 4 KJ/mol and σ = 0.3 nm. A typical value for the viscosity coefficient for a molecule in its solvent (e.g. an amino acid in water) is γ ∼ 2 × 10 3 amu ps −1 .The typical time-steps used in the numerical integration of the Langevin Eq. (2) are of the order ∆t 10 −3 − 10 −2 ps. The tree-level effective interaction associated to the potential (71) is: This function and the corresponding L = 1 and L = 2 renormalized effective interactions are plotted in Fig. 7 for κ = 1. This plot shows that, for a realistic set of parameters, the perturbative expansion remains reliable even when one integrates out a very large shell of modes, with b ∼ 10 −2 . This fact suggests that the ultra-violet dynamics is essentially free brownian motion, while the long time dynamics is dominated by very low-frequency modes, and is driven by the force field. This fact has remarkable consequences on practical numerical simulations. It implies that by using the renormalized effective potential, it should be possible to adopt integration time steps which are about 10 2 time larger than those required to simulate the dynamics in the original tree-level theory.

IX. CONCLUSIONS
In this work, we have presented a new approach to the problem of investigating the long-time out-of-equilibrium dynamics of multi-dimensional systems obeying Langevin dynamics. In the presence of decoupling of time scales, the methods based on the direct integration of the Langevin Eq. (MD) or on the time propagation of the Fokker-Planck probability density (MC) are usually inefficient, because a significant amount of computational time in invested to simulate uninteresting fast stochastic fluctuations.
We have shown that the decoupling of time scales which limits MD and MC approach can in fact be exploited to perform analytically the average over the short-time stochastic fluctuations. After the integration over the fast modes has been performed, one obtains an effective theory which describes directly the relevant dynamics, with a lower time resolution. In such an effective theory, the effective action in the path integral receives corrections, which account for the ultra-violet physics which is cut-off. We have developed a rigorous scheme which allows to organize such corrections in term of a perturbative series in which the expansion parameters are the ratio between the soft frequency scales and the hard frequency scale bΩ. Hence, sub-leading terms in the perturbative expansion come with higher inverse powers of the hard scales bΩ and become irrelevant in the limit in which the decoupling of fast and slow modes is very large.
The Feynmann diagrams which have to be calculated to obtain the corrections to any given order in this perturbation theory can be identified from their degree of slowness Diagrams with degree of slowness L generate corrections proportional to 1/(bΩ) L . In particular, we have found that the leading-order correction (i.e. L = 1) is proportional to the Laplacian of the effective potential V ef f : At the next-to-leading order, a term containing fourth-order derivatives appears: dτ ∆V ef f (x < (τ )) + 1 2 On the other hand, a space-dependent, tensor correction to the diffusion coefficient appears only as a higher-order effect (L = 5). It is important to stress the fact that, in the present approach, the ultraviolet cut-off Ω (or, equivalently, the short time scale ∆t) is kept finite at all stages. Upon taking the continuum limit ∆t → 0, all the correction terms in the effective theory vanish and one recovers the original theory, defined by the effective Schrödinger Eq. (6). The main usefulness of such an effective theory resides in the fact that it can be used to develop an improved MC approach, to compute the long-time evolution of the Fokker-Planck probability. The elementary time steps used in the RG improved MC algorithm are a factor 1/b larger those of the MC algorithm for the underlying tree-level theory. Since the dynamics in the time range (∆t, 1/b ∆t) is averaged analytically, the RG improved MC algorithm avoids investing computational time in simulating the fast-mode dynamics associated to local Brownian motion.
In the specific case of molecular interactions at room temperature, we have shown that the perturbative approach remains reliable even when one integrates large frequency shells, with b 0.01. This feature suggests that, by using the effective theory, it is possible to simulate time intervals which can be up to a factor ∼ 100 longer than in the usual MC approach.

APPENDIX A: PROPAGATOR OF THE FAST MODES
Here we derive free fast mode propagator G 0 > (ω n , ω m ) appearing in the diagrams, using the standard source technique. We first add a source term to Z 0 > : Then, we functionally differentiate twice with respect to the source: (A3) Note that since the zero mode belongs to the slow modes part of the kinetic action, the kinetic operator for the fast modes is never singular and can be inverted without troubles.

APPENDIX B: DIFFUSION MONTE CARLO ALGORITHM
Our numerical study were performed using the following diffusion Monte Carlo algorithm: