{'fraction_non_alphanumeric': 0.05007479268629188, 'fraction_numerical': 0.018414975964680426, 'mean_word_length': 4.081021251475797, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 39, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Error estimates on ergodic properties of discretized Feynman-Kac semigroups

We consider the numerical analysis of the time discretization of Feynman-Kac semigroups associated with diffusion processes. These semigroups naturally appear in several fields, such as large deviation theory, Diffusion Monte Carlo or non-linear filtering. We present errors estimates a la Talay-Tubaro on their invariant measures when the underlying continuous stochastic differential equation is discretized; as well as on the leading eigenvalue of the generator of the dynamics, which corresponds to the rate of creation of probability. This provides criteria to construct efficient integration schemes of Feynman-Kac dynamics, as well as a mathematical justification of numerical results already observed in the Diffusion Monte Carlo community. Our analysis is illustrated by numerical simulations.


Introduction
The study of Feynman-Kac semigroups for stochastic differential equations (SDEs) has been a topic of growing importance in the past two decades, since these dynamics are related to several theoretical and applied areas of mathematics. They can be seen as standard SDEs whose paths are reweighted according to the exponential of the time integral of some weight function.
Feynman-Kac semigroups naturally appear in large deviation theory, where they can be used to enhance the likelihood of observing rare fluctuations and henceforth computing cumulant generating functions [68,15]. They also have important practical applications, such as in the Diffusion Monte Carlo (DMC) method [26], which is a probabilistic way of estimating the ground state energy of Schrödinger operators; or in computational statistics, in particular in (non-linear) filtering [10,19], where relevant trajectories are selected from observations. We focus in this paper on the bias arising from the time discretization of the underlying continuous stochastic dynamics and of the time integrated weight. Our interest resides in the ergodic properties of the discretization, namely the invariant measure as well as the average rate of creation of probability. Let us briefly present our setting and results. We study a system qt ∈ D evolving in a d-dimensional space, assumed to be compact (the extension to unbounded spaces poses non-trivial issues, as discussed at various places later on). Typically, D = T d (with T = R\Z) is a d-dimensional torus. For convenience, we consider that it satisfies a stochastic differential equation with additive noise: where b : D → R d is a C ∞ (D) vector field, σ > 0 and Bt is a standard d-dimensional Brownian motion. Note that the dynamics (1) may be non-reversible. Our results can be extended to dynamics with multiplicative noise upon appropriate modifications. The infinitesimal generator of the dynamics (1), defined on the core S = C ∞ (D), reads and we denote by L † the adjoint of L on L 2 (D) endowed with the Lebesgue measure. Since D is compact and b is smooth, (1) admits a unique invariant distribution, denoted by ν, which is solution to the Fokker-Planck equation see e.g. [56,37,44]. Denoting by P(D) the set of probability measures over D, Feynman-Kac type semigroups associated with a given weight function W : D → R evolve an initial probability measure µ ∈ P(D) as follows: for any test function ϕ ∈ S, The family of mappings {Φ W t } t 0 is a measure-valued non-linear semigroup in the sense that Φ W t : P(D) → P(D) depends non-linearly on the initial condition and, for all µ ∈ P(D) and t, s ∈ R+, Φ W t (Φ W s (µ)) = Φ W t+s (µ). Such semigroups have been studied for a long time in the context of Diffusion Monte Carlo (DMC) [29,4,7,67,26] in order to estimate the principal eigenvalue of Schrödinger type operators −∆ + W , which correspond in our case to b ≡ 0. They also appear in the large deviation community [28,42,62,65,52,53] where they are related to the principal eigenvalue of L + W , which is the dual of the rate function -a result known as the Donsker-Varadhan formula [18,68,15,65,14]. Other fields such as non-linear filtering, Hidden Markov Models [36,19,20] and free energy computation [35,34,59,45] also motivate the study of such semigroups.
As discussed in Section 2.1, the semigroup (3) converges in general to the average of ϕ with respect to a tilted measure νW . More precisely, the operator L † + W has a largest eigenvalue λ which is isolated from the remainder of the spectrum and non-degenerate, with associated eigenfunction νW , and Under mild assumptions on the discretization scheme (made precise in Section 2.2), the discrete semigroup (4) converges to an invariant measure ν W,∆t in the following sense: for any test function ϕ ∈ S, The core of our work consists in making precise the difference between νW and ν W,∆t . We aim in particular at designing numerical schemes leading to the smallest possible biases. Although a series of papers study the statistical error of estimators such as (4) (see [11,13,10,57,58]), there are, to our knowledge, no available estimates on the bias of the limiting measure with respect to ∆t. However, in the context of DMC (where we recall b = 0), it was numerically observed that some discretizations provide first, second or fourth order of convergence in ∆t for the largest eigenvalue λ of L + W , see for example [4,67,49,60], and [23] for the numerical analysis in a simple case. The results presented in this paper provide a mathematical justification of such convergences, while extending them to the case b = 0. Let us also mention that Hairer and Weare have studied in [32,33] the convergence with respect to the timestep of discretized dynamics similar to the one we consider, over a finite time and for a finite population of replicas. They obtain in the limit ∆t → 0 a limiting process, the so-called Brownian fan.
We rely on the techniques developped since the works of Talay and Tubaro [64,63], taking advantage of the analytical tools developed in a series of papers [48,9,1,2,6,43,44], in order to provide a systematic framework to study the bias in the timestep. More precisely, we show in Theorem 2 that there exist an integer p 1 and a function f solution to a Poisson equation (both depending on the numerical scheme at hand and the quadrature rule for the integral), such that, for all ϕ ∈ S, D ϕ dν W,∆t = D ϕ dνW + ∆t p D ϕf dνW + ‰(∆t p+1 ). (5) This result is very similar to those of weak backward error analysis on invariant probability measures of ergodic processes, see for example [9,6,43]. Moreover, as the computation of the principal eigenvalue λ of the operator L + W is one of the main concerns in Feynman-Kac techniques, we provide in Theorem 3 the following error estimate: where Q W ∆t is the evolution operator of the discretized dynamics with weight function W . This result is interesting since it allows to justify the use of population dynamics methods for discretizations of diffusion processes, see [28,62,52] for rare events simulations, and [26] for DMC. Let us mention that, while the proof of (5) relies on previous works concerning error estimates on the invariant measure [64,63,6,9,43], the novelty of this work lies in taking into account the non-probability conserving feature of the dynamics. With this point of view and, odd as it may seem, formula (5) appears as a consequence of (6), and not conversely. An interpretation of this fact is that, in order to prove an error estimate on the invariant probability measure of this non probability-conserving dynamics, we must first show that the discretized process creates or destroys probability at a rate correct up to terms small in ∆t.
The paper is organized as follows. Section 2 is devoted to general properties of Feynman-Kac semigroups and their discretizations. We then present in Section 3 our main results concerning the numerical analysis of the error on the invariant probability measure, depending on the choice of the discretization scheme, before providing numerical applications in Section 4. Finally, Section 5 proposes possible extensions to this work. The proofs of the most technical results are gathered in Section 6.

Convergence properties of Feynman-Kac semigroups
We present in this section the setting of our study. In particular, we remind convergence results and some useful properties of continuous Feynman-Kac semigroups in Section 2.1, as well as convergence results for their discretizations in Section 2.2. Although these results are known, we believe that it is useful to gather them here to allow for a self-contained presentation of the numerical analysis framework developped in Section 3.

Continuous dynamics
We denote by Pt the evolution operator associated with the process (qt) t 0 in (1): for all µ ∈ P(D) and ϕ ∈ S,

Its weighted counterpart is
The infinitesimal generators of Pt and P W t are respectively L and L + W . We assume in the sequel that W is smooth.
The existence of a spectral gap for the generator L + W and its adjoint is a key ingredient for our study. Here and in the sequel, and otherwise explicitly mentioned, all operators are considered on the Hilbert space For a given closed operator T on L 2 (ν), we denote by T * the adjoint of T in L 2 (ν). In particular, In this functional framework, the reversibility of the dynamics is equivalent to the self-adjointness of L on L 2 (ν). We however do not assume that this is the case, and this is why we need to distinguish between eigenelements of L and L * . We can then state the following. Proposition 1. The operator L+W , considered on L 2 (ν), has a real isolated principal eigenvalue λ with associated eigenfunctionĥW ∈ S normalized as The operator L * + W then also admits λ as a real isolated principal eigenvalue, with associated eigenfunction hW ∈ S normalized as Moreover, the functionsĥW and hW are positive.
The fact thatĥW , hW ∈ S is a consequence of elliptic regularity. Let us emphasize that, as a consequence of (8), the measure νW = hW ν is the only invariant probability measure for the evolution encoded by P W −λ t . Moreover, when the underlying diffusion is reversible, i.e. b = −∇V and ν(dq) = Z −1 e −2V (q)/σ 2 dq, the operator L is self-adjoint (L * = L) so thatĥW = hW . When W = 0, it simply holds hW = 1 whatever b is.
Proof. It is shown in [27] that the operator L + W has a real isolated principal eigenvalue when considered as an operator on C 0 (D), the space of continuous functions on D. This can be proved using the Krein-Rutman theorem [22]. On the other hand, standard results of spectral theory of elliptic operators on bounded domains show that L + W on L 2 (ν) has a discrete spectrum, which is bounded above [55]. The first eigenvalue cannot be degenerate since the associated eigenvectors are smooth by elliptic regularity and are therefore also eigenvectors of L + W considered as an operator on C 0 (D). Finally, the positivity ofĥW and hW follows from the fact that the evolution semigroup P W t and its adjoint are operators with smooth and positive transition kernels (since the noise is non-degenerate), together with the equalities P W tĥW = e λtĥ W and (P W t ) * hW = e λt hW .
In what follows, we use the subspaces L 2 W (ν) and SW of functions of average 0 with respect to νW : We also introduce the measureνW =ĥW ν, the spacê and we denote by the spectral gap of L + W in L 2 (ν). The fact that the largest eigenvalue λ is a priori non-zero corresponds to a possible creation (λ > 0) or destruction (λ < 0) of probability induced by the source term W , which plays the role of an importance sampling function. The statement about the spectral gap in Proposition 1 implies the convergence of the Feynman-Kac semigroup (3), as stated in the following result.

Proposition 2.
There exists C > 0 such that, for all µ ∈ P(D) and ϕ ∈ L 2 (ν), As made clear in the proof of this result (see Section 6.1), it is possible to consider any observable ϕ ∈ L 2 (ν) even if µ is singular. This is due to the regularizing properties of the underlying diffusion for positive times, and explains why the convergence result is stated only for times t 1. The next proposition will be frequently used in this work.

Proposition 3. It holds
Proof. Integrating both sides of (8) on D, A natural corollary of Propositions 2 and 3 is that the largest eigenvalue of L + W can be obtained by a long time average of W using the Feynman-Kac semigroup (3).

Corollary 1.
There exists C > 0 such that, for any initial distribution µ ∈ P(D), Another important consequence of Proposition 1 is the invertibility of the generator and its adjoint over suitable functional spaces. Proposition 4. The operator L + W − λ is invertible on SW , in the sense that, for any g ∈ SW , the Poisson equation The proof of this result can be read in Section 6.1. Let us emphasize that the smoothness of W is crucial for this proposition to be true. Note also that the stability of the core of the operator L + W would be harder to prove for non-compact state spaces, as this is already a non-trivial statement for the Poisson equation with W = 0, see [40,41].

Discretization
We now turn to the discretization of the Feynman-Kac semigroup (3). We first define discretization schemes, and show that they are ergodic for some limiting measure under mild assumptions. We also recall the stationarity equation satisfied by this invariant probability measure, which proves crucial for the numerical analysis developped in Section 3.
The properties of discretized Feynman-Kac semigroups are related to the properties of the underlying discrete dynamics. The approximation of the continuous dynamics (1) is given, for a time time ∆t, by a Markov chain (q n ) n∈N such that q n ≃ q n∆t . This Markov chain is characterized by the evolution operator Q ∆t defined as A typical example is the Euler Maruyama scheme defined by: where (G n ) n 0 is a familly of independent and identically distributed standard d-dimensional Gaussian random variables. In order to perform our analysis in Section 3, it is convenient to rephrase discretizations of (3) such as (4) in terms of an evolution operator. For instance, we see that, defining the discretization (4) reads, for an initial measure µ0 and a test function ϕ, We use the same definition for more general discretizations of (3) characterized by an evolution operator Q W ∆t . Consistency requirements on Q W ∆t are made precise in Assumption 2 below. This allows us to take into account various integration rules, both for the underlying dynamics and the exponential weights. For instance, the choice well-known in the diffusion Monte Carlo community [60,49,46,67], defines the following semigroup: In what follows, given that the discrete semigroup defines a measure-valued dynamics, we write for simplicity µn = Φ W ∆t,n (µ0). For a given bounded operator Q on L ∞ (D) and a probability measure µ ∈ P(D), we also denote by µQ the probability measure defined as We start by recalling a one-step formulation of the non-linear dynamics (µn) n 0 , as suggested e.g. in [13]. This formulation is the basis for a stationarity property fundamental in our numerical analysis.

Lemma 1.
The sequence of probability measures µn = Φ W ∆t,n (µ0) satisfies the following dynamics: Proof. The proof relies on a simple rewriting: for all ϕ ∈ S, which gives the result.
Let us now prove that the measure-valued dynamical process (17) admits a limit measure µ∞ independent of the initial distribution µ0, and that the long time average (14) converges to the average with respect to this measure. We follow the strategy of Del Moral and collaborators [11,13,12,10], which relies on the Dobrushin ergodic coefficient of a relevant operator (see Appendix A). For this, we use the following assumption, which is typically satisfied for discretizations associated with the continuous dynamics (3) on the torus.

Assumption 1. The operator Q W
∆t satisfies a minorization and boundedness condition: there exist ε ∈ (0, 1) and η ∈ P(D) such that, for all non-negative bounded measurable function ϕ, The condition (18) is satisfied for the evolution operator (13) as soon as a condition similar to (18) is satisfied for the evolution operator Q ∆t . The latter condition is, in turn, easily seen to be true for the numerical scheme (12), with η(dq) = |D| −1 dq the normalized Lebesgue measure on D, see [44,Section 3.3.2]. Similar considerations allow to prove that (18) holds for more complicated discretization strategies [31,44].

Remark 1.
Let us however emphasize that the prefactor ε in (18) typically scales as ∆t − d 2 exp(−CL/∆t) for some constant CL > 0, as it is proportional to the probability density of reaching in one step any configuration in D starting from any other one.
We can now recall an important result which ensures the existence and uniquess of the limiting measure for the discretized Feynman-Kac dynamics. Its proof, taken from [11], is recalled in Section 6.2. To state the result, we introduce the total variation distance between two measures µ, ν ∈ P(D): where the supremum runs over measurable subsets of D. Recall that P(D) is complete for this distance. Theorem 1. Suppose that Assumption 1 holds true. Then the non-linear dynamics (17) admits a unique stationary probability measure µ∞ which is independent of the initial measure and which is a fixed point of K: Moreover, for any initial distribution µ0 ∈ P(D), Remark 2. We expect this result to hold true for a larger class of dynamics, and with more quantitative estimates. These perspectives are discussed in Section 5.
As a consequence of Theorem 1, if we define a discretization of the Feynman-Kac semigroup (3) satisfying Assumption 1, the discrete dynamics (14) admits an invariant probability measure solution to the fixed point equation (19). We denote by ν W,∆t this invariant probability measure to emphasize its dependence on both W and the timestep ∆t. In view of (19) and (17), this measure satisfies the following stationarity equation: In particular, if we define the approximate eigenvalue λ ∆t by (21) can be rewritten as: This is the stationarity equation of the discretized process upon which the analysis in Section 3 is built. Let us emphasize that it involves the approximate eigenvalue λ ∆t accounting for the rate of creation of probability of the discretized process, which differs in general from the largest eigenvalue λ of the operator L + W (which accounts for the rate of creation of probability for the continuous process). The numerical analysis of the approximation λ ∆t of λ plays an important role in Section 3.

Remark 3.
In the case W ≡ 0, the measure ν W,∆t = ν ∆t is the invariant probability measure of the discretized process without reweighting, and the evolution operator Q ∆t conserves probability. This also implies that λ ∆t = 0. Therefore (23) simplifies as which is the standard stationarity equation of the invariant probability measure for discretizations of SDEs [43,44]. Note that the largest eigenvalue of the discretized evolution operator Q ∆t is 1, as for the continuous semigroup Pt.

Numerical analysis of the discretization
We now turn to the main section of the paper, where we quantify how close ν W,∆t , the ergodic measure for the discrete Feynman-Kac dynamics, is from νW , the ergodic measure for its continous counterpart. We also make precise the difference at leading order. Following a general strategy to study the error on the invariant probability measure of discretizations of stochastic processes dating back to [64] (see also [9,43] as well as the review [44] for recent accounts), we compare the evolution operator Q W ∆t with the Feynman-Kac semigroup e ∆t(L+W ) . Although the non probability-conserving feature of the dynamics is an additional difficulty, we obtain in Section 3.1 results similar to those of [64,9,2,44] concerning the error on the invariant probability measure. Moreover, we provide in Section 3.2 error bounds for estimators of the eigenvalue λ. Finally, we show how to relate the invariant probability measures of different schemes in Section 3.3 and discuss in Section 3.4 how the Feynman-Kac discretization essentially inherits the properties of the discretization of the underlying unweighted dynamics.

Expansions of the discrete evolution operators Q W ∆t
For unweighted dynamics (W = 0), consistency assumptions on the evolution operator Q ∆t characterizing the discretization rely on an expansion in powers of ∆t of Q ∆t (see the presentation in [44]). More precisely, it is assumed that there exist an integer p 1 and differential operators (A k ) k=1,...,p+1 such that the evolution operator Q ∆t of the discrete dynamics admits the following expansion: for all ϕ ∈ S, The differential operators A k have finite order and smooth coefficients: for any k ∈ {1, . . . , p + 1}, there exist m k ∈ N and a familly of smooth functions (aα) |α| m k (with α = (α1, . . . , α d ) ∈ N d ) such that where Moreover, R ∆t is an operator uniformly bounded in ∆t in the following sense: there exist ∆t * > 0, c > 0 and m ∈ N such that where The assumptions (24) and (26) are standard for the numerical analysis of ergodic measures of SDEs [64,40,9,1,43,44], and are satisfied for a wide range of explicit and implicit schemes defined on compact domains. A scheme is of weak order p when (24) holds with see for instance [51]. Typically, A1 = L for any reasonable discretization scheme. Besides weak and strong errors, another notion of consistency is the error arising on the invariant probability measure, in situations when the Markov chain associated with Q ∆t admits an invariant probability measure ν ∆t . The error between averages with respect to ν and ν ∆t are of order at least ∆t p when the scheme is weakly consistent of order p. It can however be of higher order ∆t p ′ (with p ′ p + 1) when This condition is satisfied by operators which are proportional to powers of L. See however [3,43] for examples of situations where A k is not a power of L but the above condition is met.
In the context of Feynman-Kac averages (3) where we consider approximations Q W ∆t of e ∆t(L+W ) , we generalize the conditions (24) and (26) where R W,∆t is a uniformly bounded remainder in the sense of (26). We also assume that A W 1 is such that where A1 is a differential operator. In particular, A W 1 1 = W . Let us provide an example of such an expansion when Q W ∆t is defined by (13). Lemma 2. Assume that (24) and (26) hold with A11 = 0, and define Q W ∆t = e ∆tW Q ∆t . Then Assumption 2 holds with, for all k ∈ {1, . . . , p + 1}, Proof. The equality follows by expanding the exponential and taking the product with the semigroup expansion: there exist ∆t * and K > 0 such that with r ∆t,W C 0 K for 0 < ∆t ∆t * . Gathering the terms of order ∆t k leads to (31) plus a uniformly bounded remainder, which proves the result.
Note that, in (31), we obtain A W 1 = A1 + W where A1 is defined in (24). However, there are other ways to construct Feynman-Kac schemes Q W ∆t , using for instance a splitting strategy. Let us give an example. Assume for instance that the operator L can be split in two parts: L = L1 + L2. We can then define a splitting scheme as Q ∆t = e ∆tL 2 e ∆tL 1 , and, by discretizing the time integral of W in three parts (using Simpson's rule) intertwinned with e ∆tL 2 and e ∆tL 1 , In this case, we see that the expansion of Q W ∆t cannot be derived from the one for Q ∆t (by a statement similar to (31)). The evolution operator Q W ∆t nonetheless satisfies Assumption 2.

Statement of the main result
Before stating our main theorem, we need to introduce the following technical assumptions.
In our setting, a crucial step of the proof consists in building an approximation of the eigenvectorĥW to solve an approximate eigenvalue problem for the operator Q W ∆t . This is an important difference compared to the case W ≡ 0, and requires the following assumption.
Assumption 4 (Spectral consistency). The operator A1 + W , considered on L 2 (ν), admits λ as its largest eigenvalue, with associated eigenvectorĥW : Note that Assumptions 3 and 4 are immediately met when the schemes are weakly consistent, i.e. A1 = L, since Assumption 3 is equivalent to Proposition 4 while Assumption 4 follows from Proposition 1. However, it is possible in principle to construct numerical schemes for which A1 = L, in which case Assumptions 3 and 4 should be checked directly.
We are now in position to state our main result on the numerical discretization of Feynman-Kac semigroups, which makes precise error estimates à la Talay-Tubaro in the ergodic setting.

Theorem 2. Suppose that Assumptions 2, 3 and 4 hold. Assume also that the operators
Define also f ∈ SW as Then, there exist a timestep ∆t * > 0 and an operator R W,∆t (uniformly bounded in ∆t in the sense of (26)) such that, for any 0 < ∆t ∆t * and any ϕ ∈ S, Note that the denominator in the second term on the right-hand side of the definition ofg is positive thanks to Proposition 1. In general, in (29), we expect A W k to be (L + W ) k /k! (which corresponds to a scheme of weak order k), in which case (32) holds for a k = λ k /k! = 0 (see (47) below for a proof of the latter equality). This factor comes from the fact that Q W ∆t does not conserve probability. Indeed, for the evolution operator Q ∆t of a Markovian dynamics, one always has On the other hand, considering (29) and applying (32) where r W,∆t is a remainder term which is uniformly bounded for ∆t sufficiently small. This is the manifestation at the discrete level of the fact that, over a timestep ∆t, the dynamics increases or decreases approximately the probability mass by a factor e ∆tλ . The relation (32) should be compared to the invariance relation (28) for W = 0.

Proof of Theorem 2
The proof of Theorem 2 relies on four lemmas which allow to easily conclude the proof. We follow the same strategy as for the error analysis of the invariant probability measure proposed in [43,44] but additionnal technical difficulties arise due to the non-linearity of the stationarity equation (21). The first step (Lemma 3) is to construct the leading correction term f . We next use a projector in Lemma 4 to relate the exact stationary measure νW and its approximation ν W,∆t . An a priori estimate on the approximate eigenvalue defined in (22) is then provided in Lemma 5. Finally, an approximate inverse operator is constructed in Lemma 6. In the proofs and also in the statements below, the remainders may change from line to line in the calculation, but we do not change the notation for convenience. There are two types of remainders: terms of the form R W,∆t ϕ where R W,∆t is a differential operator satisfying (26), and functions r W,∆t such that, for any k 1, there is K > 0 and ∆t * for which r W,∆t C k K when 0 < ∆t ∆t * .
To begin with, we give the expression of the leading correction term f . It relies on an approximate reformulation of (21) which leads to an expression similar to (23) up to a remainder of order ∆t p+1 .

Lemma 3. Under the assumptions of Theorem 2, for any
where R W,∆t is a uniformly bounded remainder (in the sense of (26)) and f is defined in (33).
The proof of this lemma is presented in Section 6.3.1. Defining the approximate eigenvalueλ ∆t by (35) can be rewritten as This expression allows to identify the leading order correction term ∆t p f in ν W,∆t − νW and can be thought of as the approximate counterpart of (23). The second step is to use a projector that on the one hand stabilizes in SW the operator appearing in (37), and on the other hand relates the exact stationary measure νW and its approximation ν W,∆t . For this we introduce the following projectors: for all φ ∈ S, The operator Π is the L 2 (ν) orthogonal projector on L 2 0 (ν), while ΠW is a projector on L 2 W (ν) which is not orthogonal for the canonical scalar product on L 2 (ν). However, it is orthogonal on L 2 (νW ), so that, for all ψ, φ ∈ S, We can then show the following result, whose proof can be found in Section 6.3.2.

Lemma 4. Under the assumptions of Theorem 2, it holds, for any
where R W,∆t is a uniformly bounded remainder in the sense of (26).
Here, we see that two different operators appear inside the integrals because the factors e ∆tλ ∆t and e ∆tλ ∆t are different. The next lemma shows that these quantities are the same up to terms of order ∆t p+2 . As mentioned earlier, this is an important difference with the analysis in the case W = 0. Some a priori estimate on the eigenvalue is required to conclude the proof, whereas, for the unweighted case, the largest eigenvalue of the evolution operator is 1 with eigenvector 1 both for the continuous process and its discretization. The proof, provided in Section 6.3.3, relies on building an approximate eigenfunction for the operator Q W ∆t . Similar estimates were obtained in the Diffusion Monte Carlo context in analytically solvable cases in [49].
Once we have reached this point, it is possible to replace the eigenvalue e ∆tλ ∆t by e ∆tλ ∆t in Lemma 4. The last step is to build an approximate inverse of the operator as provided in the next lemma (see Section 6.3.4 for the proof).

Lemma 6.
Under the assumptions of Theorem 2, for any 0 < ∆t ∆t * , there is an operator S ∆t : S → S for which where R W,∆t : S → S is a uniformly bounded remainder in the sense of (26), and S ∆t admits the following uniform bounds: for any k 0, there exist K > 0 and m ∈ N (depending on k) such that We now have all the tools to prove Theorem 2. First, plugging the estimate (42) obtained in Lemma 5 in the error expansion (40) obtained in Lemma 4 leads to, for any φ ∈ S, where R W,∆t satisfies (26). We next consider the approximate inverse operator S ∆t built in Lemma 6, and set φ = S ∆t ΠW ϕ in (44). Therefore, for any ϕ ∈ S, where R W,∆t satisfies (26). Since f has average 0 with respect to νW , this gives which concludes the proof of Theorem 2.

Alternative error estimate for the principal eigenvalue
We present in this section a useful application of Theorem 2, which provides an error estimate for the approximation of the principal eigenvalue λ of the operator L + W . The choice ϕ ≡ W allows to compute this eigenvalue by ergodic averages, as shown in Proposition 3 and Corollary 1. As a result, this eigenvalue can be approximated using Theorem 2, whose application to ϕ ≡ W gives where r W,∆t is uniformly bounded for ∆t small enough. Although this formula can be used in simulations to estimate λ, we present an error estimate for an alternative approximation more commonly used in practice.
We will also see in Section 3.4 that this alternative formula can be more accurate than the estimate based on averaging W .
Theorem 3. Suppose that Assumption 2 holds, with a numerical scheme consistent at order p (that is, Then there exist ∆t * > 0 and C > 0 such that with |r ∆t,W | C for any 0 < ∆t ∆t * , and This result is important since it implies that we can approximate the eigenvalue λ by computing λ ∆t , which is proportional to the logarithm of the average creation of probability over a timestep ∆t (given by Q W ∆t 1) at stationarity. This is the reason why we need the coefficients a k to be correct up to order p (i.e. a k = λ k /k!) since they represent the creation of probability of the discretized process. The estimate (45) justifies the use of population based dynamics [28,62,52] when the underlying continuous diffusions are discretized in time. We illustrate the error estimate (45) in the numerical simulations reported in Section 4.
Proof. We use Lemma 5 to prove the theorem, which highlights the importance of this result in our context. In all this proof, r W,∆t denotes a smooth function which may change from line to line, but whose C 0 norm is always uniformly bounded for sufficiently small timesteps ∆t. From the definition (22) and the estimate (42), Expanding e ∆tλ ∆t defined in (36) in powers of ∆t and recalling that A W 1 1 = W , where we used that D f dνW = 0 and, in view of (8), Therefore, Given that e ∆tλ is uniformly bounded for 0 < ∆t ∆t * and equal to 1 at leading order in ∆t, we obtain, by expanding the logarithm, The result then follows from e −∆tλ = 1 + ∆t r λ,∆t and the definition (46) of λp+1.

TU Lemma
In the context of splitting schemes, it may be useful to relate the invariant probability measures of two numerical schemes differing by the ordering of the applied operators. This is the purpose of a result called "TU lemma" in [43], which we adapt to our context in Lemma 7. We then state a similar version of this lemma for the eigenvalues of two such schemes in Lemma 5. We will see in Section 3.4 that this last result can be combined with Theorem 3 to show that the schemes (13) and (15) both provide second order estimates of the principal eigenvalue λ using (45), when the discretization of the process Q ∆t is weakly consistent of order 2.
Then, for any ϕ ∈ S, For the TU lemma stated in [43], the typical case of application corresponds to Q ∆t = U ∆t T ∆t and Q ∆t = T ∆t U ∆t , with two Markov operators T ∆t and U ∆t . In this case, the relation (48) holds with a power n − 1 on the right-hand side. For Feynman-Kac semigroups, T W ∆t and U W ∆t are a priori such that T W ∆t 1 = 1 and U W ∆t 1 = 1. A typical case of interest is Q W ∆t = e ∆tW Q ∆t and Q W ∆t = Q ∆t e ∆tW · , in which case (48) is satisfied with T W ∆t = e −∆tW and U W ∆t = e ∆tW .
Proof. For any µ ∈ P(D) and any ϕ ∈ S, The result then follows from the ergodic limits as provided by Theorem 1.
In our framework, the approximate principal eigenvalue λ ∆t is another important feature of a discretization scheme. In fact, under an additional assumption on the operators T W ∆t and U W ∆t , schemes related by (48) share the same approximate eigenvalues in the sense of (22). This is made precise in the following proposition (see Section 6.4 for the proof).
Proposition 5. Fix a timestep ∆t > 0 and consider a numerical scheme for the Feynman-Kac dynamics corresponding to an evolution operator Q W ∆t satisfying Assumption 1, with associated invariant measure ν W,∆t given by Theorem 1, and eigenvalue λ ∆t defined by (22). Consider next a second scheme corresponding to an operator Q W ∆t related to Q W ∆t by (48), with operators U W ∆t and T W ∆t bounded on L ∞ (D) and for which there exists α > 0 such that, for any ϕ ∈ S with ϕ 0, Then, Q W ∆t satisfies Assumption 1, and its invariant probability measure is denoted byν W,∆t . Moreover, its associated eigenvalueλ ∆t defined byλ is such that:λ ∆t = λ ∆t . The eigenvalueλ ∆t should not be mistaken in this context with the definition (36), which serves as an intermediate in the proof of Theorem 2. A careful inspection of the proof shows that it would be possible to consider a slightly different assumption (50).

Remark 4. Although Proposition 5 may look odd at first sight, it has a natural interpretation in terms of matrices.
Indeed, if A ∈ R n×n and B ∈ R n×n are two square matrices with nonnegative entries, the products AB and BA share the same real principal eigenvalue. One can show this by the following argument. For any matrix M ∈ R n×n with nonnegative entries, the spectral radius is an eigenvalue of M (see [61]). This eigenvalue is the equivalent of the principal eigenvalue for the operator L+W This leads to ρ(AB) ρ(BA), and, by symmetry, ρ(BA) ρ(AB); hence ρ(AB) = ρ(BA). In the same way, evolution operators related by (48) share the same principal eigenvalue even though, a priori, they do not admit the same invariant probability measures. The proof of Proposition 5, presented in Section 6.4, follows a path similar to the one used here for matrices.

Second order schemes
We now turn to second order schemes for Feynman-Kac dynamics. They are the most interesting ones in practice, since they can provide an important improvement in the accuracy for a relatively cheap computational overhead. Moreover, in our case, they can be straightforwardly built from second order schemes for the dynamics (1), as a consequence of Theorem 2. (24) and (26) hold with the following expansion for Q ∆t :

Lemma 8. Suppose that
where A3 is a differential operator with smooth coefficients and R ∆t satisfies (26). Then the operator Q W ∆t defined by satisfies Assumption 2 with p = 2: where and R W,∆t satisfies (26).
The interpretation of this result is the following: when we have a scheme consistent at order 2 for the dynamics with W = 0, we immediately obtain a second order scheme for the Feynman-Kac dynamics using the corresponding Markov chain and a trapezoidal rule for the time integral in the exponential. Thanks to the consistency at order one (A W 1 = L + W ) and Propositions 1 and 4, the assumptions of Theorems 2 and 3 are immediately satisfied with p = 2.
Proof. The expression of A W 3 can be obtained by a direct computation or with the Baker-Campbell-Hausdorff formula (see [30]), which is a convenient way to perform the algebra allowing to make precise the various terms in expansions in powers of ∆t. Let us sketch how this is done, and refer to [43] for strategies of proof in order to make the expansions below rigorous. First, When we are interested in the computation of the principal eigenvalue with Theorem 3, we can in fact show that the left-point integration (13) is sufficient for λ ∆t to be correct at order 2 if Q ∆t is consistent at order 2 (i.e. (52) holds). In particular, the discretization scheme for the Feynman-Kac dynamics need not be consistent at order 2 for the eigenvalue to be correct at order 2 (in the same way that the invariant probability measure for discretizations of ergodic SDEs can be correct at order 2 even if the discretization itself is only weakly consistent at order 1, see [3,43,44]). This consequence of Proposition 5 is made precise in the following proposition.
Proof. The proof is a simple application of Proposition 5. Consider the scheme Q W,δ ∆t for a fixed δ ∈ [0, 1] and the scheme Q W ∆t = e ∆t W 2 Q ∆t e ∆t W 2 , which corresponds to a trapezoidal approximation of the integral. We can assume without loss of generality that Q W ∆t satisfies Assumption 1. Then, Q W,δ ∆t is related to Q W ∆t through (48) for the corresponding operators: The operators U W ∆t and T W ∆t are bounded on L ∞ (D) and satisfy (50) with α = e −∆t W L ∞ (D) /2 . Therefore, by Proposition 5, the eigenvalue λ δ ∆t associated to Q W,δ ∆t is equal toλ ∆t , the eigenvalue associated to Q W ∆t , and thus does not depend on δ. Moreover, by Lemma 8, if Q ∆t satisfies (52), Q W ∆t satisfies the assumptions of Theorem 3 with p = 2. This shows that the eigenvalue λ ∆t satisfies (45) with p = 2 whatever the integration rule (i.e. for any δ ∈ [0, 1]).

Numerical application
The goal of this section is to illustrate the error estimates presented in Section 3 on a toy example. For this, we consider (1) over the one dimensional torus D = T with possibly non-gradient drifts: where V is a smooth potential and γ ∈ R. Let us emphasize that a constant force is not the gradient of a smooth periodic function. We first make precise in Section 4.1 the Monte Carlo algorithm used to compute the Feynman-Kac averages. We next describe in Section 4.2 a Galerkin method to compute reference values for the properties of interest. Note that such a discretization method can be used only for low-dimensional systems; but, when it can be used, it typically provides more accurate results than stochastic methods. Finally, we present our numerical results in Section 4.3.

Monte Carlo discretization
Discretization of the underlying SDE. The Euler-Maruyama discretization of the dynamics (54) is given by: where G n are independent and identically distributed one-dimensional standard Gaussian variables. It is well known that this scheme is weakly consistent of order one (see for instance [51,9]. In order to test our results on a second order scheme, we use a discretization proposed e.g. in [1,69,25,66]: It can be proved that this scheme is of weak order 2. Weighted dynamics. Once the underlying SDE has been discretized, a Monte Carlo scheme for approximating the associated Feynman-Kac semigroup (14) has to be devised. Several methods have been succefully applied in order to compute Feynman-Kac averages, generally referred to as Sequential Monte Carlo or Population Monte Carlo methods [21,10,45]. For simplicity and numerical efficiency, we present here a population method with multinomial resampling. More precisions on this familly of algorithms are available in [21], see also [45,Chapter 6] in the context of free energy computation and [32] in the context of Diffusion Monte Carlo. The algorithm relies on a dynamics run over a set of replicas of the system. At each step, the replicas are updated according to the dynamics prescribed by the evolution operator Q ∆t , and are assigned an importance weight depending on the choice of discretization rule for the integral. The replicas are then resampled following a multinomial distribution with their respective weights, before computing the desired averages. This technique prevents the variance of the estimator to increase exponentially in time, a common problem when computing directly quantities such as (14). We now make precise the algorithm.
Consider a population of M replicas (qm)m=1,...,M initially distributed according to some probability measure µ ⊗ M over D M and evolving through a Markov kernel Q ∆t with timestep ∆t > 0. We denote by χ ∆t : D × D → R a weight function to be chosen later on. The algorithm consists in repeating for each time 0 n < Niter the following steps: The principal eigenvalue of the operator L + W is then estimated with (45) through while the average of ϕ is estimated by where the ≈ sign indicates the approximation arising from the finiteness of the number M of replicas and of the number Niter of steps. We do not take these errors into account and ensure numerically that they are sufficiently small in our simulations to observe the bias due to the timestep. The reader interested in the convergence rates of this type of algorithm when M → +∞ and Niter → +∞ is refered e.g. to [13,21,58,23]. Until now, we did not specify the choice of function χ ∆t , which depends on the discretization rule for the integral in (3). In practice, given a discretization of the SDE characterized by an operator Q ∆t , we use the schemes defined by the left point integration e ∆tW Q ∆t , and by the trapezoidal integration e ∆t W 2 Q ∆t e ∆t W 2 . They correspond respectively to the choices:

Galerkin discretization
We now make precise the Galerkin method that can be used to estimate λ ∆t and D ϕ dν W,∆t . This discretization provides reference values for the Monte Carlo method described in Section 4.1. In particular, when V = 0 and γ = 0, the two methods should give the same result since the Euler scheme (55) is exact in law in this specific case.

Choice of the Galerkin basis.
Since we work with periodic function, we consider the Galerkin subspace Span{e−N , . . . , eN } with en(q) = e 2iπnq .
The generator of the SDE (54) reads The The value of N is chosen sufficiently large for all results to be converged with respect to this parameter. The only source of error in the quantities we compute then arises from the finiteness of the timestep ∆t > 0, and possibly numerical quadratures to evaluate certain integrals. Our experience shows that N = 30 is already sufficient for the applications described in Section 4.3.
References quantities for ∆t = 0. The invariant probability measure νW satisfies the eigenvalue problem (L † + W )νW = λνW . We compute a reference approximation λ N 0 to λ by computing the eigenvalue of L N + B N with the largest real part: The associated eigenvector allows to construct the following approximation of νW : The normalization condition [V N W,0 ]0 = 1 ensures that ν N W,0 has a total mass 1. Averages of observables ϕ are then estimated by computing the following integral D ϕ(q)ν N W,0 (q) dq using a one-dimensional quadrature rule.
Reference quantities for ∆t > 0. We next approximate the evolution operators of the first order scheme e ∆tW Q ∆t and of the second order one e ∆t W 2 Q ∆t e ∆t W 2 , respectively as For each value of ∆t, we construct the above matrices, and compute their respective principal eigenvalues Λ N ∆t,1 , Λ N ∆t,2 and eigenvectors V N,1 W,∆t , V N,2 W,∆t ∈ C 2N+1 by diagonalization (still with the normalization condition [V N,j W,∆t ]0 = 1 for j = 1, 2). We then consider the following approximations of the principal eigenvalue λ of the Feynman-Kac operator L + W , based on (45): Averages of ϕ with respect to the invariant probability measure are approximated by the following quantity, using the eigenvectors V N,1 W,∆t and V N,2 W,∆t : for j = 1, 2, In view of Theorem 2, we expect the average of ϕ to converge linearly in ∆t for the first order scheme when ∆t → 0, and quadratically for the second order scheme. We also use the TU-lemma to show that, by appropriately correcting the first order scheme, we recover the same results as for the second order scheme. More precisely, we apply (49) with U W ∆t = e ∆t W 2 , which leads to the following approximation of the average (estimated in practice using a numerical quadrature): On the contrary, from Proposition 6, the eigenvalues λ N ∆t,1 and λ N ∆t,2 should be equal, and therefore λ N ∆t,1 need not be corrected.

Numerical results
Zero-potential case. We first choose V = 0, σ = √ 2, W (q) = (cos(2πq)) 2 and ϕ(q) = exp(cos(2πq)). As mentioned earlier, in this case, the Euler scheme (55) is exact in law, so that the only source of error arises from the integration of the exponential weight. We consider the dynamics represented by the operator e ∆tW Q ∆t and e ∆t W 2 Q ∆t e ∆t W 2 with Q ∆t = e ∆tL , and first compare the results of the Galerkin discretization discussed in Section 4.2. The results reported in Figure 1 confirm our predictions: the averages of ϕ converge at first and second order for the first order and second order Galerkin schemes respectively; while the eigenvalues are the same, as expected from Proposition 6, and so both converge at second order. In this case, the numerical method based on (45) is therefore more accurate than the one based on averaging W with (34) to compute the principal eigenvalue λ, which would lead to errors of order 1 in the timestep (numerical results not shown here).   We next consider the Monte Carlo scheme presented in Section 4.1, taking M = 5 × 10 4 and an integration time T = 5 × 10 2 , with Niter = ⌊ T ∆t ⌋ for each timestep ∆t. We use half of the time for burn-in, and average in time over the second half of the simulation. Moreover, for each value of ∆t, we run 30 realizations in order to reduce the variance of the estimator and to estimate error bars on the Monte Carlo estimates (the error bars are not displayed on the pictures when they are smaller than the symbol size). The choice of the function χ ∆t depends on the scheme through (59). We compare in Figure 2 the results of the Monte Carlo algorithm with the Galerkin approximation, which serves as a reference. The agreement is very good, up to small errors arising from the finiteness of the population and of the simulation time. This result was expected since, given that the integration by the Euler scheme is exact in law in this case, the Monte Carlo method must match exactly the Galerkin approximation (provided N , Niter and M are all sufficiently large).   (57) and (58) respectively. The first and second order schemes are relative to the choice of the weight χ ∆t in (59). The dashed lines correspond to linear and quadratic fits of the Monte Carlo points for small timesteps.
Situation with a strong potential. We next show an application with a non-zero drift by setting V (q) = cos(2πq) and γ = 1. Let us recall that this dynamics is non-reversible since a constant function is not the gradient of a smooth periodic potential. The other parameters are left unchanged. Concerning the Galerkin approximation, we consider the two schemes described in Section 4.2, and characterized by the matrices defined in (60). For these schemes, the eigenvalues are the same and converge at second order (so we only consider one scheme), while the averages of ϕ converge at first and second order respectively. For the Monte Carlo algorithm described in Section 4.1, we consider the three following schemes: • Q ∆t is discretized with the Euler scheme (55), and χ ∆t (q, q ′ ) = ∆tW (q) is chosen as the left point integration; in this case, the eigenvalue and the average of ϕ converge at order one, so the scheme is referred to as first order.
• Q ∆t is discretized with the second order scheme (56), and we set χ ∆t (q, q ′ ) = ∆t (W (q)/2 + W (q ′ )/2), which corresponds to a trapezoidal rule for the time integral; in this case, both the eigenvalue and the average of ϕ converge at order two, so we refer to this scheme as second order.
We present the numerical results obtained with the various schemes we consider in Figures 3 (eigenvalues) and 4 (averages of ϕ): • Concerning the eigenvalues computed with the Monte Carlo method, we indeed observe first order convergence for the first order scheme, and second order convergence for the hybrid and second order schemes.
In particular, the results of the hybrid and the second order scheme are exactly the same. The Galerkin method also converges at second order, but with a much smaller prefactor. This is due to the fact that in this case most of the error is due to the discretization of the dynamics rather than the discretization of the time integral.
• Concerning the average of ϕ, the first order scheme converges at order one, while the hybrid and second order scheme converge at order two. We would have expected the hybrid scheme to converge at first order but, once again, this is due to the fact that most of the error is due to the discretization of the dynamics, and not to the time integral -as shown by the results of the Galerkin method, which amounts to observing the error due to the discretization of the time integral only. We indeed observe first and second order convergence for the Galerkin approximation, but we see that the error is orders of magnitude smaller than the one of the Monte Carlo approximation. This explains why the Monte Carlo hybrid and second order schemes seem to provide the same results.  Situation with a weak potential. In order to obtain a better trade-off between the error due to the discretization of the dynamics and of the time integral, we run the same simulations as in the previous situations  but with a smaller potential energy V (q) = 0.02 cos(2πq), leaving the other parameters unchanged. The results are the following: • All the eigenvalues now seem to converge at second order (see Figure 5). This is due to the fact that the error due to the discretization of the dynamics is very small, and that the discretization of the time integral, which gives the dominant error term, always leads to an effective second order convergence.
• The behaviour of the average of ϕ is more interesting (see Figure 6). The Galerkin first and second order schemes provide first and second order convergence respectively. The hybrid scheme exhibits a first order convergence, that matches the Galerkin first order scheme for small timesteps. This result can be expected since the two schemes match at order one. The first order scheme also converges at first order but with a larger prefactor, which is due to the discretization of the dynamics. On the other hand, the second order Monte Carlo scheme converges at second order, like the Galerkin second order scheme, but with a different prefactor.
Conclusion. These numerical applications show the validity of our analysis on a simple test case. However, we observe numerically that the prefactor of the leading error term depends on the choice of parameters. This has the consequence that some schemes may effectively seem to exhibit an improved order of convergence than expected, while they actually have a small prefactor at leading order, depending on the discretization at hand.

Possible extensions
The analysis and simulations we performed in this work were done for SDEs with a non-degenerate noise on a torus. We however believe that most of our results could be extended to more general settings, by considering unbounded configuration spaces and/or SDEs with degenerate noise such as inertial Langevin dynamics. This requires however  both conditions on the underlying continuous dynamics (in terms of growth conditions on b and W , for instance, and Hörmander-type conditions for dynamics with degenerate noise), and on discretization schemes (it is indeed well known by now that the existence of invariant probability measures for numerical discretizations of SDEs with non globally Lipschitz drifts on unbounded spaces requires some care, see for instance the example provided in [47]). Theoretical extensions of our results along these lines are however open problems, to the best of our knowledge.
In particular, it should first be proved that the stability result stated in Assumption 3 indeed holds true. When the state space is not compact, one typically has to work on spaces of regular functions in weighted L ∞ spaces [39,31,50]. In this functional setting, the stability property should be rephrased as the invariance of S under the action of h −1 W (L * + W − λ) −1 (hW ·). In the case W = 0, which does not involve the eigenvector hW (since h0 = 1), this is already a quite technical result to obtain (see [40,41]). Here, the presence of the eigenvector hW adds a significant difficulty.
Another technical issue concerns the convergence of the numerical scheme, as addressed in Section 2.2. Theorem 1 holds true when Assumption 1 (minorization and boundedness) is satisfied for the evolution semigroup Q W ∆t . Such estimates cannot hold for discretizations of processes on R d . An interesting question is therefore to extend Theorem 1 to the case of unbounded state spaces, for instance by using Lyapunov techniques as in [31]. This is suggested by perturbative arguments in [38,39], in the situation when then underlying dynamics admits a Lyapunov function and W is not too large with respect to this function.
Another point concerning the convergence of numerical schemes is the scaling of the factor ε with ∆t, as noted in Remark 1. It would be interesting to prove that the convergence rate of the numerical scheme is uniform in the timestep, in the sense that there would exist C, δ, ∆t * > 0 such that for an appropriate norm. Such results can be proved when W = 0, relying on the approach provided in [31], see for instance [5,43,25].  Let us first give a result which shows that it suffices to prove Proposition 2 for probability measures which admit a positive and bounded density with respect to the Lebesgue measure. This results relies on the regularizing properties of the underlying diffusion.

Lemma 9.
For any α > 0, denote by Pα(D) the subspace of probability measures which admit a smooth density with respect to the Lebesgue measure, and whose density is bounded below by α > 0 and bounded above by 1/α. Then there exists α * > 0 such that Φ1(µ) ∈ Pα * (D) for any µ ∈ P(D).
Proof. Note that, for any ϕ ∈ S, where pt is the integral kernel of the semigroup e t(L+W ) . By parabolic regularity (see for instance [24]), the integral kernel is smooth for any t > 0. It is also positive when W = 0 by the irreducibility properties of the underlying non-degenerate diffusion: there exists η > 0 such that (setting t = 1) Given that W is bounded, a similar property holds for p W t : there exists α > 0 such that Since, for any bounded measurable function ϕ, it follows that Φ1(µ) has a smooth density with respect to the Lebesgue measure, denoted by F1,µ: which gives the claimed result.
We can now provide the proof of Proposition 2.
Proof of Proposition 2. In view of the semigroup property Φt(µ) = Φt−1(Φ1(µ)) when t 1, it is sufficient by Lemma 9 to prove the result for measures µ ∈ Pα(D), where α > 0. The proof is conducted in two steps: we first prove a convergence result for the linear semigroup P W −λ t in L 2 (ν) and any times t 0, and then rely on the fact that any probability measure in Pα(D) is equivalent to ν to obtain (9).
Introduce the projector (different from the one defined in (38)) A simple computation shows that ΠW commutes with L + W and P W −λ t . It is easily seen that the spectrum of the operator ΠW (L + W − λ) ΠW is and that the associated semigroup satisfies e t Π W (L+W −λ) Π W = P W −λ t ΠW . By the Hille-Yosida theorem (see for instance [54,8]), there exists therefore a constant C > 0 such that, for any ϕ ∈ L 2 (ν), We now show that (64) implies the convergence result (9) for the class of probability measures Pα(D). For a given µ ∈ Pα(D) and ϕ ∈ L 2 (ν), where we used a Cauchy-Schwarz inequality on L 2 (ν) to go from the first to the second line. The latter computation shows that, for any µ ∈ Pα(D) and ϕ ∈ L 2 (ν), there are functions at, bt for which with |at| K ϕ L 2 (ν) e −δ W t and |bt| Ke −δ W t for some constant K > 0 independent of µ and ϕ. Moreover, there exists ε > 0 such that ε ĥ W 1/ε. Note also that |bt| ε 2 for t ln K/ε 2 /δW and that it follows that, for t ln(2K/ε)/δW , The inequality allows to obtain the desired conclusion.
Let us conclude this section with the proof of Proposition 4.
Proof of Proposition 4. The exponential convergence result (64) implies that the operator ΠW (L + W − λ) ΠW is invertible on Ran( ΠW ) = L 2 W (ν) with inverse given by The solution to (L + W − λ)u = g with g ∈ L 2 W (ν) then admits a unique solution in L 2 W (ν). By elliptic regularity, u ∈ SW when g ∈ SW . The result for L * + W − λ can be obtained by a similar reasoning.
Note that, alternatively, it would have been possible to resort to the Fredholm alternative to prove Proposition 4.

Proof of Theorem 1
Theorem 1 is a rewriting of [11,Corollary 2.5], which is stated in the context of a finite state space. In order for the paper to be self-contained, we prove Theorem 1 in our setting of continuous but compact state space, and in the simplified case of a time-homogeneous Markov chain, adapting the arguments of [11]. The idea is to prove some contraction property using the Dobrushin coefficient defined in Appendix A and the reformulation (66) below of the semigroup. We work on the space of probability measures P(D) endowed with the total variation distance.
Define the weights gn = (Q W ∆t ) n 1, and the Markov operator Sn as The dynamics (14) can then be rephrased as This equality can be proved by induction. The result is clear for n = 0. For n = 1, we have that Assuming that µn satisfies (66) at rank n, using (14) and recalling the definition (16), so that, using the recursion hypothesis and Q W ∆t gn = gn+1, it follows which concludes the recurrence. We next introduce the familly of operators Tn : P(D) → P(D) defined by: so that from (66) we have µn = Φ W ∆t,n (µ) = µTnKn. Using the definitions of Appendix A, we obtain, for two initial measures µ, ν ∈ P(D), µn − νn TV = µTnKn − νTnKn TV |||Kn||| µTn − νTn TV .
The last step consists in using Assumption 1 in order to obtain a lower bound on α(S k ) independent of k. First, for all q ∈ D and A ⊂ D, Then, it follows from definition (87) that since the infimum is taken over partitions (Ai) m i=1 of D. As a result, we obtain that, for all measures µ, ν ∈ P(D), Setting ν = µm for m ∈ N and using the semigroup property, we get so that (µn) n 1 is a Cauchy sequence in P(D). By completeness of P(D) for the total variation norm, we can conclude that, for any initial measure µ0, there exists µ∞ such that µn → µ∞ in total variation norm. Then using the one step formulation of the dynamics (17) and the semigroup property, we obtain with the choice ν = Kµ0, so that, taking n → ∞ and using the continuity of K on P(D) endowed with the total variation norm, it follows that µ∞ = Kµ∞. Passing to the limit m → +∞ in (68), Finally, it follows from (67) that the limit µ∞ does not depend on the initial distribution µ0.

Proof of Lemma 3
The idea is to approximate at leading order the stationary measure ν W,∆t as (1 + ∆t p f )νW , since we expect the invariant probability measure to be correct at order p. We start from the stationarity equation (35) and search for a function f ∈ S and a remainder R W,∆t : S → R satisfying (26) such that, for all φ ∈ S, In view of the expansion (29) of Q W ∆t and of the invariance relation (32), the first term of the left hand side is where R W,∆t gathers the terms of order at least p + 2, and is uniformly bounded in ∆t for 0 < ∆t ∆t * in the sense of (26) when f ∈ S. On the other hand, the second term on the left hand side of (69) can be written as, using again (32), where R W,∆t is uniformly bounded in ∆t in the sense of (26) when f ∈ S. We can now equate the different orders in powers of ∆t on both sides of (69) and choose f such that only a remainder of order p + 2 remains. The terms a k ∆t k D φ dνW cancel, so the first non-trivial condition to be satisfied to eliminate terms of order ∆t p reads This equality is satisfied for all φ ∈ S if and only if (take e.g. φ = f ) The condition arising from the equality of terms of order ∆t p+1 is Using that A W 1 1 = W along with condition (32), we have a1 = λ. In addition, taking adjoints in L 2 (ν) and recalling A W Moreover, in view of (70), one can subtract λ D f hW dν D φ dνW from the right hand side of last equation. Finally, we obtain the following equation (with unknown f ): for all φ ∈ S, By Assumption 3, the operator A * 1 + W − λ is invertible onŜW and leaves this space invariant. We can therefore define a solution f0 to the following equation: The function hW f0 is uniquely defined inŜW by Assumption 3 sinceg has average 0 with respect toνW , and one can check that it is indeed solution of (71). Since the eigenvector hW is regular with hW > 0, the function f0 belongs to S. However, f0 is not a priori of average 0 with respect to νW , so that condition (70) is not satisfied. We can however consider the function fα = f0 + α, which is still such that (71) holds. The choice α = − D f0 dνW ensures that (70) is satisfied. This provides the solution (33) and concludes the proof.

Proof of Lemma 4
We start by considering (23) and (37) for ϕ = ΠW φ with φ ∈ S: and We next stabilize the operator in SW by another application of the projector ΠW . First, thanks to (73). Second, since f has average 0 with respect to νW , In view of (74), the first term of the right hand side of the above equation is a remainder of order ∆t p+1 . Therefore, since ΠW φ has average 0 with respect to νW . Combining (76) with (75), where R W,∆t satisfies (26). This concludes the proof of the lemma.

Proof of Lemma 5
The first part of the proof of Lemma 5 consists in constructing an approximate eigenvectorĥ W,∆t ofĥW for the evolution operator Q W ∆t . We use to this end Assumption 4 and (32), as well as the definition of the leading order correction f in (33). More precisely, we considerĥ W,∆t = u0 + ∆t u1 + . . . + ∆t p up ∈ and look for functions u1, . . . , up ∈ SW and u0 ∈ S with D u0 dν = 1 such that Q W ∆tĥW,∆t = e ∆tλ ∆tĥ W,∆t + ∆t p+2 r W,∆t , with r W,∆t L ∞ (D) C for 0 < ∆t ∆t * . Recall that, by (36), Expanding the left hand side of (77) using (29) leads to with the convention A W 0 = Id and up+1 = 0. The right hand side of (77) can be expanded as e ∆tλ ∆tĥ W,∆t where we introduced λ0 = 1, and λp+1 is defined in (46): We see from (32) that λm = am for m ∈ {1, . . . , p}, with in particular λ1 = D W dνW = λ. We now build the functions um by induction. Let us show the first steps of the recurrence, before proceeding to the general argument. Plugging (78) and (79) in (77), the equality of terms of order 1 leads to the trivial equality u0 = u0. Equating terms of order ∆t gives In view of Assumption 4, we can conclude that u0 =ĥW . The identification of terms of order ∆t 2 in (78)-(79) leads to A W 2 u0 + A W 1 u1 + A W 0 u2 = λ2u0 + λ1u1 + λ0u2, which can be rewritten as where the expression of λ2 is given by (80) when p 2 and by (46) when p = 1. In order to prove that (81) is well-posed, it is sufficient to show that g1,0 belongs to SW . We show in fact in the sequel that each function u k is solution to a Poisson equation similar to (81) with a right-hand side that always belongs to SW . Let us now present the inductive construction to any order, until the terminal case k = p, showing in particular the well-posedness of the equations defining each mode u k . This construction is reminiscent of techniques used to build the expansion of the invariant probability measure in ∆t in related works, in particular [9]. Suppose that we have built functions u0, . . . , u k ∈ SW for some k 1. Inserting again (78) and (79) into (77) and equating terms of order ∆t k+1 then leads to For m = k + 1, we have A W 0 u k+1 = u k+1 on the left hand side and λ0u k+1 = u k+1 on the right hand side, so that the terms of order k + 1 compensate. Taking aside the terms of order m = k leads to the equation: A sufficient condition for the solution u k to exist in SW is that g k,m ∈ SW for m ∈ {0, . . . , k − 1}. For m ∈ {1, . . . , k − 1}, a sufficient condition for that is that A W k+1−m um has average 0 with respect to νW , which is clear from (32) and the fact that um ∈ SW . It therefore only remains to show that g k,0 = −A W k+1ĥW + λ k+1ĥW belongs to SW . Two cases have to be distinguished here: (a) if k < p, then k + 1 p and we can still use the invariance relation (32) applied to φ ≡ĥW , along with the fact that λ k+1 = a k+1 : Dĥ W dνW + a k+1 Dĥ W dνW = 0.
(b) in the terminal case k = p, we cannot use (32) and λp+1 has a different expression (recall (46)). Let us compute this expression explicitly. In view of (33), and, given that W f0hW =g + λhW f0 − A * 1 (hW f0) and D W dνW = λ, Finally, using the expression ofg in (33) and D hW dν = 1, From this calculation, we obtain, with (46) At this stage, it suffices to prove that the remainder term is uniformly of order ∆t p+2 for ∆t sufficiently small. We note to this end thatĥ W,∆t =ĥW + ∆t u1 + . . . + ∆t p up, where the functions u1, . . . , up are regular and hW > 0. Given that the state space D is compact, there exists ε > 0 such thatĥW ε > 0. This implies in particular that there exists ∆t ′ > 0 such that, for any 0 < ∆t ∆t ′ , it holdsĥ W,∆t ε/2 > 0. We also know that there exists ∆t * > 0 and C > 0 such that, for any 0 < ∆t ∆t * , it holds r W,∆t L ∞ (D) C. As a result, for 0 < ∆t min(∆t ′ , ∆t * ),

Proof of Lemma 6
We follow the strategy outlined in [43,44], which uses a truncated inverse series expansion. The first step is to use the expansion of the eigenvalue e ∆tλ ∆t as in the proof of  [43,44] and it suffices to write the formal series expansion of the inverse of A + ∆tB ∆t = (Id + ∆tB ∆t A −1 )A up to order p by setting S ∆t = A −1 p n=0 (−1) n B ∆t A −1 n , and then only retaining the terms of order at most ∆t p+1 in this expression. More precisely, denoting C k = ΠW (A W k − λ k )ΠW , we find where the operators C k are defined implicitly using the operators C k and A −1 . The operator S ∆t is well defined and leaves SW invariant since each C k consists in a finite number of applications of operators of the form C k A −1 and a final application of A −1 = ΠW (A1 + W − λ) −1 ΠW . It is then easy to check that, by construction, the equality (43) is satisfied.

Proof of Proposition 5
We first show that, if Q W ∆t satisfies Assumption 1 with a reference probability measure η, then Q W ∆t satisfies Assumption 1 with the same measure η. By Assumption 1, there exist ε > 0 and a measure η ∈ P(D) such that, for any bounded measurable nonnegative function ϕ, so that, applying U W ∆t on the right of Q W ∆t and T W ∆t on the left, Using (50) leads to εα 2 η(ϕ) Q W ∆t ϕ α −2 ε −1 η(ϕ), so that Q W ∆t satisfies Assumption 1. In view of Theorem 1, the scheme Q W ∆t admits a unique invariant probability measureν W,∆t and an eigenvalueλ ∆t defined by (51) . Now, integrating (84) with respect to ν W,∆t and using (21) gives εη(ϕ) e ∆tλ ∆t ν W,∆t (ϕ) ε −1 η(ϕ).
The same reasoning holds forν W,∆t . There exists therefore ε ′ > 0 for which the following inequalities hold in the sense of positive measures: We are now in position to prove the equality of the eigenvalues λ ∆t andλ ∆t defined respectively by (22) and (51). From (21) (86) It remains to note that the right hand side of (86) is uniformly bounded in n. Denoting by ψn = Q W ∆t n U W ∆t ϕ for a positive ϕ ∈ S, we obtain using (50) and (85) Therefore, the right-hand side of (86) is uniformly bounded for all n 0, which proves that λ ∆t λ ∆t by taking the limit n → +∞. A similar reasoning leads toλ ∆t λ ∆t , hence λ ∆t =λ ∆t .