# On the Mathematical foundations of Diffusion Monte Carlo

Michel Caffarel, Pierre del Moral, Luc de Montella

► To cite this version:

Michel Caffarel, Pierre del Moral, Luc de Montella. On the Mathematical foundations of Diffusion Monte Carlo. Journal of Mathematical Physics, 2025, 66 (1), pp.013301. 10.1063/5.0202800. hal-04409602v6

HAL Id: hal-04409602

https://hal.science/hal-04409602v6

Submitted on 4 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On the Mathematical foundations of Diffusion Monte Carlo

Michel Caffarel¹, Pierre Del Moral² and Luc de Montella²,³

¹University of Toulouse, Lab. de Chimie et Physique Quantiques, Toulouse, 31062, FR. E-Mail: caffarel@irsamc.ups-tlse.fr

²Centre de Recherche Inria Bordeaux Sud-Ouest, Talence, 33405, FR. E-Mail: pierre.del-moral@inria.fr, luc.de-montella@inria.fr

³Naval Group, Bouguenais, 44340, FR. E-Mail: luc.demontella@naval-group.com

December 4, 2024

## Abstract

The Diffusion Monte Carlo method with constant number of walkers, also called Stochastic Reconfiguration as well as Sequential Monte Carlo, is a widely used Monte Carlo methodology for computing the ground-state energy and wave function of quantum systems. In this study, we present the first mathematically rigorous analysis of this class of stochastic methods on non necessarily compact state spaces, including linear diffusions evolving in quadratic absorbing potentials, yielding what seems to be the first result of this type for this class of models. We present a novel and general mathematical framework with easily checked Lyapunov stability conditions that ensure the uniform-in-time convergence of Diffusion Monte Carlo estimates towards the top of the spectrum of Schrödinger operators. For transient free evolutions, we also present a divergence blow up of the estimates w.r.t. the time horizon even when the asymptotic fluctuation variances are uniformly bounded. We also illustrate the impact of these results in the context of generalized coupled quantum harmonic oscillators with non necessarily reversible nor stable diffusive particle and a quadratic energy absorbing well associated with a semi-definite positive matrix force.

## 1 Introduction

The many-body Schrödinger equation describes interacting quantum particles. Depending on the domain of application, these particles may represent electrons in solid-state physics or quantum chemistry, nucleons in nuclear physics, atoms in quantum liquid physics, or coupled modes of oscillators in molecular spectroscopy, among the main applications. Except for trivial quantum systems, it is impossible to solve this equation analytically. The diffusion Monte Carlo method (abbreviated DMC) provides a powerful stochastic approach to numerically approximate the ground state energy and wave function of Schrödinger operators.
The DMC methodology has a long and rich history, dating back to its first mention in 1949 by Ulam and Metropolis in [1]. The idea was first implemented by Donsker and Kac [2], and by Kalos [3] in the early 1960s. Over the years, the physics community has proposed numerous variants of Diffusion Monte Carlo, known by various names such as Green's function Monte Carlo, [3, 4] Fixed-Node Diffusion Monte Carlo, [5], Pure Diffusion Monte Carlo, [6, 7] Stochastic Reconfiguration Monte Carlo, [8, 9, 10, 11] and Reptation Monte Carlo, [12] to cite the main ones. Despite their apparent diversity, all these approaches are fundamentally based, in one way or another, on the stochastic simulation of a specific implementation of the Feynman-Kac formula with importance sampling.

For a more detailed discussion on the origins and the applications of these Monte Carlo techniques in physics we refer the reader to the recent review article [13] as well as to [14, 15] and references therein.

The version of interest employed here is the DMC method with a fixed number of walkers, commonly known in physics as *Stochastic Reconfiguration Monte Carlo*; see the pioneering article by Hetherington [8], followed by Sorella and co-authors [9, 10] and by the first author and his co-workers in [11].

In mathematics, the methodology may also be referred to by different names, such as genetic algorithm with selection and mutation, population Monte Carlo or sequential Monte Carlo [16, 17, 18, 19, 20]. For a more thorough discussion on these application model areas we refer to the books [21, 22] and references therein.

These sequential Monte Carlo methods do not rely on biased variational techniques. They can be seen as a sophisticated genetic-type Monte Carlo methodology to simulate interacting quantum many-body systems. Various asymptotic results have been derived, including central limit theorems and large deviation principles, see for instance [23, 24] and [25, 26], as well as the books [27, 21, 22] for an overview.

Our work concerns less studied non-asymptotic and time-uniform problems. Recall-
ing that the estimation of ground state energies relies on the limiting behavior of the
walkers' evolution in the DMC method, it is therefore crucial to obtain uniform-in-time
convergence estimates. Despite its importance, there is a notable gap in the literature
and very few results have been proven in this respect. To the best of our knowledge,
such uniform controls are mainly valid for compact state space models, see for instance
[28, 21, 22] as well as [29]. Surprisingly, the theoretical efficiency of the DMC method has
never been verified rigorously even in basic linear-Gaussian scenarios such as the simple
and well known harmonic oscillator. In this paper, we address this gap by establishing
the first uniform-in-time convergence estimates that apply to general state space models
including the coupled harmonic oscillators presented in [30].

Our approach is partly based on recent developments on the stability of positive
semigroups presented in [31, 32], see also the analysis of generalized coupled harmonic
oscillators presented in [30]. In the present article, we provide a natural Lyapunov con-
dition that ensures the exponential stability of possibly time varying positive semigroups
on non necessarily compact state spaces (cf. (6) and the local conditions (7)). In the
context of time homogeneous positive semigroups, these conditions ensure the existence
of an unique leading eigen-triple (see for instance (12)). We underline that these results
do not rely on any reversibility-type condition, nor on some spectral theorem. They
can be seen as an extended version of Perron-Frobenius and Krein-Rutman theorems for
possibly time varying positive operators.

We present a nonlinear Markov chain interpretation of the DMC methodologies. In this interpretation, the genetic type evolution of the walkers can be seen as a mean field particle simulation of a nonlinear Markov chain (see Section 2.4). In this context, we present an auxiliary Lyapunov condition that depends on the potential function and the free evolution of the walkers that ensures the time uniform performance of the DMC methodology (cf. condition (6), as well as Theorem 1 and Corollary 2). We illustrate this condition in the context of generalized coupled harmonic oscillators for a linear diffusive-type particle and a quadratic energy absorbing well associated with a semi-definite positive matrix force. In this context, we also show that the DMC methodology may diverge when the free evolution of the walkers is unstable for any fixed number of walkers, even if the asymptotic variance of the Central Limit Theorem is uniformly bounded with respect to the time parameter (see Proposition 1, as well as Section 5.4 and Proposition 5).

In the context of absorbing wells centered at the origin, this study leads us to con-
jecture that stable free evolution transitions is a necessary and sufficient condition for
the DMC method to be uniformly convergent w.r.t. the time horizon.

Additionally, we propose and, to some extent, establish the validity of an importance
sampling transformation to overcome this difficulty. This type of technique is related but
not identical to the use of guiding wave functions in physics to direct the Monte Carlo
moves to improve the efficiency of the DMC method [33, 34]. In contrast with con-
ventional guiding waves techniques our approach is based on conditional free evolutions
transitions and survival weight potential functions (cf. Section 2.4 and Section 5.2).

The rest of the article is organized as follows: In Section 2, we provide a detailed description of the general framework in which our study is set, as well as the theoretical foundations on which our proof will be based.

Section 3 is devoted to the presentation of our main results. Section 4 is mainly concerned the detailed proofs of time-uniform estimates.

Section 5 is devoted to the application of our convergence result to generalized coupled harmonic oscillators [30, 35, 36]. These models arise in various fields such as molecular spectroscopy[37], quantum optics [38], quantum cryptography [39] and photosynthesis [40]. In signal processing, the harmonic oscillator and the DMC methods coincides with the Kalman and the particle filter [41, 42].

2 Description of the models

2.1 Free evolution semigroups

Consider a Markov chain $X_n$ indexed by $n \in \mathbb{N}$ and taking values in a locally compact Polish space $(E, \mathcal{E})$, where $\mathcal{E}$ is the Borel $\sigma$-field on $E$. Let $C(E)$ be the algebra of
continuous measurable functions on E. We also define $C_b(E) \subset C(E)$ as the sub-algebra of bounded measurable continuous functions endowed with the supremum norm $||.||$. With a slight abuse of notation, we denote by 0 and 1 the null and unit scalars as well as the null and unit functions on E and we denote by $I: x \in E \mapsto I(x) = x$ the identity function on E. Throughout, we will use $a, b, c, c_1$, or $c_2$ to represent positive constants, whose values may vary from line to line.

For $n \in \mathbb{N}^*$, we consider the Markov transitions $P_n$ associated with $X_n$, and assume that they are Feller; in the sense that for any $f \in C_b(E)$ we have $P_n(f) \in C_b(E)$, with the function $P_n(f)$ defined for any $x \in E$ by the integral operator

$$P_n(f)(x) := \int_E P_n(x, dy)f(y) = \mathbb{E}(f(X_n) | X_{n-1} = x).$$

Let $C_\infty(E) \subset C(E)$ be the sub-algebra of uniformly positive continuous functions $V$ that grow at infinity; that is, for any $r \ge V_* := \inf_E V > 0$, the $r$-sub-level set $V(r) := \{V \le r\} \subset E$ is a non-empty compact subset. We further assume that there exists a $P$-Lyapunov function $V \in C_\infty(E)$; in the sense that $V(E) \subset [1, \infty)$ and there exists $\epsilon \in [0, 1)$ and $c \in \mathbb{R}$ such that for any $n \in \mathbb{N}^*$ we have

$$P_n(V) \le \epsilon V + c. \quad (1)$$

Let $C_V(E) \subset C(E)$ be the sub-space of functions $f \in C(E)$ such that $f/V$ is bounded, equipped with the norm $||f||_V := ||f/V||$. The Markov semigroup associated with the Markov chain $X_n$ is defined for any $f \in C_V(E)$ by

$$P_{k,n}(f)(x) := \mathbb{E}(f(X_n) | X_k = x).$$

Condition (1) ensures that $P_{k,n}$ is $V$-Feller in the sense that for $f \in C_V(E)$ we have $P_{k,n}(f) \in C_V(E)$. To ensure the semigroup $P_{k,n}$ is exponentially stable [31], we assume the integral operator

$$P_n(x, dy) = p_n(x, y)\nu(dy)$$

has a density $p_n$ w.r.t. some Radon measure $\nu$ satisfying for some $r_1 > 0$ and for any $r \ge r_1$ the local minorization condition

$$\left\{ \begin{array}{l} 0 < \inf_{n \in \mathbb{N}^*} \inf_{V(r)^2} p_n \le \sup_{n \in \mathbb{N}^*} \sup_{V(r)^2} p_n < \infty \\ 0 < \nu(V(r)) < \infty. \end{array} \right. \quad (2)$$

The $V$-norm semigroup contraction techniques developed in Section 8.2 in [43] (see also Lemma 2.3 in [32] and Theorem 2.2 in [31]), ensure that for any $\mu \in P_V(E)$, there exists some parameters $a, b > 0$ such that for any $k \le n$, and any $\mu_1, \mu_2 \in P_V(E)$ we have

$$\|\mu_1 P_{k,n} - \mu_2 P_{k,n}\|_V \le ae^{-b(n-k)} \|\mu_1 - \mu_2\|_V. \quad (3)$$

Note that the r.h.s condition in (2) is met as soon as $V$ has compact sub-level sets with non empty interior and $\nu$ is a Radon measure with full support; that is $\nu$ is finite
on compact sets and strictly positive on non-empty open sets. For time-homogeneous models, the l.h.s. minorization condition is satisfied as soon as $(x, y) \in (E^\circ)^2 \mapsto p_n(x, y)$ is a continuous positive function on the interior $E^\circ$ of the set $E$.

## 2.2 Feynman-Kac semigroups

We denote by $C_0(E) := \{1/V : V \in C_\infty(E)\} \subset C_b(E)$ the sub-algebra of bounded continuous positive functions $h$ that vanish at infinity; that is, for any $0 < \epsilon \le \|h\| < \infty$ the $\epsilon$-super-level set $\{h \ge \epsilon\} \subset E$ is a non empty compact subset.

Consider a family of strictly positive bounded continuous functions $G_n$ such that

$$G_n P_{n+1}(V)/V \in C_0(E) \quad (4)$$

Condition (1) ensures that this is met for any family of functions in $C_0(E)$.

We associate with it the discrete time generation Feynman-Kac semigroups

$$Q_{k,n}(f)(x) = \mathbb{E} \left( f(X_n) \prod_{p=k}^{n-1} G_p(X_p) \mid X_k = x \right)$$

and

$$\hat{Q}_{k,n}(f) := G_k^{-1} Q_{k,n}(G_n f),$$

with $G_k^{-1} := 1/G_k$. Note that condition (4) implies that for any positive function $f \in C_V(E)$ we have $Q_{n,n+1}(f)/V \le Q_{n,k}(V)/V \in C_0(E)$. Therefore, $Q_{n,n+1}$ is V-Feller. In the same vain, we check that $Q_{k,n}$ is V-Feller.

To simplify notation, for $k = (n-1)$ we sometimes write $Q_n$ and $\hat{Q}_n$ instead of $Q_{n-1,n}$ and $\hat{Q}_{n-1,n}$. In this notation, we have

$$Q_n(f)(x) = G_{n-1}(x) P_n(f)(x) \quad \text{and} \quad \hat{Q}_n(f)(x) = P_n(G_n f)(x).$$

We also use the convention $Q_{n,n} = P_{n,n} = \text{Id}$, the identity operator.

Let $M_b(E)$ be the set of bounded signed measures on $E$. Also, let $P(E) \subset M_b(E)$ be the convex subset of probability measures on $E$ and denote by $P_V(E)$ the convex set of probability measures $\mu \in P(E)$ such that $\mu(V) < +\infty$. The left action of $Q_n$ on $P_V(E)$ is given for any $(\eta, f) \in (P_V(E), C_V(E))$ by the formula

$$(\eta Q_n)(f) := \eta(Q_n(f)) = \int \eta(dx) Q_n(f)(x) \\ = \int \eta(dx) Q_n(x, dy) f(y). \quad (5)$$

By Fubini's theorem, the integration order doesn't matter. Thus to simplify notation, we sometimes write $\eta Q_n(f)$ instead of $(\eta Q_n)(f)$ or $\eta(Q_n(f))$.

We further assume the Lyapunov function $V$ introduced in (1) is a Q-Lyapunov function in the sense that (1) holds and there exists $\Theta \in C_0(E)$ and a compact subset
$K \subset E$ such that for any $n \ge 1$ we have

$$ \left\{ \begin{array}{l} Q_n(V)/V \le \Theta \\ (G_{n-1}(x) - G_{n-1}(y))(1_{E \setminus K}(x)V(x) - 1_{E \setminus K}(y)V(y)) \le 0. \end{array} \right. \qquad (6) $$

Note that the second condition in (6) holds as soon as there exists $G \in C_0(E)$ such that for any $G_n \le G$, for any $n \ge 0$. This condition ensures that for any positive function $f \in C_V(E)$ and $n \ge 1$ we have

$$ Q_n(f)/V \le Q_n(V)/V \le \Theta. $$

By (2) the integral operator $Q_n(x, dy) = q_n(x, y) \nu(dy)$ also has a density given by $q_n(x, y) = G_{n-1}(x)p_n(x, y)$ and for any $r \ge r_1$ we have the local condition

$$ 0 < \inf_{V(r)^2} q_n \le \sup_{V(r)^2} q_n < \infty. \qquad (7) $$

Consider the normalized measure valued process $\eta_n \in P_V(E)$ starting at $\eta_0 \in P_V(E)$ defined for any $n \ge 1$ by

$$ \eta_{n+1} = \phi_{n+1}(\eta_n) := \psi_{G_n}(\eta_n)P_{n+1} \quad \text{and} \quad \hat{\eta}_n := \psi_{G_n}(\eta_n), \qquad (8) $$

with the updated Boltzmann-Gibbs transformations $\psi_{G_n}$ associated with the potential function $G_n$ defined by

$$ \psi_{G_n}(\eta_n)(dx) := \frac{1}{\eta(G_n)} G_n(x) \eta(dx). $$

We readily check that the evolution semigroup $\phi_{k,n} = \phi_{k+1,n} \circ \phi_k$ associated with the flow of measure $\eta_n$ is given for any $k \le n$ by the formula

$$ \phi_{k,n}(\eta_k) = \frac{\eta_k Q_{k,n}}{\eta_k Q_{k,n}(1)}. $$

Note that for any $\mu \in P_V(E)$ we have the updating formula

$$ \phi_{k,n}(\eta_k) = \psi_{H_{k,n}^{\mu}}(\eta_k) \bar{Q}_{k,n}, \qquad (9) $$

with the Markov operator

$$ \bar{Q}_{k,n}(f) := Q_{k,n}(f)/Q_{k,n}(1) \quad \text{and} \quad H_{k,n}^{\mu}(x) := \frac{Q_{k,n}(1)(x)}{\phi_{0,k}(\mu)Q_{k,n}(1)}. \qquad (10) $$

By Theorem 4.2 in [32] (see also Theorem 1 in [44] in the context of nonlinear filtering), for any $\mu \in P_V(E)$, there exists some parameters $a, b > 0$ such that for any $k \le n$, and any $\mu_1, \mu_2 \in P_V(E)$ we have

$$ \|\mu_1 \bar{Q}_{k,n} - \mu_2 \bar{Q}_{k,n}\|_V \le a e^{-b(n-k)} \|\mu_1 - \mu_2\|_{V/H_{k,n}^{\mu}}. \qquad (11) $$
Note that for constant potential functions $G_n(x) = G_n(y)$ we have $\bar{Q}_{k,n} = P_{k,n}$ and the above contraction estimates resume to (3). Moreover, we stress that our focus will be on establishing exponential convergence rate without focusing on the estimation constant. Therefore, we are not concerned with whether certain constants in (3) and (10), such as $a$, might be large due to factors like high dimensionality. Our primary concern lies solely in their existence.

For time homogeneous models $Q_n = Q$, Theorem 4.4 in [32] ensures the existence of a leading eigen triple $(h, E_0, \eta_\infty) \in (C_V(E) \times \mathbb{R}_+^* \times P_V(E))$, such that

$$Q(h) = E_0 h, \quad \eta_\infty Q = E_0 \eta_\infty \quad \text{and} \quad \eta_\infty(h) = 1. \qquad (12)$$

## 2.3 Schrödinger semigroups

The objects defined in the previous subsection are core to a variety of physics problem. Indeed, consider an Hamiltonian differential operator $\mathcal{H}$ given by the formula

$$\mathcal{H} := -\mathcal{L} + U,$$

where $U$ is a potential energy function from $E$ to $\mathbb{R}_+$, and $\mathcal{L}$ is a kinetic energy operator acting on a subset $\mathcal{D}(\mathcal{L})$ of $\mathcal{C}(E)$. The time dependent Schrödinger equation and the imaginary time version associated with the hamiltonian $\mathcal{H}$ are given, respectively, by the equations

$$i \partial_t \psi_t(x) = \mathcal{H}(\psi_t)(x) \quad \text{and} \quad -\partial_t \varphi_t(x) = \mathcal{H}(\varphi_t)(x),$$

with prescribed initial conditions $(\psi_0, \varphi_0)$. In the above display, $i \in \mathbb{C}$ stands for the imaginary unit. The right-hand side equation is obtained via a formal time change by setting $\varphi_t(x) = \psi_{-it}(x)$, and can be equivalently written in the following form

$$\partial_t \varphi_t(x) = \mathcal{L}(\varphi_t)(x) - U(x)\varphi_t(x), \qquad (13)$$

with initial condition $\varphi_0$.

For a twice differentiable function $\varphi_0$, the solution of (13) is given by the Feynman-Kac path integral formula

$$\begin{align}
\varphi_t(x) &= Q_t(\varphi_0)(x) := \int Q_t(x, dy) \varphi_0(y) \tag{14} \\
&= \mathbb{E} \left( \varphi_0(\mathcal{X}_t) \exp \left( -\int_0^t U(\mathcal{X}_s) ds \right) \mid \mathcal{X}_0 = x \right). \nonumber
\end{align}$$

In the above display, $\mathcal{X}_t$ stands for a time homogeneous stochastic process $\mathcal{X}_t$ on $E$, with generator $\mathcal{L}$. To facilitate the interpretation of the theoretical and numerical physics in the measure theoretical framework used in this article, we note that the Feynman-Kac propagator defined by the integral operator (14) is sometimes written in terms of the exponential of the Hamiltonian operator with the exponential-type symbol

$$Q_t := e^{-t\mathcal{H}},$$
or in the bra-kets formalism

$$Q_t(\varphi_0) = |e^{-t\mathcal{H}}|\varphi_0\rangle.$$

The integral operator $Q_t$ is sometimes called the Feynman-Kac propagator. For any $s, t \ge 0$ the integral operators $Q_t$ satisfy the semigroup property

$$Q_{s+t}(x, dz) = (Q_s Q_t)(x, dy) := \int Q_s(x, dz) Q_t(z, dy) \\ \implies \varphi_{s+t} = Q_s(\varphi_t).$$

In terms of left action bra-kets, defining $\mu_{\varphi}(dx) := \varphi(x)dx$, Fubini's theorem yields

$$\begin{aligned} \langle\varphi|e^{-s\mathcal{H}}|\varphi_t\rangle &= \int dx \, \varphi(x) Q_s(x, dy) \varphi_t(dy) = (\mu_{\varphi} Q_s)(\psi_t) \\ &= \mu_{\varphi}((Q_s Q_t)(\varphi_0)) \\ &= \mu_{\varphi} Q_{s+t}(\varphi_0) = \langle\varphi|e^{-(s+t)\mathcal{H}}|\varphi_0\rangle. \end{aligned}$$

The exponential notation is compatible with finite space models and the matrix notation of the continuous one-parameter semigroup for time homogeneous models. The bra-ket notation (a.k.a. Dirac notation) is also used to represent linear projection forms acting on Hilbert spaces associated with some reversible or some stationary measure, such as the Lebesgue measure for the harmonic oscillator.

The present article deals with different types of non necessarily stationary stochastic processes, including the free evolution process $X_t$ discussed in (14). Apart in the reversible situation in which spectral theorems are stated on the Hilbert space associated with a reversible measure, the use of the exponential symbol or the use of the bra-kets formalism is clearly not adapted to represent different expectations with respect to different types of stochastic and non-necessarily reversible processes.

To analyze these general stochastic models, we have chosen to only use elementary and standard measure theory notation such as (5). The integral actions of a given integral operator $Q_t$ on the right for functions and on the left for measures are clearly compatible with finite space models and matrix notation. The left action $\mu \mapsto \mu Q_t$ maps measures into measures, while the right action $f \mapsto Q_t(f)$ maps functions into functions.

$$\left\{ \begin{aligned} (\mu Q_t)(dy) &:= \int_E \mu(dx) Q_t(x, dy) \\ Q_t(f)(x) &:= \int_E Q_t(x, dy)f(y). \end{aligned} \right.$$

In this context the normalized measures are defined for any $s \le t$ by the flow of measures

$$\mu_t(f) := \Phi_{s,t}(\mu_s)(f) := \mu_s Q_{t-s}(f) / \mu_s Q_{t-s}(1), \quad (15)$$
where $\mu_0$ stands for the distribution of the initial random state $\mathcal{X}_0$. Note that

$$ 
\begin{aligned}
\partial_t \mu_0 Q_t(1) &= -\mathbb{E} \left( U(\mathcal{X}_t) \exp \left( -\int_0^t U(\mathcal{X}_s) \, ds \right) \right) \\
&= -\mu_0 Q_t(U) \\
\implies \partial_t \log \mu_0 Q_t(1) &= -\mu_t(U) \\
\implies \mu_0 Q_t(1) &= \exp \left( -\int_0^t \mu_s(U) \, ds \right).
\end{aligned} 
\quad (16) $$

The operator $Q_t$ defined in (14) is sometimes called a Feynman-Kac propagator. However, despite its mathematical elegance, it can rarely be solved analytically. Under some regularity conditions (cf. for instance [32]) the flow of measures $\mu_t$ converge as $t \to \infty$ to some limiting fixed point measure $\mu_\infty = \Phi_{s,t}(\mu_\infty)$ (a.k.a. quasi-invariant measure). In this case, choosing $\mu_0 = \mu_\infty$ in (16) we have

$$ \mu_\infty Q_t(1) = e^{-\lambda_0 t} \quad \text{with} \quad \lambda_0 := \mu_\infty(U). $$

Whenever it exists, the ground state $h_0$ is the leading eigen-function associated with $\lambda_0$; that is, we have

$$ Q_t(h_0) = e^{-\lambda_0 t} h_0. $$

On a given time mesh $t_n = n\delta$, with time step $t_n - t_{n-1} = \delta > 0$ we clearly have

$$ \mu_{t_n}(f) = \frac{\mu_{t_{n-1}} Q_{\delta}(f)}{\mu_{t_{n-1}} Q_{\delta}(1)} \quad \text{and} \quad Q_{\delta}(x, dy) = G_{\delta}(x) P_{\delta}(x, dy), $$

with the function

$$ \mathcal{G}_{\delta}(x) := Q_{\delta}(1)(x), $$

and the Markov transition

$$ P_{\delta}(x, dy) := \frac{Q_{\delta}(x, dy)}{Q_{\delta}(1)(x)}. $$

Note that

$$ \mu_\infty(\mathcal{G}_\delta) = \mu_\infty Q_\delta(1) = e^{-\lambda_0 \delta}, $$

as well as

$$ Q_{\delta}(h_0) = e^{-\lambda_0 \delta} h_0 \quad \text{and} \quad \mu_{\infty} Q_{\delta} = e^{-\lambda_0 \delta} \mu_{\infty}. $$

In other words, choosing $(G, P) = (\mathcal{G}_\delta, \mathcal{P}_\delta)$ in the time-homogeneous Feynman-Kac model discussed in (8) and (12), we obtain the leading triple

$$ (\eta_\infty, E_0, h) = (\mu_\infty, e^{-\lambda_0 \delta}, h_0). $$

This yields the formula

$$ -\frac{1}{\delta} \log \mu_\infty(\mathcal{G}_\delta) = \lambda_0 = \mu_\infty(U). $$
Unfortunately, with the notable exception of coupled harmonic models (cf. [30] as well as Proposition 7.1 in [31]), the potential function $G_\delta$ can rarely be evaluated and the Markov transition $P_\delta$ cannot be sampled. The Feynman-Kac measure $\eta_n$ introduced in (8) can also be interpreted as the solution of a discrete-time approximation of the formula (14). Indeed, consider a discrete time approximation $X_{t_n}$ of the process $X_{t_n}$ and let

$$Q(x, dy) = G(x)P(x, dy),$$

with

$$\left\{ \begin{array}{l} G(x) = \exp(-U(x)\delta) \\ P(x, dy) = \mathbb{P}(X_{t_n}^{\delta} \in dy \mid X_{t_{n-1}}^{\delta} = x). \end{array} \right. \qquad (17)$$

In this situation, choosing $n = \lfloor t/\delta \rfloor$ we have

$$\begin{aligned} \eta_0 Q_{0,n}(f) &= \mathbb{E} \left( f(X_{t_n}^{\delta}) \prod_{0 \le k < n} G(X_{t_k}^{\delta}) \right) \\ &= \mathbb{E} \left( f(X_{t_n}^{\delta}) \exp \left\{ - \sum_{0 \le k < n} U(X_{t_k}^{\delta}) (t_{k+1} - t_k) \right\} \right) \\ &\approx_{\delta \downarrow 0} \eta_0 Q_t(f). \end{aligned}$$

The leading triple $(\eta_\infty, E_0, h)$ associated with the discrete time Feynman-Kac approximation model now depends on the time step $\delta$ and we have

$$(\eta_\infty, h) \approx_{\delta \downarrow 0} (\mu_\infty, h_0),$$

and

$$-\frac{1}{\delta} \log \eta_\infty(G) = -\frac{1}{\delta} \log \eta_\infty(e^{-U\delta}) \approx_{\delta \downarrow 0} \lambda_0 = \mu_\infty(U).$$

It is clearly out of the scope of this article to analyze the bias introduced by the discrete time approximation discussed above. We refer to [45] for a thorough discussion of this matter.

## 2.4 Diffusion Monte Carlo

The Diffusion Monte Carlo methodology relies on the fact that the flow of measures $\eta_n$ introduced in (8) can be interpreted as the probability distributions $\eta_n = \text{Law}(\bar{X}_n)$ of the random states $\bar{X}_n$ of a nonlinear Markov chain $\bar{X}_n$. The choice of the Markov chain is far from unique. For instance, we have

$$\eta_{n+1} = \phi_{n+1}(\eta_n) = \eta_n K_{n+1, \eta_n},$$
with the local Markov transition

$$
\begin{align*}
K_{n+1, \eta_n}(x, dz) &:= (S_{n, \eta_n} P_{n+1}) (x, dz) \\
&:= \int S_{n, \eta_n}(x, dy) P_{n+1}(y, dz) \\
&:= \mathbb{P}(\bar{X}_{n+1} \in dz | \bar{X}_n = x).
\end{align*}
$$

In the above display, $S_{n,\eta_n}$ stands for the Markov transition

$$
S_{n,\eta_n}(x, dy) := \epsilon_n(\eta_n)G_n(x)\delta_x(dy) \\
\qquad + (1 - \epsilon_n(\eta_n)G_n(x)) \psi_{G_n}(\eta_n)(dy).
$$

for some tuning parameter $\epsilon_n(\eta_n) \in [0, 1]$ chosen such that $\epsilon_n(\eta_n)G_n(x) \in [0, 1]$. For instance, for $]0, 1]$-valued potential functions, we can choose $\epsilon_n(\eta_n) = 0$ as well as $\epsilon_n(\eta_n) = 1$. For more general models, we can also choose the inverse of the $\eta_n$-essential supremum of $G_n$.

Note that the transition $\bar{X}_n \rightsquigarrow \bar{X}_{n+1}$ depends on the probability distributions $\eta_n$ of the random states $\bar{X}_n$. In reference with similar nonlinear Markov chain models arising in fluid mechanics, the Markov chain $\bar{X}_n$ is called a McKean interpretation of the flow of measures (8).

The mean field particle interpretation associated with a given McKean model is defined by a discrete-time system of $N$ walkers $\xi_n = (\xi_n^i)_{1\le i\le N}$. The system starts with $N$ independent copies of a random variable $\bar{X}_0 = X_0$ with distribution $\eta_0$. Given the system $\xi_n$ at some time $n \ge 0$, we sample $N$ conditionally independent walkers $\xi_{n+1}^i$ with their respective distribution

$$
K_{n+1, \eta_n^N}(\xi_n^i, dx) \quad \text{with} \quad \eta_n^N := \frac{1}{N} \sum_{1 \le i \le N} \delta_{\xi_n^i}.
$$

In other words, the DMC method consists of approximating the measure $\eta_n$ by using the occupation measure $\eta_n^N$ associated with a system of $N$ walkers. The initial positions of the walkers are randomly chosen from the distribution $\eta_0$. The evolution of each walker follows then the following selection/mutation steps:

* Selection: We evaluate the current position $\xi_n^i$ of a walker and its potential value $G_n(\xi_n^i)$. With probability $(1 - \epsilon_n(\eta_n^N)G_n(\xi_n^i))$, $\xi_n^i$ is killed and instantly replaced by another walker say $\hat{\xi}_n^i = \xi_n^j$ with a probability proportional to $G_n(\xi_n^j)$ and $j \in \{1, \dots, N\}$; otherwise we keep it and set $\hat{\xi}_n^i := \xi_n^i$.

* Mutation: We move the selected walker $\hat{\xi}_n^i = x$ to a new location $\xi_{n+1}^i = y$ using the transition kernel $P_{n+1}(x, dy)$.

The selection transition associated with the choice $\epsilon_n(\eta_n^N) = 0$ coincides with the so-called proportional selection/reconfiguration. Note that the walker with the highest
potential value is always selected when $\epsilon_n(\eta_n^N)$ is the inverse of the $\eta_n^N$-essential supremum of $G_n$.

For [0, 1]-valued potential functions $G_n$ we can also choose $\epsilon_n(\eta_n^N) = 1$. In this situation, particle are killed at a geometric clock that depends on the potential function. In the context of the discrete time approximating Feynman-Kac models discussed in (17), when the time step $\delta$ tends to 0, these geometric killing-rates converges to an exponential killing rate. The limiting DMC scheme in continuous time consists of a system of $N$ walkers. Between killing times, walkers explore the space with a free evolution with generator $L$; at rate $U$ the walkers are killed and instantly another walker in the pool duplicates (see for instance [16, 27, 46, 47, 48]). We underline that without the geometric killing rates, the variance of the proportional reconfiguration associated with the choice $\epsilon_n(\eta_n^N) = 0$ blows up when the time step tends to 0.

We expect that the occupation measures of the system approximate the solution of the measure-valued process (8); that is, in a sense to be given, for any time horizon $n \ge 0$ we have

$$ \eta_n^N = \frac{1}{N} \sum_{i=1}^{N} \delta_{\xi_n^i} \underset{N \to \infty}{\longrightarrow} \eta_n, $$

as well as

$$ \psi_{G_n}(\eta_n^N) = \sum_{i=1}^{N} \frac{G_n(\xi_n^i)}{\sum_{1 \le j \le N} G_n(\xi_n^j)} \delta_{\xi_n^i} \underset{N \to \infty}{\longrightarrow} \psi_{G_n}(\eta_n) = \hat{\eta}_n. $$

The Lyapunov drift condition (1) combined with the local minorization condition (2) ensures that the free evolution of the walkers is stable, in the sense that it forgets exponentially fast its initial condition (see (3)). This condition is not satisfied for linear gaussian processes with an unstable drift matrix. In this context, as shown in Proposition 1 and Proposition 5 (see also Section 5.4) the DMC method diverges for any fixed number of walkers, even if the asymptotic variance of the Central Limit Theorem stated in Lemma 7 is uniformly bounded w.r.t. the time horizon.

Note that the stability regions and the Lyapunov functions are connected to the potential function by (6). Importance sampling techniques and twisted guiding waves can be used to design stable-like free evolutions. For instance, using the decomposition

$$ \hat{Q}_{n+1}(f)(x) = \hat{G}_n \hat{P}_{n+1}(f), $$

with the potential function

$$ \hat{G}_n := P_{n+1}(G_{n+1}) \quad \text{and} \quad \hat{P}_{n+1}(f) := \frac{P_{n+1}(G_{n+1}f)}{P_{n+1}(G_{n+1})}, $$

we readily check the evolution equation

$$ \hat{\eta}_{n+1} = \hat{\phi}_{n+1}(\hat{\eta}_n) := \psi_{\hat{G}_n}(\hat{\eta}_n) \hat{P}_{n+1}. \quad (18) $$

This shows that the updated measures evolve as in (8) by replacing $(G_n, P_n)$ by $(\hat{G}_n, \hat{P}_n)$. The DMC associated with these objects is a genetic-type Monte Carlo sampler
with selection fitness functions $\hat{G}_n$ and mutation transitions $\hat{P}_n$. From the mathematical viewpoint, this model coincides with the one discussed in (8). Nevertheless, the free evolution of the walkers associated with $(\hat{G}_n, \hat{P}_n)$ is now driven by the potential function. This local conditioning importance sampling strategy if often used to turn infinite energy absorbing wells (a.k.a. hard obstacles) into soft ones [17, 49].

More generally, for any given time mesh $k_n \le k_{n+1}$ we have

$$ \hat{Q}_{k_n, k_{n+1}}(f) = \hat{G}_{k_n} \hat{P}_{k_n, k_{n+1}}(f), $$

with the potential function

$$ \hat{G}_{k_n} := \hat{Q}_{k_n, k_{n+1}}(1) \quad \text{and} \quad \hat{P}_{k_n, k_{n+1}}(f) := \frac{\hat{Q}_{k_n, k_{n+1}}(f)}{\hat{Q}_{k_n, k_{n+1}}(1)}. $$

This yields the formula

$$ \hat{\eta}_{k_{n+1}}(f) = \hat{\phi}_{k_n, k_{n+1}}(\hat{\eta}_{k_n}) = \psi_{\hat{G}_{k_n}}(\hat{\eta}_{k_n}) \hat{P}_{k_n, k_{n+1}}. $$

This shows that the updated measures $\hat{\eta}_{k_n}$ evolve as in (8) by replacing $(G_n, P_n)$ by $(\hat{G}_{k_n}, \hat{P}_{k_n, k_{n+1}})$. The DMC associated with these objects is a genetic-type Monte Carlo sampler with selection fitness functions $\hat{G}_{k_n}$ and mutation transitions $\hat{P}_{k_n, k_{n+1}}$.

As shown in Section 5.2 (see also Corollary 4) in the context of coupled harmonic oscillators there exists a time mesh for which the mutation transitions $\hat{P}_{k_n, k_{n+1}}$ and the potential functions $\hat{G}_{k_n}$ satisfy the required stability properties. For more general models, these objects do not have an analytic form. In this context, we can use the unbiased Monte Carlo methodologies discussed in [50], see also Section 2.3.2 in [27], and Section 11.5 in [21].

# 3 Statement of the main results

## 3.1 Some regularity conditions

For $f \in C_b(E)$, time-uniform $L_p$-convergence of the error made by the DMC method in estimating $\eta_n(f)$ have been obtained (see for example [27, 28], as well as Chapter 4 in [21], Chapter 12 in [22] and the more recent article [32]) under the strong mixing assumption that there exist $\epsilon_P \in \mathbb{R}_+$ and $\epsilon_G \in \mathbb{R}_+$ such that

$$ \left\{ \begin{array}{l} \forall (x_1, x_2, n) \in E^2 \times \mathbb{N}, \quad G_n(x_1) \ge \epsilon_G G_n(x_2) \\ P_n(x_1, dy) \ge \epsilon_P P_n(x_2, dy). \end{array} \right. $$

A significant consequence of this assumption is a time uniform bound on the potential function defined for some $\gamma \in \mathcal{P}_V(E)$ and $k \in \mathbb{N}$ by

$$ G_{k, k+n}^{\gamma} : x \in \mathbb{E} \mapsto G_{k, k+n}^{\gamma}(x) := 1/H_{k, k+n}^{\gamma}(x). $$
Unfortunately, these uniform minorization and majorization conditions are rarely satisfied when $E$ is non-compact.

In order to guarantee a time-uniform $L_p$-convergence for more general models including coupled harmonic oscillators, our framework requires to estimate uniformly in time the inverse moments of $\eta_q^N(H_{q,n}^\gamma)$. To do this, we first assume that there exists $\gamma \in \mathcal{P}_V(E)$ such that

$$ \sup_{n \in \mathbb{N}} \phi_{0,n}(\gamma) (G_n) < +\infty. \qquad (19) $$

For time-homogeneous models, without any further conditions on the potential, condition (19) is easily checked with $\gamma = \eta_\infty$. For this scenario, we will then consider in the rest that $\gamma = \eta_\infty$. Moreover, this hypothesis trivially holds if the functions $G_n$ are bounded by some constant independent of $n$.

We assume that there exists $W \in \mathcal{C}_V(E)$, $\alpha \in [0, 1]$ and a $Q$-Lyapunov function $\bar{W} \in \mathcal{C}_V(E)$ such that

$$ Q_n(W) \geq \chi \times W \quad \text{and} \quad W^{-\alpha} \leq \bar{W}, \qquad (20) $$

where

$$ \chi := \sup_{(n,x) \in (\mathbb{N}^* \times E)} P_n(G_n)(x). $$

For time homogeneous models, this condition can be relaxed into the following

$$ Q(W) \geq \min\{\chi, E_0\} \times W \quad \text{and} \quad W^{-\alpha} \leq \bar{W}. \qquad (21) $$

Note that, for time-homogeneous models, the set of functions $W \in \mathcal{P}_V(E)$ such that $Q(W) \geq E_0 W$ is non-empty as it contains at least the ground state $h$. It is also worth noting that it is not necessary to know the exact value of $\chi$ nor the one of $E_0$ in order to prove that (20) or (21) hold. Indeed, if one of these constants is less than some $C \in \mathbb{R} \cup \{+\infty\}$, then it is sufficient to prove that for any $c < C$, there exist $W_c \in \mathcal{C}_V(E)$ and $\alpha_c \in (0, 1]$ such that

$$ Q(W_c) \geq c W_c \quad \text{and} \quad W_c^{-\alpha_c} \leq V. $$

Finally, we assume that there exists a $Q$-Lyapunov function $\bar{V} \in \mathcal{C}(E)$ and $\lambda \in 2\mathbb{N}$ such that

$$ V^{\lambda} \leq \bar{V}. \qquad (22) $$

Without further mention, we assume that $V, \bar{V}$ and $\bar{W}$ are integrable with respect to $\eta_0$, i.e., $\eta_0 \in \mathcal{P}_V(E) \cap \mathcal{P}_{\bar{W}}(E) \cap \mathcal{P}_{\bar{V}}(E)$. Under conditions (19) and (20) or (21), it is then possible to obtain a time uniform bound on the random potential function $G_{k,n}^\gamma$.

**Lemma 1.** Let $\gamma$ and $(\eta, \mu)$ be defined as in (19). There exists $\bar{\beta} \in \mathbb{R}_+^*$ such that for any $\beta \leq \bar{\beta}$ we have
$$ \left\{ \begin{array}{ll} \sup_{(q,n,N) \in \mathbb{N}^3} & \mathbb{E} \left[ \phi_q(\eta_{q-1}^N)(H_{q,n}^\gamma)^{-\beta} \right] < +\infty \\ \\ \sup_{\substack{(q,n,N) \in \mathbb{N}^3 \\ q \le n}} & \mathbb{E} \left[ \eta_q^N (H_{q,n}^\gamma)^{-\beta} \right] < +\infty. \end{array} \right. \qquad (23) $$

The proof of this pivotal Lemma is postponed to the appendix.

## 3.2 A time-uniform convergence Theorem

The main goal of this paper is to establish that in the context we described, the $L_p$-norm of the error made by the DMC method in approximating the Feynman-Kac measures $\eta_n$ remains bounded in time and converges to zero as the number of particles increases. Our main result can be stated as follows

**Theorem 1.** For any $p \in \mathbb{N}^*$, there exists $c \in \mathbb{R}$ and $\beta \in (0, 1]$ such that for any $f \in C_{V^{4p}}(E)$ we have

$$ \sup_{n \in \mathbb{N}} \mathbb{E} \left( |\eta_n(f) - \eta_n^N(f)|^p \right)^{\frac{1}{p}} \le c N^{-\frac{\beta}{2}}. $$

The proof of this theorem is provided in subsection 4.1. Note that, due to the $\beta$ factor, the convergence rate obtained here may be slower than other existing convergence results cited previously. The presence of this factor arises from technical reasons in the proof of Theorem 1 and is directly related to the choice of the Lyapunov functions $V$ and $W$. It could be interesting to conduct a numerical study to determine whether there are cases exhibiting a convergence rate proportional to $N^{-\frac{\beta}{2}}$ with $\beta$ lower than one.

For time-homogeneous models, a direct consequence of Theorem 1 is a control over the estimation of the limiting quasi-invariant measure $\eta_\infty$.

**Corollary 1.** For any $p \in \mathbb{N}^*$, there exists $(c_1, c_2, \omega) \in \mathbb{R}_+^{*3}$ and $\beta \in (0, 1]$ such that for any $f \in C_{V^{4p}}(E)$ we have

$$ \mathbb{E}(|\eta_n^N(f) - \eta_\infty(f)|^p)^{\frac{1}{p}} \le c_1 N^{-\beta/2} + c_2 e^{-\omega n}. $$

Moreover, there exists $c \in \mathbb{R}$ such that, with $a = \frac{1}{\omega} \ln(c_2/c_1)$ and $b = \frac{\beta}{2\omega}$, have

$$ \sup_{n \ge a + b \ln(N)} \mathbb{E} \left( |\eta_n^N(f) - \eta_\infty(f)|^p \right)^{\frac{1}{p}} \le c N^{-\frac{\beta}{2}}. $$

The proof of the above Corollary is provided in subsection 4.2.

Assuming that there exists a measure $\mu \in \mathcal{P}(E)$ that is reversible for $P$, it becomes possible to obtain a re-normalized weak form of the ground state $h$ and its associated eigenvalue from the limit measures $\eta_\infty$ and $\hat{\eta}_\infty$ of $\eta_n$ and $\hat{\eta}_n$.
Indeed, referring to Section 9.5.5 in [43] (see also [16] as well as Section 12.4 in [21]),
we have

$$
\eta_{\infty}(G) = E_0,
$$

and

$$
\eta_{\infty}(f) = \frac{\mu(P(h)f)}{\mu(P(h))} = \frac{\mu(hP(f))}{\mu(h)} = \psi_h(\mu)P(f).
$$

Note that

$$
\eta_{\infty}(f) = \frac{\mu(Q(h)f/G)}{\mu(h)} = E_0 \psi_h(\mu)(f/G).
$$

In the reverse angle, we have the updated limiting measures

$$
\hat{\eta}_{\infty}(f) = \psi_G(\eta_{\infty})(f) = \psi_h(\mu)(f) = \frac{\mu(hf)}{\mu(h)} \implies \eta_{\infty} = \hat{\eta}_{\infty}P.
$$

The existence of a reversible measure is actually not required to express the ground-
state energy using the limit measure; indeed, we always have

$$
\eta_{\infty}(h) = \phi(\eta_{\infty})(h) = \frac{\eta_{\infty}(GP(h))}{\eta_{\infty}(G)} = E_0 \frac{\eta_{\infty}(h)}{\eta_{\infty}(G)} \Rightarrow \eta_{\infty}(G) = E_0.
$$

Those equalities, combined with the convergence stated in Corollary **1**, guarantee
the efficiency of the DMC method for approximating the ground-state energy and wave
function of quantum systems.

**Corollary 2.** Let $p \in \mathbb{N}^*$, and assume that $G \in C_{V^{\frac{\lambda}{4p}}}(E)$. There exists $(a, b, c) \in \mathbb{R}_+^{*3}$ and $\beta \in (0, 1]$ such that for any $f \in C_{V^{\frac{\lambda}{4p}}}(E)$ we have

$$
\left\{
\begin{array}{ll}
\displaystyle \sup_{n \ge a+b \ln(N)} & \mathbb{E} (|E_0 - \eta_n^N(G)|^p)^{\frac{1}{p}} \le c N^{-\frac{\beta}{2}} \\
\\
\displaystyle \sup_{n \ge a+b \ln(N)} & \mathbb{E} (|\psi_h(\mu)(P(f)) - \eta_n^N(f)|^p)^{\frac{1}{p}} \le c N^{-\frac{\beta}{2}} \\
\\
\displaystyle \sup_{n \ge a+b \ln(N)} & \mathbb{E} (|E_0 \psi_h(\mu)(f/G) - \eta_n^N(f)|^p)^{\frac{1}{p}} \le c N^{-\frac{\beta}{2}}.
\end{array}
\right.
$$

**3.3 Coupled harmonic oscillators**

To illustrate the practical applications of Theorem **1**, we carry out an in-depth study
of the generalized coupled harmonic oscillator [30]. Given $d \ge 1$, consider the flow
of measures $\mu_t$ on $\mathbb{R}^d$ defined by (15) with the the kinetic energy operator $\mathcal{L}$ and the
potential function $U$ defined for some $d \times d$ real matrices $(C, D, F)$ and $x \in \mathbb{R}^d$ by:

$$
U(x) := \exp \left( -\frac{1}{2} x' F x \right)
$$

and
$$ \mathcal{L} := \sum_{1 \le k,l \le d} C_{k,l} x_l \partial_{x_k} + \sum_{1 \le k,l \le d} D_{k,l} \partial_{x_k, x_l}, $$

with $x'$ standing for the transposition of vector $x$ (or matrix, where appropriate).

We define the normalized Markov integral operator $\bar{Q}_t$, associated with $Q_t$ from equation (14), using the ratio formula:

$$ \bar{Q}_t(f)(x) := Q_t(f)(x)/Q_t(1)(x) $$

We know from [30] that, for $t > 0$, there exists two positive-definite matrices $B$ and $S$, a matrix $A$ and $C \in \mathbb{R}_+^*$ such that:

$$ \left\{ \begin{array}{l} \bar{Q}_t(f)(x) \sim \mathcal{N}(Ax, B) \\ Qa_t(1)(x) = C \exp \left(-\frac{1}{2}x'Sx\right). \end{array} \right. $$

For any $n \in \mathbb{N}$, the sequence of measure $\gamma_{tn}$ is then equal to the flow of measures $(\eta_n)_{n \in \mathbb{N}}$ starting from $\gamma_0$ and evolving according to (8) with the the transitions kernel $P_{A,B}$ and the positive function $G_S$ defined by

$$ \left\{ \begin{array}{l} P_{A,B}(x, dy) = \frac{e^{-\frac{1}{2}(Ax-y)'B^{-1}(Ax-y)}}{(2\pi)^{k/2}|B|^{1/2}} dy \\ G_S(x) := \exp\left(-\frac{x'Sx}{2}\right) \end{array} \right. . \qquad (24) $$

For the harmonic oscillator, studying the discrete-time models instead of the continuous-time models involves then no approximation. Therefore, we will focus on discrete-time models in the following analysis.

Note that, strictly speaking, in physical terms, the harmonic oscillator framework implies that $\eta_0$ is a normal distribution and that $A$ is the null matrix. However, for the purposes of this study, within the coupled harmonic oscillator framework, we assume that the initial measure is not restricted to a normal distribution and that $A$ can be any real matrix, not necessarily the null matrix. Additionally, the matrix $A$ may not be symmetric nor a stable (Hurwitz) matrix. Furthermore, the transition $P_{A,B}$ is not necessarily reversible, unless $AB = BA'$.

In order to state a convergence result on the error of the DMC method that is not restricted to the harmonic oscillator, we now consider a family of Feller Markov transitions $(P_n)_{n \in \mathbb{N}}$ with positive densities $p_n$, an initial distribution $\eta_0 \in \mathcal{P}(E)$ and a family of positive functions $(G_n)_{n \in \mathbb{N}} \in \mathcal{C}_0(\mathbb{R}^d)^{\mathbb{N}}$ bounded by 1.

Let $(A, \mathfrak{A})$ be real matrices and $(B, \mathfrak{B}, S)$ be real positive definite matrices. Denote by $p_{A,B}$ the density of $P_{A,B}$ and by $E_{A,B,S}$ the ground-state energy associated with the operators $P_{A,B}$ and $G_S$. We assume then the following

* There exists $c_1 \in \mathbb{R}$ such that for any $(n, x, y) \in \mathbb{N} \times \mathbb{R}^d \times \mathbb{R}^d$ we have

$$ p_n(x, y) \le c_1 p_{A,B}(x, y). \qquad (25) $$
• For any $(n, x, y) \in \mathbb{N} \times \mathbb{R}^d \times \mathbb{R}^d$ we have

$$G_n(x) p_{n+1}(x, y) \geq \frac{\chi}{E_{\mathfrak{A}, \mathfrak{B}, S}} \times G_S(x) p_{\mathfrak{A}, \mathfrak{B}}(x, y). \quad (26)$$

• There exists a compact $K \subset \mathbb{R}^d$ such that for any $(n, x, y) \in \mathbb{N} \times \mathbb{R}^d \times \mathbb{R}^d$ we have

$$[G_S(x) - G_S(y)][G_n^{-1}(x)1_{\mathbb{R}^d \setminus K}(x) - G_n^{-1}(y)1_{\mathbb{R}^d \setminus K}(y)] \geq 0. \quad (27)$$

• There exists $c_3 \in \mathbb{R}_+^*$ such that $G_S^{-c_3}$ is integrable with respect to $\eta_0$.

If $P_n$ and $G_n$ are time-independent, it is possible to replace $\chi$ in (26) by the ground state energy associated with $P$ and $G$.

These conditions hold trivially for the coupled harmonic oscillator, i.e if $P_n = P_{A,B}$, $G_n = G_S$ and $\eta_0$ is a normal distribution.

In this context, Theorem 1 leads to a simple sufficient matrix condition that guarantees the uniform convergence of the DMC method. The proof of the following corollary can be found in subsection 5.1.

**Corollary 3.** Consider a family of Feller Markov transitions $(P_n)_{n \in \mathbb{N}}$ and a family of positive functions $(G_n)_{n \in \mathbb{N}} \in C_0(\mathbb{R}^d)^{\mathbb{N}}$ bounded by 1, verifying (25), (26) and (27) with $A'$SA < S$. Then, for any $p \in \mathbb{N}^*$, there exist $(\beta, \alpha, c) \in (0, 1] \times \mathbb{R}_+^{*2}$ such that for any function $f \in C_V(\mathbb{R}^d)$ we have

$$\sup_{n \in \mathbb{N}} \mathbb{E} \left( |\eta_n(f) - \eta_n^N(f)|^p \right)^{\frac{1}{p}} \leq c N^{-\frac{\beta}{2}},$$

with

$$V: x \in \mathbb{R}^d \mapsto \exp \left( \frac{\alpha}{2} x^T S x \right).$$

Returning to the coupled harmonic oscillator, we consider $P = P_{A,B}$ and $G = G_S$. Assuming the same conditions on $\eta_0$ as in Corollary 3, we aim to design an algorithm that estimates $\eta_n$ and such that its error converges to 0 for any value of $A$, $B$ and $S$. To satisfy the convergence condition of Corollary 3, we propose a modification of the DMC method by introducing a change in the transition kernel and the selection function. We define, for all $k \geq 1$, the functions $G^{(k)} \in C_V(E)$ and the Markov kernels $P^{(k)}$ on $E$, such that for all $f \in C_V(E)$, we have:

$$Q_{0,k}(f) = G^{(k)} P^{(k)}(f) \quad \text{with} \quad \begin{cases} G^{(k)} := Q_{0,k}(1) \\ P^{(k)}(f) := Q_{0,k}(f)/Q_{0,k}(1). \end{cases} \quad (28)$$

We use the convention $G^{(0)} = 1$ and $P^{(0)} = \text{Id}$ for $k = 0$.

For $k \in \mathbb{N}^*$, consider a system of walkers $\xi_n^{(k)} = (\xi_n^{(k),i})_{1 \leq i \leq N}$ associated to the DMC method with initial distribution $\psi_G(\eta_0)$, transitions $P^{(k)}$ and selection function $G^{(k)}$ as well as the empirical measures
$$ \eta_{n}^{(k),N} := \frac{1}{N} \sum_{1 \le i \le N} \delta_{\xi_{n}^{(k),i}}. \qquad (29) $$

This system of walkers offers an approximation of the measures $\eta_n$ for any $n \in k\mathbb{N}$. This type of change in the approximation, based on an importance sampling transformation is analogous to using a guiding waves function to direct the Monte Carlo moves. Without any additional condition, Theorem 1 ensures the uniform convergence of the model. The details of the proof can be found in Subsection 5.1.

**Corollary 4.** Let $p \in \mathbb{N}^*$, there exists $\bar{k} \in \mathbb{N}$ such that for any $k \ge \bar{k}$, there exist $(\beta, \alpha, c) \in (0, 1] \times \mathbb{R}_+^2$ satisfying for any $f \in C_V(\mathbb{R})$

$$ \sup_{n \in \mathbb{N}} \mathbb{E} \left( |\hat{\eta}_{nk}(f) - \eta_{n}^{(k),N}(f)|^p \right)^{\frac{1}{p}} \le cN^{-\frac{\beta}{2}}, $$

with

$$ V: x \in \mathbb{R}^d \mapsto \exp \left( \frac{\alpha}{2} x^T S x \right). $$

Although the method does not provide an approximation for every time step, several strategies can be used to fill the gaps left by the approximation. A simple approach, though more computationally intensive, is to run independent systems of walkers for each time step in the interval $[0, k-1]$. This method not only fills in the gaps, but also maintains the convergence property.

So far, our focus has been on approximating the flow of measures $\eta_n$. However, as shown in (18), the flow of measures $\hat{\eta}_n$ is defined similarly to $\eta_n$ but with $(G, P)$ replaced by $(\hat{G}, \hat{P})$. In the context of the coupled harmonic oscillator discussed earlier, denoting by $\mathcal{N}(A, B)$ the Gaussian distribution centered at $A$ with co-variance matrix $B$, we have:

$$ \left\{ \begin{array}{l} \delta_x \hat{P} = \mathcal{N}(\hat{A}x, \hat{B}) \\ \hat{G}(x) \propto \exp \left(-\frac{1}{2} x' \hat{S} x \right) \end{array} \right. \qquad (30) $$

with

$$ \left\{ \begin{array}{l} \hat{B} := (B^{-1} + S)^{-1} \le B \\ \hat{A} := \hat{B} B^{-1} A \\ \hat{S} := A' (S - S(B^{-1} + S)^{-1} S) A \ge 0. \end{array} \right. \qquad (31) $$

This shows that, within this framework, the convergence condition stated in Corollary 3 can be extended to the updated measures $(\hat{\eta}_n)_{n \in \mathbb{N}}$ by replacing $A$ and $S$ with $\hat{A}$ and $\hat{S}$. Similarly, Corollary 4 can be extended in the same way.

Note that all corollaries in this subsection can be extended to a control on the estimation of the limit measures, ground state, and eigenvalue using the same approach as presented in Corollary 2.

Our study concludes with a focus on the divergence of the DMC method when approximating the one-dimensional harmonic oscillator. This confirms that the stability
condition stated in Corollary 3 is necessary and that, in some cases, the set of assumptions presented can closely approximate a sufficient and necessary condition. Additionally, it emphasizes the significance of the importance sampling method introduced in the previous corollary. Specifically, in the one-dimensional context, the sufficient condition for uniform convergence of the DMC method is expressed by $A^2 < 1$. Proposition 1 establishes the divergence of the error made by the DMC method when $A^2 > 1$, leaving open only the case $A = 1$.

**Proposition 1.** Assume that $A^2 > 1$ and $P_0 > 0$. For any $p \in \mathbb{N}^*$ we have

$$ \sup_{n \in \mathbb{N}} \mathbb{E} \left( |\eta_n(I) - \eta_n^N(I)|^p \right)^{\frac{1}{p}} = +\infty. $$

We complement the divergence result with an observation on the asymptotic convergence. Our next proposition proves that, regardless of the value of $A$, the univariate harmonic oscillator can always be uniformly controlled over time as the number of walkers approaches infinity. This underscores the importance of achieving uniform control for a finite number of walkers, rather than focusing solely on the asymptotic scenario.

**Proposition 2.** Consider the univariate harmonic oscillator. For any $n \ge 1$, we have the following convergence in law as $N$ tends toward $+\infty$

$$ \sqrt{N}[\eta_n(I) - \eta_n^N(I)] \xrightarrow[N \to \infty]{\mathcal{L}} \mathcal{N}(0, \sigma_n^2), $$

with $\sigma_n^2$ the asymptotic variance such that

$$ \sup_{n \in \mathbb{N}} \sigma_n^2 < \infty. $$

The proof of those propositions can be found in Subsection 5.4.

# 4 Stochastic interpolation

## 4.1 Time varying semigroups

In this subsection, we focus on proving Theorem 1. To take advantage of the conditional independence of the walkers, we structure our approach around the following decomposition of the difference between the Feynman-Kac measure and its empirical approximation, using the convention $\eta_{-1}^N = \eta_0$. Following [27, 28], we use the following stochastic interpolation formula

$$ \eta_n^N - \eta_n = \sum_{q=0}^{n} [\phi_{q,n}(\eta_q^N) - \phi_{q,n}(\phi_q(\eta_{q-1}^N))]. \quad (32) $$

Each term on the right-hand side represents the error that occurs when using the DMC approximation instead of the real propagator for a single extra time step. Combining the uniform bound given in Lemma 1 with the contraction property (11), the following Lemma establishes an exponentially decreasing control for these local errors.
**Lemma 2.** For any $p \in \mathbb{N}^*$, there exists $(c, \rho, \beta) \in \mathbb{R}_+^{*2} \times (0, 1]$ such that for any function $f \in C_{V^{\frac{\lambda}{4p}}}(E)$ and any $(N, q, n) \in \mathbb{N}^3$ with $q \le n$ we have

$$ \mathbb{E} \left[ \left| \left[ \phi_{q,n}(\eta_q^N) - \phi_{q,n}(\phi_q(\eta_{q-1}^N)) \right] (f) \right|^p \right]^{\frac{1}{p}} \le c e^{-(n-q)\rho} N^{-\frac{\beta}{2}}. \quad (33) $$

**Proof:**

Let $(\eta, \mu) \in \mathcal{P}(E)$ and let $\gamma \in \mathcal{P}_V(E)$ be defined as in (19). Consider $H_{q,n} := H_{q,n}^{\gamma}$ as defined in (10). Applying the updating formula (9), we obtain

$$ \begin{aligned} & \phi_{q,n}(\eta)(f) - \phi_{q,n}(\mu)(f) \\ &= (\psi_{H_{q,n}}(\eta)\bar{Q}_{q,n} - \psi_{H_{q,n}}(\mu)\bar{Q}_{q,n})(f) \\ &= \frac{1}{\eta(H_{q,n})}(\eta - \mu)(H_{q,n}\bar{Q}_{q,n}[f - \psi_{H_{q,n}}(\mu)\bar{Q}_{q,n}(f)]). \end{aligned} $$

This yields the formula

$$ \phi_{q,n}(\eta)(f) - \phi_{q,n}(\mu)(f) = \frac{1}{\eta(H_{q,n})}(\eta - \mu)(F_{q,n}^{\mu}), \quad (34) $$

with the function

$$ F_{q,n}^{\mu}(x) := H_{q,n}(x) \int_{E} \psi_{H_{q,n}}(\mu)(dy)[\bar{Q}_{q,n}(f)(x) - \bar{Q}_{q,n}(f)(y)]. $$

Then, applying Hölder's inequality for any $\beta \in [0, 1)$, we obtain the estimate

$$ \begin{aligned} & \mathbb{E} \left( |\phi_{q,n}(\eta_q^N)(f) - \phi_{q,n}(\phi_q(\eta_{q-1}^N))(f)|^p \right)^{\frac{1}{p}} \\ & \le \mathbb{E} \left( \eta_q^N (H_{q,n})^{-2p} |(\eta_q^N - \phi_q(\eta_{q-1}^N))(F_{q,n}^N)|^{2p(1-\beta)} \right)^{\frac{1}{2p}} \\ & \times \mathbb{E} \left( |(\eta_q^N - \phi_q(\eta_{q-1}^N))(F_{q,n}^N)|^{2\beta p} \right)^{\frac{1}{2p}}, \end{aligned} $$

with $F_{q,n}^N := F_{q,n}^{\phi_q(\eta_{q-1}^N)}$. Using (34) this yields the estimate

$$ \begin{aligned} & \mathbb{E} \left( |\phi_{q,n}(\eta_q^N)(f) - \phi_{q,n}(\phi_q(\eta_{q-1}^N))(f)|^p \right)^{\frac{1}{p}} \\ & \le \mathbb{E} \left( |(\eta_q^N - \phi_q(\eta_{q-1}^N))(F_{q,n}^N)|^{2\beta p} \right)^{\frac{1}{2p}} \\ & \times \mathbb{E} \left( \eta_q^N (H_{q,n})^{-2\beta p} |\phi_{q,n}(\eta_q^N)(f) - \phi_{q,n}(\phi_q(\eta_{q-1}^N))(f)|^{2p(1-\beta)} \right)^{\frac{1}{2p}}. \end{aligned} $$
Recalling that $f \in C_{V^{4p}}(E)$, this implies that

$$
\begin{align*}
\mathbb{E} (|\phi_{q,n}(\eta_q^N)(f) - \phi_{q,n}(\phi_q(\eta_{q-1}^N))(f)|^p)^{\frac{1}{p}} & \\
& \leq \left[ \mathbb{E} (\phi_{q,n}(\eta_q^N)(V^\lambda))^{\frac{1}{4p}} + \mathbb{E} (\phi_{q,n}(\phi_q(\eta_{q-1}^N))(V^\lambda))^{\frac{1}{4p}} \right] \\
& \times \mathbb{E} (\eta_q^N (H_{q,n})^{-4\beta p})^{\frac{1}{4p}} \mathbb{E} \left( |(\eta_q^N - \phi_q(\eta_{q-1}^N))(F_{q,n}^N)|^{2\beta p} \right)^{\frac{1}{2p}}.
\end{align*}
$$

From Lemmas 1 and 9, we deduce that, to conclude, it is enough to prove that, for some constant $c \in \mathbb{R}_+^*$ independent of $n$, $q$ and $N$, we have

$$
\mathbb{E} \left[ \left[ |\eta_q^N - \phi_q(\eta_{q-1}^N)| (F_{q,n}^N)^{|2\beta p|} \right]^{\frac{1}{2p}} \right] < c e^{-c_2(n-q)} N^{-\frac{\beta}{2}}. \quad (35)
$$

Let $\beta' = \frac{2p\beta}{\lambda}$ and assume $\beta$ small enough so that $\beta' < 1/4$.

For $q > 0$, the walkers $(\xi_q^i)_{1\le i\le N}$ are independent conditionally to $\eta_{q-1}^N$ we have

$$
[\eta_q^N - \phi_q(\eta_{q-1}^N)](f) = \eta_q^N \left( \frac{1}{N} \sum_{1 \le i \le N} h_i \right),
$$

with

$$
h_i := f - S_{q-1, \eta_{q-1}^N} P_q(f)(\xi_{q-1}^i).
$$

Moreover, we have

$$
\phi_q(\eta_{q-1}^N) = \frac{1}{N} \sum_{1 \le i \le n} \mu_i \quad \text{with} \quad \mu_i = \delta_{\xi_{q-1}^i} S_{q-1, \eta_{q-1}^N} P_q.
$$

Since for any $i \in [1, N]$, we have $\mu_i(h_i) = 0$, we can apply Lemma 7.3.3 from [21], and deduce that there exists $C \in \mathbb{R}$ such that

$$
\begin{align*}
\mathbb{E} \left[ \left[ |\eta_q^N - \phi_q(\eta_{q-1}^N)| (F_{q,n}^N)^{|2\beta p|} \right]^{\frac{1}{2p}} \right] \\
&= \mathbb{E} \left[ \mathbb{E} \left( |[\eta_q^N - \phi_q(\eta_{q-1}^N)] (F_{q,n}^N)^{|\lambda\beta'|} | \xi_{q-1} \right) \right]^{\frac{1}{2p}} \\
&\le \mathbb{E} \left[ \mathbb{E} \left( |[\eta_q^N - \phi_q(\eta_{q-1}^N)] (F_{q,n}^N)^{|\lambda|} | \xi_{q-1} \right)^{\beta'} \right]^{\frac{1}{2p}} \\
&\le \frac{C}{N^{\beta/2}} \mathbb{E} \left[ \phi_q(\eta_{q-1}^N) (|F_{q,n}^N|^\lambda)^{\beta'} \right]^{1/2p}.
\end{align*}
$$

For $q=0$, the walkers are iid with common distribution $\eta_0$. The previous reasoning
therefore holds using the convention $E(X | \xi_{-1}) = \mathbb{E}(X)$.
Applying the contraction property (11) with $\mu = \delta_x$ and $\eta = \delta_y$ we get the existence of $(a, \rho) \in \mathbb{R}_+^2$ such that

$$ \begin{aligned} & |\bar{Q}_{q,n}(f)(x) - \bar{Q}_{q,n}(f)(y)| \\ & \leq ae^{-\rho(n-q)} \left(1 + \frac{V(x)}{H_{q,n}(x)}\right) \left(1 + \frac{V(y)}{H_{q,n}(y)}\right). \end{aligned} \quad (36) $$

By substituting (36) into the definition of $F_{q,n}$ and applying Hölder's inequality along with Jensen's inequality, we obtain, for some $(a', \rho') \in \mathbb{R}_+^{*2}$

$$ \begin{aligned} & a' e^{\rho'(n-q)} \mathbb{E} \left[ \phi_q(\eta_{q-1}^N) (|F_{q,n}^N|^\lambda)^{\beta'} \right]^{1/2p} \\ & \leq \mathbb{E} \left[ \phi_q(\eta_{q-1}^N) \left\{ (H_{q,n} + V)^\lambda \right\}^{\beta'} \phi_q(\eta_{q-1}^N) \left\{ (H_{q,n} + V)^\lambda \right\}^{\lambda\beta'} \right. \\ & \qquad \left. \times \phi_q(\eta_{q-1}^N) (H_{q,n})^{-\lambda\beta'} \right]^{\frac{1}{2p}} \\ & \leq \mathbb{E} \left[ \phi_q(\eta_{q-1}^N) \left\{ (H_{q,n} + V)^\lambda \right\} \right]^{\frac{1}{4p}} \mathbb{E} \left[ \phi_q(\eta_{q-1}^N) (H_{q,n})^{-2\lambda\beta'} \right]^{\frac{1}{4p}}. \end{aligned} $$

From our hypothesis on $Q_n$, we deduce from Lemma 3.2 in [32] that there exists a constant $c$ such that for any $(q', n') \in \mathbb{N}^2$, $H_{q',n'} \leq cV$. We can then conclude by choosing a small enough $\beta'$ and using Lemmas 1 and 9. $\square$

The proof of Theorem 1 is now relatively straightforward.

### Proof of Theorem 1:

Let $f \in C_{V^{\frac{\lambda}{4p}}}(E)$. From the sub-additivity of the $L_p$-norm applied in (32), we have

$$ E(|\eta_n^N(f) - \eta_n(f)|^p)^{\frac{1}{p}} \leq \sum_{q=0}^{n} \mathbb{E} \left( |\phi_{q,n}(\eta_q^N)(f) - \phi_{q,n}(\phi_q(\eta_{q-1}^N))(f)|^p \right)^{\frac{1}{p}}. $$

By applying Lemma 2, we deduce that there exists $(c, \rho, \beta) \in \mathbb{R} \times \mathbb{R}_+^* \times (0, 1]$ such that for any $N \in \mathbb{N}^*$ we have

$$ \mathbb{E}(|\eta_n^N(f) - \eta_n(f)|^p)^{\frac{1}{p}} \leq \frac{c}{N^{\frac{\beta}{2}}} \sum_{0 \leq l \leq n} e^{-(n-l)\rho} \leq \frac{c}{N^{\frac{\beta}{2}}(1 - e^{-\rho})}. $$

This ends the proof of the theorem. $\square$
## 4.2 Ground state estimates

This subsection concentrates on proving Corollary 1. We consider thus the time-homogeneous model.

Let $f \in C_{V^{\frac{\lambda}{4p}}} (E)$. Notice that we can decompose the error made by the DMC method in the following way:

$$ \mathbb{E}(|\eta_n^N(f) - \eta_\infty(f)|^p)^{\frac{1}{p}} \le \mathbb{E}(|\eta_n^N(f) - \eta_n(f)|^p)^{\frac{1}{p}} + |\eta_n(f) - \eta_\infty(f)|. $$

Theorem 1 implies that there exists $(c_1, \beta) \in \mathbb{R}^* \times (0, 1]$ such that

$$ \sup_{n \in \mathbb{N}} \mathbb{E}(|\eta_n^N(f) - \eta_n(f)|^p)^{\frac{1}{p}} \le \frac{c_1}{N^{\beta/2}}. $$

According to Theorem 4.3 in [32], there exists $(C_2, \omega) \in \mathbb{R}_+^{*2}$ such that

$$ |\eta_n(f) - \eta_\infty(f)| \le c_2 e^{-\omega n}. $$

Hence

$$ \mathbb{E}(|\eta_n^N(f) - \eta_\infty(f)|^p)^{\frac{1}{p}} \le \frac{c_1}{N^{\beta/2}} + c_2 e^{-\omega n}. $$

Thus, letting $a = \frac{1}{\omega} \ln(c_2/c_1)$ and $b = \frac{\beta}{2\omega}$, we deduce that there exists $c \in \mathbb{R}$ such that

$$ \sup_{n \ge a+b \ln(N)} \mathbb{E}(|\eta_n^N(f) - \eta_\infty(f)|^p)^{\frac{1}{p}} \le \frac{c}{N^{\beta/2}}. $$

This concludes the proof. □

# 5 Coupled harmonic oscillators

This section is devoted to a detailed study of the coupled harmonic oscillator. We recall that in this framework, we consider the flow of measures (8) starting from some initial measure $\eta_0$ on $\mathbb{R}^d$, and that the transition kernel and the potential function considered are defined for some positive definite matrices $B$ and $S$ and some matrix $A$ by (24).

This study will be carried out in four steps. First, we consider the general framework described at the beginning of Subsection 3.3. There, by performing a complete study of the Lyapunov functions associated with $P_{A,B}$ and $G_S$, we establish that under the simple condition $A'SA < S$, the DMC method can be controlled uniformly in time.

Note that in the coupled harmonic oscillator framework, if the initial measure is normal, then the entire flow of measures consists of normal distributions. More specifically, the measures $\eta_n = N(m_n, \Omega_n)$ are defined with
$$m_n = \mathcal{E}(\Omega_{n-1}) m_{n-1} \quad \text{and} \quad \Omega_n = \Phi(\Omega_{n-1}),$$

with the mappings

$$\mathcal{E}(\Omega) := A(I + \Omega S)^{-1} \quad \text{and} \quad \Phi(\Omega) := A(I + \Omega S)^{-1} \Omega A' + B.$$

Building on this, we will show that the condition $A'SA < S$ is a stability condition ensuring that $m_n$ is a decreasing sequence converging to the null vector.

Returning to the standard harmonic oscillator, we then study a variation of the DMC method. This allows for uniform control of the method without any imposed condition on the matrices $A$ and $S$. This part also leads us to fully describe the ground state $h$ as well as the ground state energy associated with $P_{A,B}$ and $G_S$.

We continue with a result on the stability of the flow of measure $\eta_n$ with respect to the initial measure for the coupled harmonic oscillator. We prove that the harmonic oscillator tends to forget its initial condition exponentially fast. This result is of both theoretical and practical interest. It justifies that the constants in Corollary 3 should not be too dependent on the initial measure, and it proves that a bad initialisation of the walkers when running a simulation should have negligible consequences.

To complete our study, we examine the univariate harmonic oscillator. In this context, we prove that the expected error of the DMC method diverges to $\infty$ as soon as $A > 1$. Additionally, we establish the asymptotic convergence of this error for any value of $A$. Together, these results not only confirm the validity of our convergence theorem but also validate our choice to study the uniform-in-time convergence of the error rather than restricting our analysis to its asymptotic convergence.

## 5.1 Lyapunov functions

This subsection is dedicated to the proof of Corollary 3. Therefore, we place ourselves within the framework associated with this corollary. We will conduct the proof assuming (20). If (21) holds with $E_0 < \chi$ then the demonstration is completely analogous.

Consider a family of Feller Markov transitions $(P_n)_{n \in \mathbb{N}}$ and a family of positive functions $(G_n)_{n \in \mathbb{N}} \in C_0(\mathbb{R}^d)^N$ bounded by 1, verifying (25), (26) and (27). It is clear from Subsection 2.2 that proving the existence of a continuous $P$-Lyapunov function $V \in C_\infty(\mathbb{R}^d)$ makes (25) and (26) sufficient conditions for (2) to hold. To guarantee the existence of an appropriate $Q$-Lyapunov function, we need a result obtained by Kato in [51]. We present it here using the formulation provided in [52].

**Lemma 3.** Suppose that $D \subset \mathbb{R}$ is an interval, and let $A$ be a continuous function from $D$ to the space of real $d \times d$ matrices. In this case, there exist $d$ eigenvalues (counted with algebraic multiplicities) of $A(t)$ which can be parameterized as continuous functions $\lambda_1(t), \dots, \lambda_d(t)$ from $D$ to $\mathbb{R}$.

We can now ensure the existence of a $Q$-Lyapunov function under a simple matrix condition.
**Lemma 4.** Assume that $A^TSA < S$. There exists $\alpha \in \mathbb{R}_+^*$ such that the function

$$V: x \in \mathbb{R}^d \mapsto \exp\left(\frac{\alpha}{2}x^T Sx\right), \qquad (37)$$

is a Q-Lyapunov function and it is integrable w.r.t $\eta_0$.

**Proof:**

From (27), we deduce that the r.h.s of (6) holds for $V$ and $G_n$. Then, together with (25), we deduce that it is enough to prove that $V$ is a Lyapunov function for $P_{A,B}$ with $\epsilon < 1/c$. Let's then compute $P_{A,B}(f)$ for any function of the form

$$f: x \in \mathbb{R}^d \mapsto \exp\left(\frac{1}{2}x^T Fx\right),$$

where $F$ is an invertible matrix such that $B^{-1} - F$ is positive definite.

In this setting, the Woodbury matrix identity provides the following equality:

$$(B - F^{-1})^{-1} = B^{-1} - B^{-1} (B^{-1} - F)^{-1} B^{-1}.$$

We have then for any $x \in \mathbb{R}^d$, with $\bar{B}_F := (B^{-1} - F)^{-1}$:

$$\begin{aligned} P_{A,B}(f)(x) &= \frac{1}{(2\pi)^{d/2} \det(B)^{1/2}} \int_{\mathbb{R}^d} e^{-\frac{1}{2}[(Ax-y)^\top B^{-1}(Ax-y) - y^T F y]} dy \\ &= f(x) \frac{\exp(\frac{1}{2}x^T(A^T B^{-1} \bar{B}_F B^{-1} A - A^T B^{-1} A - F)x)}{(2\pi)^{d/2} \det(B)^{1/2}} \\ & \quad \times \int_{\mathbb{R}^d} \exp\left(-\frac{1}{2}\left[(\bar{B}_F B^{-1} A x - y)^\top \bar{B}_F^{-1} (\bar{B}_F B^{-1} A x - y)\right]\right) dy \end{aligned}$$

This yields the formulae

$$\begin{aligned} P_{A,B}(f)(x) &= \sqrt{\frac{\det(\bar{B}_F)}{\det(B)}} f(x) \\ & \quad \times \exp\left(\frac{1}{2}x^T(A^T[B^{-1}(B^{-1} - F)^{-1}B^{-1} - B^{-1}]A - F)x\right) \\ &= \frac{1}{\sqrt{\det(I_d - BF)}} f(x) \exp\left(-\frac{1}{2}x^T(F + A^T(FB - I_d)^{-1}FA)x\right). \end{aligned} \qquad (38)$$
From those calculations, we deduce that $V$ is a Lyapunov function for $P_{A,B}$ if the matrices $B^{-1} - \alpha S$ and $S - A^T(I_d - \alpha SB)^{-1}SA$ are positive definite.

Let $\lambda_B$ be the greatest eigenvalue of $B$ and $\lambda_S$ be the greatest eigenvalue of $S$. It is clear that for $\alpha \in (0, \frac{1}{\lambda_B\lambda_S})$, $B^{-1} - \alpha S$ is positive definite.

Consider now the function

$$ \psi: \alpha \in \left[0, \frac{1}{\lambda_B \lambda_S}\right) \mapsto \text{sp}(S - A^T(I_d - \alpha SB)^{-1}SA) \in \mathbb{R}^d. $$

Here, $\text{sp}(M)$ represents the spectrum of a matrix $M$ with multiplicity taken into account.

Given the hypotheses on $A$ and $S$, we can conclude that $\psi(0) \subset \mathbb{R}_+^{*d}$. Furthermore, by Lemma 2, it is clear that $\psi$ is a continuous function. Since $\mathbb{R}_+^{*d}$ is an open set, there exists $\bar{\alpha} \in \mathbb{R}$ such that for any $\alpha \in (0, \bar{\alpha})$, $\psi(\alpha) \subset \mathbb{R}_+^{*d}$.

By choosing a sufficiently small value for $\alpha$ to ensure that $V$ is integrable w.r.t $\eta_0$, we can conclude. $\square$

From this Lemma and the hypothesis on $G_n$, we deduce that the l.h.s of (6) holds as well. We can now focus on verifying that (21) holds by proving the following Lemma.

**Lemma 5.** There exists a positive definite matrix $H$ such that

$$ \forall n \in \mathbb{N}, \quad Q_n(W) \geq \chi \times W \quad \text{with} \quad W(x) := \exp(-\frac{1}{2}x^T H x) $$

**Proof:**

From (26), we have

$$ Q_n(W) \geq \frac{\chi}{E_{A',B',S}} \times G_S P_{A',B'}(W). $$

Using (38), we derive the following expression for $G_S(x)P_{A',B'}(W)(x)$ with $x \in \mathbb{R}^d$:

$$ \frac{W(x) \exp \left( \frac{1}{2} x^T (H - A'^T (HB' + I_d)^{-1} H A' - S) x \right)}{\sqrt{\det(I_d + B'H)}} $$

Hence, choosing $H$ as the solution to the Riccati equation

$$ H - A'^T(HB' + I_d)^{-1}HA' - S = 0, $$

we deduce that

$$ E_{A',B',S} = \frac{1}{\sqrt{\det(I_d + B'H)}}. $$

Thus:

$$ Q_n(W) \geq \chi \times W. $$
For $\alpha$ sufficiently small, $W^{-\alpha}$ is lower than $V$. The right-hand side of (21) is then also verified.

Under the condition $A^T S A < S$, we have confirmed that all the assumptions concerning $P_n$, $G_n$, and $\eta_0$ in Theorem 1 are satisfied. We can therefore apply it to conclude on the proof of Corollary 3.

Although justified by previous calculations, the meaning and implication of this matrix condition may not be immediately apparent. To remedy this, we consider the framework of the coupled harmonic oscillator to provide a property that illustrates how this condition is a stability condition that ensures that the mean of the flow of measures $\eta_n = \mathcal{N}(m_n, \Omega_n)$ is a strictly decreasing sequence converging to zero.

**Proposition 3.** Let $B$ and $S$ be positive definite matrices and $A$ be a matrix. Consider a Gaussian distribution $\eta_0 = \mathcal{N}(m_0, \Omega_0)$ and consider the flow of measures $\eta_n = \mathcal{N}(m_n, \Omega_n)$ given in (8) with $P = P_{A,B}$ and $G = G_S$. The following implication holds

$$A'SA < S \implies \exists C \in [0,1) \text{ s.t. } \forall k < n, \|m_n\| \le C^{n-k} \|m_k\|.$$

**Proof :**

The flow of measures $\eta_n = \mathcal{N}(m_n, \Omega_n)$ are defined with

$$m_n = \mathcal{E}(\Omega_{n-1}) m_{n-1} \quad \text{and} \quad \Omega_n = \Phi(\Omega_{n-1}),$$

with the mappings

$$\mathcal{E}(\Omega) := A(I + \Omega S)^{-1} \quad \text{and} \quad \Phi(\Omega) := A(I + \Omega S)^{-1} \Omega A' + B.$$

We set

$$\mathcal{E}_{k,n}(\Omega) := \mathcal{E}(\Phi^{n-1}(\Omega)) \dots \mathcal{E}(\Phi^{n-k}(\Omega))$$

as well as

$$\mathcal{E}_n(\Omega) = \mathcal{E}_{0,n}(\Omega)$$

and

$$\Phi^n = \Phi \circ \Phi^{n-1}.$$

In this notation, we have

$$m_n = \mathcal{E}_{k,n}(\Omega_0) m_k.$$

Note that

$$A'SA < S \iff \overline{A'A} < I \quad \text{with} \quad \overline{A} := S^{1/2}A S^{-1/2}.$$

In the same vein, we have

$$\mathcal{E}(\Omega) = S^{-1/2} \bar{\mathcal{E}}(\Omega) S^{1/2}$$

with

$$\bar{\mathcal{E}}(\Omega) := \bar{A} (I + S^{1/2} \Omega S^{1/2})^{-1}.$$
Let $\|M\|_2 := \sqrt{\lambda_{max}(M'M)}$ be the spectral norm, where $\lambda_{max}(S)$ stands for the maximal eigenvalue of a symmetric matrix $S$. In this notation, we have the estimate

$$ \begin{aligned} \|\bar{\mathcal{E}}(\Omega)\|_2 &\le \|\bar{A}\|_2 \lambda_{max}((I + S^{1/2} \Omega S^{1/2})^{-1}) \\ &\le \frac{\|\bar{A}\|_2}{\lambda_{min}(I + S^{1/2} \Omega S^{1/2})} \le \|\bar{A}\|_2 < 1. \end{aligned} $$

Using the easily checked decomposition

$$ \mathcal{E}_{k,n}(\Omega) = S^{-1/2} (\bar{\mathcal{E}}(\Phi^{n-1}(\Omega)) \dots \bar{\mathcal{E}}(\Phi^{n-k}(\Omega))) S^{1/2}, $$

this yields the rather crude exponential decay

$$ \|\mathcal{E}_{k,n}(\Omega)\|_2 \le \|\bar{A}\|_2^{n-k}. $$

We conclude that

$$ A'SA < S \implies \|m_n\| \le \|\bar{A}\|_2^{n-k} \quad \|m_k\| \to_{n \to \infty} 0. $$

## 5.2 Conditional free evolutions

This subsection focuses on the study of the importance sampling described in Section 2.4, Section 3.3 and on the proof of Corollary 4.

We consider the coupled harmonic oscillator, i.e, for some matrices $A$, $B$ and $S$, with $B$ and $S$ symmetric definite positive, we consider $P = P_{A,B}$ and $G = G_S$. Up to this point, we have established that the $L_p$-norm of the error made by the DMC method is uniformly bounded in time, with a convergence rate of $\frac{1}{N^{\beta/2}}$ for some $\beta \in (0, 1]$ when $A^T S A < A$.

Our aim is now to use Theorem 1 to prove that the approximation of the measures $\hat{\eta}_n$ made by the DMC method - enhanced by the importance sampling scheme described in (29) - remains uniformly bounded in time, regardless of the value of $A$ and $S$. However, there's a trade-off involved: we will only have access to the measures at specific times. Indeed, despite the converging property that we are about to prove, the sequence of empirical measures $\hat{\eta}_n^{(k)}$ only approximates the measures $\hat{\eta}_l$ for $l \in k \mathbb{N}$.

Before proving the main result of this section, we conduction a short and necessary study of the functions $G^{(k)}$.

**Lemma 6.** There exists a sequence of positive definite matrices $S_n$ such that

$$ G^{(n)}(x) \propto \exp \left( -\frac{1}{2} x' S_n x \right). $$

Moreover, there exists $S_{inf}, S_{sup}$ symmetric positive definite such that $S_{inf} < S_n < S_{sup}$.

**Proof :**
Using the same notations as the one used in the proof of Proposition 3, the Markov operators $\bar{Q}_{0,n} := Q_{0,n}/Q_{0,n}(1)$ are then given by

$$ \delta_x \bar{Q}_{0,n} = \phi_{0,n}(\delta_x) = \mathcal{N}(A_n x, B_n) \quad (39) $$

with

$$ A_n := \mathcal{E}_n(0) \quad \text{and} \quad B_n = \Phi^n(0). $$

Using the product formula

$$ Q_{0,n}(1)(x) = \prod_{0 \le l < n} \phi_{0,l}(\delta_x)(G), $$

we check that

$$ Q_{0,n}(1)(x) \propto \exp \left( -\frac{1}{2} x' S_n x \right) $$

with

$$ \begin{aligned} S_n &:= \sum_{0 \le l < n} \mathcal{E}_l(0)' (S - S (\Phi^l(0)^{-1} + S)^{-1} S) \mathcal{E}_l(0) \\ &= \sum_{0 \le l < n} \mathcal{E}_l(0)' (S^{-1} + \Phi^l(0))^{-1} \mathcal{E}_l(0) \ge 0 \end{aligned} $$

The matrices $\mathcal{E}_l(0)' (S^{-1} + \Phi^l(0))^{-1} \mathcal{E}_l(0)$ being positive for all $l \in \mathbb{N}^*$, $S_n$ is lower bounded by $S_1$.

Moreover, recalling that $S > 0$ and using the results from Theorem 1.3. and Corollary 1.4 in [53] as well as the proof of Proposition 2.3, we know that there exists a positive definite matrix $P_\infty$, a matrix $E$ such that $\lambda_{max}(E) < 1$ as well as a uniformly bounded sequence of matrices $L_n$ such that we have

$$ \mathcal{E}_n(0) = E^n L_n \quad \text{and} \quad \Phi^n(0) = P_\infty - E^n L_n P_\infty(E^n)' $$

The sequence $\mathcal{E}_l(0)' (S - S (\Phi^l(0)^{-1} + S)^{-1} S) \mathcal{E}_l(0)$ is then converging, thus there exists $C_1 \in \mathbb{R}$ such that

$$ S_n \le C_1 \sum_{0 \le l < n} \mathcal{E}_l(0)' \mathcal{E}_l(0) $$

Since $E^n$ is converging to zero, there exists $n_1 \in \mathbb{N}$ such that $\|E^{n_1}\| < 1$, with $\|\cdot\|$ being the spectral norm. Let $C_2 = C_1 \times \max_{k \in \llbracket 0, n_1 \rrbracket} \|E^k\| \times \sup_{k \in \mathbb{N}} \|L_k\|$. The spectral norm being sub-multiplicative, we have then

$$ \|S_n\| \le C_2 \sum_{0 \le l < n} \|E^{n_1 \lfloor l/n_1 \rfloor}\| \le n_1 C_2 \sum_{0 \le l < \lfloor n/n_1 \rfloor} \|E^{n_1}\|^l. $$
The norm of $S_n$ being bounded, we can conclude that there exists a positive definite matrix $S_{sup}$ such that

$$ \forall n \in \mathbb{N}, S_n \leq S_{sup} $$

□

The proof of Corollary 4 is now relatively straightforward.

**Proof of Corollary 4 :**

From Corollary 3, Lemma 6 and (39), we deduce that to prove the uniform convergence toward the Feynman-Kac measure of the DMC method associated with $G^{(n)}$ and $P^{(n)}$, it suffices to prove that

$$ A'_{n}S_{n}A_{n} < S_{n}. $$

Let $\lambda_{max}(S_n)$ be the greatest eigenvalue of $S_n$. The matrix $S - \lambda_{max}(S_n)I_d$ being negative, we get from Sylvester's law of inertia that for any non-singular matrix $A$, $A'(S - \lambda_{max}(S_n)I_d)A$ is negative. By density of the invertible matrices, this holds for any matrix $A$. Thus

$$ A'_{n}S_{n}A_{n} \leq \lambda_{max}(S_{n})A'_{n}A_{n}. $$

Denoting by $||M||_2$ be the spectral norm, we have then

$$ A'_{n}S_{n}A_{n} \leq \lambda_{max}(S_{n}) ||A_{n}||_{2}^{2} I. $$

From Lemma 6 and its proof, we deduce that $A'_{n}S_{n}A_{n}$ is converging to zero. Thus, for $n$ large enough, we have

$$ A'_{n}S_{n}A_{n} < S_{inf} < S_{n}. $$

This concludes the proof.

□

We conclude this section with an observation regarding the ground state energy associated to the coupled harmonic oscillator for any value of $A$, $B$ and $S$. Notice first that all of the hypotheses required for the existence of an eigen-triple (12) are trivially met by taking $V: x \in \mathbb{R}^d \mapsto x'x + 1$. Furthermore, with $x_0$ being the null vector, it is possible to express the triple as

$$ h(x) \propto h_{\infty}(x) := \exp \left( -\frac{1}{2} x' S_{\infty} x \right) , \quad E_0 = \frac{1}{\sqrt{\det(I_d + BS_{\infty})}} $$

and

$$ \eta_{\infty} \sim \mathcal{N}(0, P_{\infty}), $$
with

$$S_{\infty} := \lim_{n \to \infty} S_n = \sum_{l=0}^{\infty} \varepsilon_l(0)' (S^{-1} + \Phi^l(0))^{-1} \varepsilon_l(0)$$

and

$$P_{\infty} = \Phi^{\infty}(0).$$

Indeed, we have

$$\begin{align}
h_n(x) &:= \exp\left(-\frac{1}{2}x'S_nx\right) = \frac{Q_{0,n}(1)(x)}{Q_{0,n}(1)(x_0)} \nonumber \\
&\Rightarrow Q(h_n) = \frac{Q_{0,n+1}(1)(x_0)}{Q_{0,n}(1)(x_0)} h_{n+1}. \tag{40}
\end{align}$$

Using the calculations from (38) we deduce that

$$\begin{align*}
Q_{0,n+1}(1)(x_0) &= Q[Q_{0,n}(1)](x_0) \\
&= Q_{0,n}(1)(x_0) \times Q(h_n)(x_0) \\
&= \frac{Q_{0,n}(1)(x_0)}{\sqrt{\det(I_d + BS_n)}}.
\end{align*}$$

Hence, going to the limit in (40) we obtain

$$Q(h_{\infty}) = \frac{1}{\sqrt{\det(I_d + BS_{\infty})}} h_{\infty}.$$

The value of $\eta_{\infty}$ directly follows from the computations of the flow of Gaussian measures $\eta_n$ in the proof of Proposition 3.

It is also worth noting that, by analogy with the calculations done in the proof of Lemma 5, $S_{\infty}$ is the solution to the Riccati equation

$$S_{\infty} - A'(S_{\infty}B + I_d)^{-1}S_{\infty}A - S = 0.$$

## 5.3 Stability w.r.t to the initial measure

So far, in the context of the coupled harmonic oscillator, we have concentrated on examining the convergence of the DMC method with a fixed initial distribution. This naturally raises the question: does the choice of the initial measure significantly affect the control of the method ? To answer this question, we establish a stability result for the Feynman-Kac flow of measure (8) in the harmonic oscillator framework. Our analysis shows that the system tends to forget its initial condition exponentially fast. Consequently, regardless of the initial condition, the flow of the Feynman-Kac measure converges exponentially to the stationary measure $\eta_{\infty}$.
**Proposition 4.** Let S, B be two symmetric positive definite matrices on $\mathbb{R}^d$ and A be a matrix. Consider $G = G_S$ and $P = P_{A,B}$. Let $\nu_1$ and $\nu_2$ be two probability measures on $\mathbb{R}^d$. There exists $(a, b) \in \mathbb{R} \times [0, 1)$ such that

$$ \|\phi_n(\nu_1) - \phi_n(\nu_2)\|_{tv} \le ab^n, $$

with the total variation distance between the measures $\nu_1$ and $\nu_2$ defined by

$$ \|\nu_1 - \nu_2\|_{tv} := \sup_{\|f\|_{\infty} \le 1} |\nu_1(f) - \nu_2(f)| $$

**Proof:**

Let $f$ such that $\|f\|_{\infty} \le 1$, then

$$ \begin{aligned} & \phi_n(\nu_1)(f) - \phi_n(\nu_2)(f) \\ &= \frac{\nu_1 Q_n(f)}{\nu_1 Q_n(1)} - \frac{\nu_2 Q_n(f)}{\nu_2 Q_n(1)} \\ &= \psi_{Q_n(1)}(\nu_1) \bar{Q}_n(f) - \psi_{Q_n(1)}(\nu_2) \bar{Q}_n(f) \\ &= \int_{\mathbb{R}^d} \psi_{Q_n(1)}(\nu_1)(dx) \psi_{Q_n(1)}(\nu_2)(dy) [\bar{Q}_n(f)(x) - \bar{Q}_n(f)(y)]. \end{aligned} $$

We have then

$$ \begin{aligned} & |\phi_n(\nu_1)(f) - \phi_n(\nu_2)(f)| \\ & \le \int_{\mathbb{R}^d} \psi_{Q_n(1)}(\nu_1)(dx) \psi_{Q_n(1)}(\nu_2)(dy) \|\phi_n(\delta_x) - \phi_n(\delta_y)\|_{tv}. \end{aligned} \quad (41) $$

From Pinsker's inequality, we have, with $Ent$ being the relative entropy - or Kullback-Leibler divergence

$$ \|\phi_n(\delta_x) - \phi_n(\delta_y)\|_{tv} \le \sqrt{\frac{1}{2} Ent(\phi_n(\delta_x) | \phi_n(\delta_y))}. $$

Using the results from [54], we know that for two Gaussian distributions $\mu_1 = \mathcal{N}(m_1, \Sigma_1^2)$ and $\mu_2 = \mathcal{N}(m_2, \Sigma_2^2)$, we have

$$ \begin{aligned} Ent(\mu_1|\mu_2) = \frac{1}{2}[\log|\Sigma_2\Sigma_1^{-1}| + \operatorname{tr}(\Sigma_2^{-1}\Sigma_1) \\ + \{(m_1 - m_2)'\Sigma_2^{-1}(m_1 - m_2) - n\}]. \end{aligned} $$

Thus, using the notations from the proof of Property 3, $\phi_n(\delta_x)$ and $\phi_n(\delta_y)$ being Gaussians with the same covariance $\Omega_n$ and respective mean $\mathcal{E}_n(0)x$ and $\mathcal{E}_n(0)y$, we deduce that
$$Ent(\phi_n(\delta_x) | \phi_n(\delta_y)) = (x - y)' \mathcal{E}_n(0)' \Omega_n \mathcal{E}_n(0)(x - y).$$

From Theorem 1.3. and Corollary 1.4 in [53], we know that $\Omega_n$ is uniformly bounded in time and that $\mathcal{E}_n(0)$ converges exponentially fast to zero. Thus, there exists $(a, b) \in \mathbb{R} \times [0, 1)$ such that

$$Ent(\phi_n(\delta_x) | \phi_n(\delta_y)) \leq ab^n \|x - y\|^2.$$

Combining this with Pinsker's inequality and (41), we get, for some $(a, b) \in \mathbb{R} \times [0, 1)$

$$|\phi_n(\nu_1)(f) - \phi_n(\nu_2)(f)| \leq ab^n[\psi_{Q_n(1)}(\nu_1)(||\cdot||) + \psi_{Q_n(1)}(\nu_2)(||\cdot||)],$$

where $||\cdot||$ represents the function that maps $x \in \mathbb{R}^d$ to its norm. Using the computations from Lemma 6, we know that

$$Q_n(1) \propto \exp\left(-\frac{1}{2} x'S_n x\right) := q_n,$$

and that there exists two symmetric positive definite matrices $S_{\text{sup}}$ and $S_{\text{inf}}$ such that for any $n \in \mathbb{N}$, $S_{\text{inf}} \leq S_n \leq S_{\text{sup}}$.

Considering the function $q_{\text{inf}} := \exp(-\frac{1}{2} x'S_{\text{inf}}x)$ and $q_{\text{sup}} := \exp(-\frac{1}{2} x'S_{\text{sup}}x)$, we have then

$$\Psi_{Q_n(1)}(\nu_1)(||\cdot||) = \frac{\nu_1(||\cdot||q_n)}{\nu_1(q_n)} \leq \frac{\nu_1(||\cdot||q_{\text{inf}})}{\nu_1(q_{\text{sup}})}.$$

The function $||\cdot||_{\text{inf}}$ and $q_{\text{sup}}$ being bounded, we deduce that the previous upper bound is finite.

Hence, for any $f$ such that $||f||_{\infty} \leq 1$, we have

$$|\phi_n(\nu_1)(f) - \phi_n(\nu_2)(f)| \leq ab^n \left[ \frac{\nu_1(||\cdot||q_{\text{inf}})}{\nu_1(q_{\text{sup}})} + \frac{\nu_2(||\cdot||q_{\text{inf}})}{\nu_2(q_{\text{sup}})} \right].$$

We can thus conclude.

□

## 5.4 Divergence and fluctuation estimates

In the previous subsections, we presented a simple sufficient condition for controlling the DMC method and introduced an importance sampling technique that satisfies this criterion. However, it is natural to question the robustness of this condition and whether it is necessary to use importance sampling. Specifically, for the multi-dimensional harmonic oscillator, the convergence condition reduced to $A^2 < 1$, and we will prove divergence of the DMC method when this stability condition is not met.
Within this framework, we can break down the evolution of the walkers into two distinct steps, a mutation transition and a selection transition

$$ (\xi_n^i)_{i \in [1,N]} \in \mathbb{R}^N \xrightarrow{\text{selection}} (\hat{\xi}_n^i)_{i \in [1,N]} \in \mathbb{R}^N \\ \xrightarrow{\text{mutation}} (\xi_{n+1}^i)_{i \in [1,N]}. $$

The initial configuration $(\xi_0^i)_{i \in [1,N]}$ is determined by sampling $N$ independent random variables from the distribution $\eta_0$. The selection transition involves the sampling of $N$ independent random variables $(\hat{\xi}_n^i)_{i \in [1,N]}$ using the weighted distributions

$$ \epsilon_n(\eta_n^N) G_S(\xi_n^i) \delta_{\xi_n^i} + (1 - \epsilon_n(\eta_n^N) G_S(\xi_n^i)) \sum_{k \in [1,N]} \frac{e^{-\frac{S}{2}\xi_n^{k^2}}}{\sum_{j \in [1,N]} e^{-\frac{S}{2}\xi_n^{j^2}}} \delta_{\xi_n^k}. $$

The mutation transition is defined using a family of Gaussian random variables with zero-mean and unit variance $(V_n^i)_{i \in [1,N]}$ such that

$$ \xi_n^i = A \hat{\xi}_{n-1}^i + \sqrt{B} V_n^i. $$

The measures $(\eta_n)_{n \in \mathbb{N}}$ can be described exhaustively using the Kalman filter equations. It provides us with the mean and variances $(m_n, \sigma_n^2)$ of the Gaussian random variables $\eta_n$ with the recurrent equations

$$ \left\{ \begin{array}{l} m_{n+1} = \frac{A}{1 + S\sigma_n^2} m_n \\ \\ \sigma_{n+1}^2 = \frac{A^2\sigma_n^2}{1 + S\sigma_n^2} + B. \end{array} \right. \qquad (42) $$

In this scenario, when the condition $A^2 > 1$ is met, it is possible to prove that the DMC's error in approximating the Feynman-Kac measure does not admit a uniform-in-time bound. It is properly stated in Property 1, and we can now conduct its proof.

**Proof of Proposition 1:**

For any $n \ge 2$, we know from (42) that

$$ \eta_n(I) := m_n = \frac{A^2}{(1 + S\sigma_{n-1}^2)(1 + S\sigma_{n-2}^2)} m_{n-2} \le A^2 m_{n-2}. $$

For any $n \in \mathbb{N}^*$, let $\xi_n^* = \min_{i \in [1,N]} \xi_n^i$ and define the random variables $V_n^*$ in the following way:
$$V_n^* = - \max_{i \in [1, N]} |V_n^i|.$$

By definition of the evolution of the walkers, there exists $(i, j) \in [1, N]$ such that

$$\begin{aligned} \xi_{2n}^* &= A^2 \xi_{2n-2}^j + \sqrt{B} V_{2n-1}^i + A \sqrt{B} V_{2n-2}^j \\ &\ge A^2 \xi_{2n-2}^* + \sqrt{B} V_{2n-1}^* + |A| \sqrt{B} V_{2n-2}^* \end{aligned}$$

Thus

$$\eta_{2n}^N(I) - \eta_{2n}(I) \ge A^2(\xi_{2(n-1)}^* - m_{2(n-1)}) + \sqrt{B} V_{2n-1}^* + |A| \sqrt{B} V_{2n-2}^*.$$

Iterating the process, we obtain

$$\frac{\eta_{2n}^N(I) - \eta_{2n}(I)}{A^{2n}} \ge (\xi_0^* - m_0) + \sqrt{B} \sum_{1 \le k \le n} \frac{V_{2k-1}^*}{A^{2k}} \\ + |A| \sqrt{B} \sum_{1 \le k \le n} \frac{V_{2(k-1)}^*}{A^{2k}}.$$

For any sequence of $N$ independent centred Gaussian random variables $U_i$ with unit variance, we have

$$\mathbb{E} \left[ \max_{1 \le i \le N} |U_i| \right] \le \sqrt{2 \log(2N)}.$$

This inequality is obtained by using Jensen's inequality as follows, with $t = \sqrt{2 \log(2N)}$

$$\begin{aligned} \exp \left[ t \mathbb{E} \left( \max_{1 \le i \le N} |U_i| \right) \right] &\le \mathbb{E} \left[ \exp \left( t \max_{1 \le i \le N} |U_i| \right) \right] \\ &\le \sum_{i=1}^{N} \mathbb{E} \left[ \exp (t |U_i|) \right], \end{aligned}$$

and noticing that

$$\begin{aligned} \mathbb{E} [\exp (t |U_i|)] &= 2 \int_0^{+\infty} \exp \left( - \frac{(x-t)^2 + t^2}{2} \right) dx \\ &\le 2 \exp(t^2/2). \end{aligned}$$
Then, on the event

$$ \Omega_{\epsilon} := \left\{ \xi_{0}^{*} \geq \epsilon + m_{0} + \frac{2\sqrt{B}(1 + |A|)}{A^{2} - 1} \sqrt{2 \log(2N)} \right\}, \quad (43) $$

with $\epsilon \in \mathbb{R}_{+}^{*}$, we have

$$ \mathbb{E}[\eta_{2n}^{N}(I) - \eta_{2n}(I)|\xi_{0}^{*}] \geq \epsilon A^{2n} \xrightarrow{n \to +\infty} +\infty. \quad (44) $$

Integrating over $\xi_{0}^{*}$ we deduce

$$ \mathbb{E}[|\eta_{2n}^{N}(I) - \eta_{2n}(I)|] \geq \epsilon A^{2n} \mathbb{P}(\Omega_{\epsilon}). \quad (45) $$

We can then conclude by noticing

$$ \mathbb{P}(\Omega_{\epsilon}) = \eta_{0} \left\{ \left[ \epsilon + X_{0} + \frac{2\sqrt{B}(1 + |A|)}{A^{2} - 1} \sqrt{2 \log(2N)}, +\infty \right) \right\}^{N} > 0. $$

□

For the case $A = 1$, we are unable to determine whether a uniform bound exists. To the best of our knowledge, the strongest divergence-type result established so far is a linear bound on the variance of the unnormalized measure for $R = S = 1$ obtained in [55]. In that work, the authors employed a stronger exponential drift condition compared to our linear condition (1). As a result, they used a quadratic Lyapunov function, whereas we adopted an exponential Lyapunov function, which enabled us to apply Corollary 3 to a broader range of test functions.

The divergence result in Proposition 1 highlights the importance of studying non-asymptotic uniform convergence results rather than relying solely on central limit theorems (CLTs). Extensive research has been devoted to CLTs and, under appropriate assumptions. We quote the first studies in this field [25, 27] mainly based on uniformly bounded potential and test functions.

More general fluctuation theorems that apply to more general models including diffusion-type processes with Lipschitz drift and diffusion functions as well as test functions with at most polynomial growth are discussed in [56, 57]. We can formulate the following result :

**Lemma 7.** In the context of proportional selection/reconfiguration, for any $n \geq 1$, we have the following convergence in law as $N$ tends toward $+\infty$

$$ \sqrt{N}[\eta_n(I) - \eta_n^N(I)] \xrightarrow{\mathcal{L}} \mathcal{N}(0, \sigma_n^2), $$
with the asymptotic variance

$$ \sigma_n^2 := \sum_{p=0}^{n} \eta_p \left[ \left( \frac{Q_{p,n}(1)}{\eta_p Q_{p,n}(1)} \bar{Q}_{p,n}(I - \eta_n(I)) \right)^2 \right]. $$

Related asymptotic variance formulae and comparisons are discussed in [49] in the context of random walks with absorbing barriers, including geometric killing rates and local reflection moves.

In [44], Nick Witheley also obtains a uniform time bound on the asymptotic variance. However, our next proposition shows that, for a given system, the DMC method can have a uniform time bound on its asymptotic variance, despite the fact that its non-asymptotic variance is unbounded.

**Proposition 5.** Let $G: x \in \mathbb{R} \mapsto e^{-\frac{x^2}{2}} S$, and $P$ such that $\delta_x P \sim N(Ax, B)$ for some $(A, B, S) \in (1, +\infty] \times \mathbb{R}_+^2$. There exists an initial distribution $\eta_0$ such that

$$ \sup_{n \in \mathbb{N}} \sigma_n^2 < \infty \quad \text{and} \quad \sup_{n \in \mathbb{N}} \mathbb{E} \left( |\eta_n(I) - \eta_n^N(I)|^p \right)^{\frac{1}{p}} = +\infty, \quad (46) $$

where $\sigma_n^2$ is defined as in Lemma 7.

**Proof:**

Considering $\eta_0 \sim N(m_0, \sigma_0^2)$ and using the Kalman filter's equations, we are able to fully describe the measures $\eta_n$. For $n \in \mathbb{N}^*$ we have

$$ \eta_n \sim N(m_n, \sigma_n^2) \quad \text{with} \quad \left\{ \begin{array}{l} m_{n+1} = \frac{A}{1 + S\sigma_n^2} m_n \\ \sigma_{n+1}^2 = \frac{A^2}{1 + S\sigma_n^2} + B. \end{array} \right. $$

The limit measure is given by $\eta_\infty = N(0, \sigma_\infty^2)$ where $\sigma_\infty^2$ is the fixed point of the function $x \in \mathbb{R}_+ \mapsto \frac{A^2}{1+Sx} + B$. From now on, we assume that $\eta_0 = \eta_\infty$. We have then

$$ \sigma_n^2 = \sum_{p=0}^{n} \frac{\eta_\infty (Q_{p,n}(I)^2)}{(\eta_\infty Q_{p,n}(1))^2}. $$

Using the calculations and notations from Lemma 10, we have then

$$ \begin{aligned} \eta_\infty (Q_{p,n}(1)) &= \frac{\lambda_{n-p}}{\sqrt{2\pi\sigma_\infty^2}} \int_{\mathbb{R}} e^{-\frac{y^2}{2\sigma_\infty^2}} e^{-\frac{y^2}{2}q_{n-p}dy} \\ &= \frac{\lambda_{n-p}}{\sqrt{q_{n-p}\sigma_\infty^2 + 1}}, \end{aligned} $$
as well as

$$ 
\begin{aligned}
\eta_{\infty} (Q_{p,n}(I)^2) &= \frac{\mu_{n-p}^2}{\sqrt{2\pi}\sigma_{\infty}^2} \int_{\mathbb{R}} y^2 e^{-\frac{y^2}{2\sigma_{\infty}^2}} e^{-y^2 q_{n-p}} dy \\
&= \frac{\mu_{n-p}^2 \sigma_{\infty}^2}{(2q_{n-p}\sigma_{\infty}^2 + 1)^{3/2}}.
\end{aligned}
$$

Then

$$ \sigma_n^2 = \sigma_\infty^2 \sum_{p=0}^{n} \frac{\mu_{n-p}^2 (q_{n-p}\sigma_\infty^2 + 1)^{3/2}}{\lambda_{n-p}^2 (2q_{n-p}\sigma_\infty^2 + 1)^{3/2}}. $$

As a solution to a Riccati equation, $q_n$ converges towards some $q_\infty \in \mathbb{R}_+^*$ such that

$$ q_{\infty} = \frac{A^2 q_{\infty}}{1 + q_{\infty} B} + S. $$

Since $A > 1$, we obtain that $1 + q_\infty B = A^2 + S/q_\infty > A$. Thus, there exists $C \in (0, 1)$ such that, for $n$ large enough, we have

$$ \mu_{n+1} \leq C \frac{\mu_n}{\sqrt{1 + q_n B}}. $$

Thus, comparing the definition of $\lambda_n$ and the previous bound, we deduce that there exists $\alpha \in \mathbb{R}$ such that we have

$$ \frac{\mu_{n-p}^2 (q_{n-p}\sigma_\infty^2 + 1)^{3/2}}{\lambda_{n-p}^2 (2q_{n-p}\sigma_\infty^2 + 1)^{3/2}} \leq \alpha C^{2(n-p)}. $$

We can thus conclude on the right-hand side. of (46). The left-hand side is a direct application of Proposition 1. □

## Acknowledgments

MC thanks the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement no. 863481) for financial support.

Luc de Montella would like to thank Naval Group for their financial support through the CIFRE PhD program.

We also thank the anonymous reviewer for her/his excellent suggestions for improving the paper. Her/his detailed comments greatly improved the presentation of the article.
# Appendix

## Some technical Lemmas

**Lemma 8.** Let $(G, V) \in \mathcal{C}(E) \times \mathcal{C}_{\infty}(E)$ be positive functions and $K \subset E$ be such that the right-hand side of (6) holds. There exists $c \in \mathbb{R}_+$ such that for any probability measure $\mu$ on $E$:

$$ \psi_G(\mu)(V) \leq \mu(V) + c. $$

**Proof:** Let $\bar{V} := 1_{\mathbb{R}\setminus K}V$. For any $(x, t) \in E^2$, we have:

$$ \begin{aligned} & 0 ge(G(x) - G(y))(\bar{V}(x) - \bar{V}(y)) \\ & = G(x)\bar{V}(x) + G(y)\bar{V}(y) - G(y)\bar{V}(x) - G(x)\bar{V}(y). \end{aligned} $$

Integrating with respect to the probability measure $\mu$ over both $x$ and $y$, we obtain:

$$ \mu(G\bar{V}) - \mu(G)\mu(\bar{V}) \leq 0. $$

Hence

$$ \psi_G(\mu)(\bar{V}) - \mu(\bar{V}) = \frac{\mu(G\bar{V}) - \mu(G)\mu(\bar{V})}{\mu(G)} \leq 0. $$

Thus

$$ \begin{aligned} \psi_G(\mu)(V) \leq \psi_G(\mu)(\bar{V}) + \sup_K V &\leq \mu(\bar{V}) + \sup_K V \\ &\leq \mu(V) + \sup_K V. \end{aligned} $$

□

**Lemma 9.** Let $V \in \mathcal{C}_{\infty}(E)$ be a Q-Lyapunov function, then

$$ \left\{ \begin{array}{l} \sup_{(q,n,N) \in \mathbb{N}^3} \mathbb{E} \left[ \phi_{q,n}(\eta_q^N)(V) \right] < +\infty \\ \\ \sup_{(q,n) \in \mathbb{N}^2} \phi_{q,n}(\eta_q)(V) < +\infty. \end{array} \right. \qquad (47) $$

**Proof:** First, let's prove that for any $\gamma \in \mathcal{P}_V(E)$, we have, with $(\epsilon, c) \in [0, 1) \times \mathbb{R}$ defined in (1),

$$ \sup_{\substack{(q,n) \in \mathbb{N}^2 \\ q \le n}} \phi_{q,n}(\gamma)(V) \le \gamma(V) + \frac{c}{1-\epsilon}. \qquad (48) $$
To do so, we begin by using the Q-Lyapunov property of V as well as Lemma 8 for some $l \in \mathbb{N}^*$ to deduce that

$$ \phi_{q,n}(\gamma)(V) = \frac{\phi_{q,n-1}(\gamma)(G_n P_n(V))}{\phi_{q,n-1}(\gamma)(G_n)} \le \epsilon \phi_{q,n-1}(\gamma)(V) + c'. $$

By iterating the process, we obtain (48). We now prove that

$$ \sup_{n \in \mathbb{N}} \mathbb{E} [\eta_n^N(V)] \le \eta_0(V) + \frac{c'}{1-\epsilon}. \quad (49) $$

Notice first that

$$ \begin{aligned} \mathbb{E} [\eta_n^N(V) | \eta_{n-1}^N] &= \frac{1}{N} \sum_{i=1}^N \eta_{n-1}^N S_{n-1, \eta_{n-1}^N} P_n(V) \\ &= \psi_{G_{n-1}}(\eta_{n-1}^N)(P_n(V)). \end{aligned} $$

Then, using as previously the Q-Lyapunov property of V and Lemma 8 we obtain

$$ \psi_{G_{n-1}}(\eta_{n-1}^N)(P_n(V)) \le \epsilon \eta_{n-1}^N(V) + c'. $$

By iterating the process, we obtain

$$ \mathbb{E} [\eta_n^N(V)] \le \epsilon^n \eta_0(V) + c' \sum_{i=0}^{n-1} \epsilon^i \le \eta_0(V) + \frac{c'}{1-\epsilon}. $$

Thus, by combining (48) and (49), we can conclude regarding the first part of (47). The second part is obtained by proceeding in a strictly analogous way. $\square$

## Proof of Lemma 1

We first consider the case where only (20) holds and prove that

$$ \phi_{0,q}(\gamma)(Q_{q,n}(1)) \le C \chi^{n-q-1}. \quad (50) $$

For $q \le n + 2$, we have

$$ \begin{aligned} Q_{q,n}(1)(x) &= Q_{q,n-2}(G_{n-1}P_n(G_n))(x) \\ &\le \chi Q_{q,n-2}(G_{n-1})(x) \le \chi Q_{q,n-1}(1)(x). \end{aligned} $$

Iterating the process, we deduce
$$Q_{q,n}(1) \leq \chi^{n-q-1} G_{q-1}.$$

Using (19) we deduce (50).

For $q \geq 1$ and $\mu \in \{\eta_q^N, \phi_q(\eta_{q-1}^N)\}$, we have then

$$\begin{aligned} \frac{1}{\mu(H_{q,n}^{\gamma})} &= \phi_{q,n}(\mu)(W) \times \frac{\phi_{0,q}(\gamma)Q_{q,n}(1)}{\mu(Q_{q,n}(W))} \\ &\leq C \frac{\chi^{n-q-1}\phi_{q,n}(\mu)(W)}{\mu(Q_{q,n}(W))}. \end{aligned}$$

From Holder's inequality, Jensen's inequality and the hypothesis on $W$, we get:

$$\begin{aligned} \mathbb{E}[\mu(H_{q,n})^{-\beta}] &\leq C^{\beta} \mathbb{E}[\phi_{q,n}(\mu)(W)^{2\beta}]^{\frac{1}{2}} \\ &\quad \times \mathbb{E}[\chi^{2\beta(n-q-1)}\mu(Q_{q,n}(W))^{-2\beta}]^{\frac{1}{2}} \\ &\leq \frac{C^{\beta}}{\chi^{\beta}} \mathbb{E}[\phi_{q,n}(\mu)(V)^{2\beta}]^{\frac{1}{2}} \mathbb{E}[\mu(W^{-2\beta})]^{\frac{1}{2}}. \end{aligned}$$

Then, by choosing $\beta$ small enough such that $W^{-2\beta} \leq \bar{W}$ and $V^{2\beta} \leq \bar{V}$, where $\bar{W}$ and $\bar{V}$ are $Q$-Lyapunov functions, we obtain

$$\begin{aligned} \sup_{\substack{(q,n,N) \in \mathbb{N}^3 \\ q \leq n}} \mathbb{E}[\mu(H_{q,n}^{\gamma})^{-\beta}] &\leq \frac{C^{\beta}}{\chi^{\beta}} \sup_{\substack{(q,n,N) \in \mathbb{N}^3 \\ q \leq n}} \mathbb{E}[\mu(\bar{V})]^{\frac{1}{2}} \\ &\quad \times \sup_{\substack{(q,n,N) \in \mathbb{N}^3 \\ q \leq n}} \mathbb{E}[\mu(\bar{W})]^{\frac{1}{2}}. \end{aligned}$$

We can conclude using Lemma 9.

The demonstration for time-homogeneous models with $E_0 < \chi$ is analogous. Indeed the equivalent of (50) is obtained by noticing that

$$\begin{aligned} \phi_{0,q}(\eta_{\infty})Q_{q,n}(1) &= \eta_{\infty}Q_{q,n-1}(G) = \eta_{\infty}Q_{q,n-1}(1)\eta_{\infty}(G) \\ &= E_0 \times \eta_{\infty}Q_{q,n-1}(1) = E_0^{n-q}. \end{aligned}$$

The rest of the proof follows the same arguments, thus it is skipped. This ends the proof of the lemma. $\square$
**Lemma 10.** Let $G: x \in \mathbb{R} \mapsto e^{-\frac{x^2}{2}}S$, and $P$ such that $\delta_x P \sim \mathcal{N}(Ax, B)$ for some $(A, B, S) \in \mathbb{R} \times \mathbb{R}_+^{*2}$. Denoting by $Q^n$ the n-times composition of the operator $Q$, we have

$$Q^n(1)(x) = \lambda_n e^{-\frac{x^2}{2}q_n} \quad \text{and} \quad Q^n(I)(x) = \mu_n x e^{-\frac{x^2}{2}q_n}, \qquad (51)$$

with the parameters $(q_0, \lambda_0, \mu_0) = (S, 1, A)$ and

$$\left\{ \begin{array}{l} q_{n+1} = \frac{A^2 q_n}{1 + q_n B} + S \\ \\ \lambda_{n+1} = \frac{\lambda_n}{\sqrt{1 + q_n B}} \quad \text{and} \quad \mu_{n+1} = \frac{A \mu_n}{(1 + q_n B)^{3/2}}. \end{array} \right.$$

**Proof :**

For $n = 0$ the result is immediate. Assume that it holds for some $n \in \mathbb{N}$. In this situation, we have

$$\begin{aligned} Q^{n+1}(1)(x) &= Q(Q^n(1))(x) = \lambda_n Q \left[ y \mapsto e^{-\frac{y^2}{2}q_n} \right] (x) \\ &= \frac{\lambda_n}{\sqrt{2\pi B}} e^{-\frac{x^2}{2}S} \int_{\mathbb{R}} e^{-\frac{y^2}{2}q_n - \frac{(Ax-y)^2}{2B}} dy \\ &= \frac{\lambda_n}{\sqrt{2\pi B}} e^{-\frac{x^2}{2} \left( \frac{A^2}{B} - \frac{A^2}{B(1+q_n B)} + S \right)} \\ & \quad \times \int_{\mathbb{R}} e^{-\frac{(1+q_n B) \left( \frac{A}{1+q_n B} x - y \right)^2}{2B}} dy \\ &= \frac{\lambda_n}{\sqrt{1+q_n B}} e^{-\frac{x^2}{2} \left( \frac{A^2 q_n}{1+q_n B} + S \right)}. \end{aligned}$$

Similarly, we have

$$\begin{aligned} Q^{n+1}(I)(x) &= Q(Q^n(I))(x) = \mu_n Q \left[ y \mapsto ye^{-\frac{y^2}{2}q_n} \right] (x) \\ &= \frac{\mu_n}{\sqrt{2\pi B}} e^{-\frac{x^2}{2} \left( \frac{A^2 q_n}{1+q_n B} + S \right)} \\ & \quad \times \int_{\mathbb{R}} ye^{-\frac{(1+q_n B) \left( \frac{A}{1+q_n B} x - y \right)^2}{2B}} dy \\ &= \frac{A\mu_n}{(1+q_n B)^{3/2}} xe^{-\frac{x^2}{2} \left( \frac{A^2 q_n}{1+q_n B} + S \right)}. \end{aligned}$$

This ends the proof.

<figure><img src="image_2.png" /></figure>
References

[1] Nicholas Metropolis and Stanislaw Ulam. The Monte Carlo method. *Journal of the American statistical association*, 44(247):335-341, 1949.

[2] Monroe D. Donsker and Mark Kac. A sampling method for determining the lowest eigenvalue and the principal eigenfunction of schrödinger's equation. *Journal of Research of the National Bureau of Standards*, 44(50):551-557, 1950.

[3] Malvin H Kalos. Monte Carlo calculations of the ground state of three-and four-body nuclei. *Physical Review*, 128(4):1791, 1962.

[4] Malvin H. Kalos D. M. Ceperley. *Monte Carlo Methods in Statistical Physics*. Springer-Verlag Berlin, 1986.

[5] Peter J. Reynolds, David M. Ceperley, Berni J. Alder, and Jr. Lester, William A. Fixed-node quantum Monte Carlo for molecules. *The Journal of Chemical Physics*, 77(11):5593-5603, 12 1982.

[6] Michel Caffarel and Pierre Claverie. Treatment of the Schrödinger equation through a Monte Carlo method based upon the generalized Feynman-Kac formula. *Journal of Statistical Physics*, 43:797-801, 1986.

[7] Michel Caffarel and Pierre Claverie. Development of a pure diffusion quantum Monte Carlo method using a full generalized Feynman-Kac formula. i. formalism. *The Journal of chemical physics*, 88(2):1088-1099, 1988.

[8] Jack H. Hetherington. Observations on the statistical iteration of matrices. *Physical Review A*, 30(5):2713, 1984.

[9] Sandro Sorella. Green function Monte Carlo with stochastic reconfiguration. *Physical review letters*, 80(20):4558, 1998.

[10] Sandro Sorella and Luca Capriotti. Green function Monte Carlo with stochastic reconfiguration: An effective remedy for the sign problem. *Physical Review B*, 61(4):2599, 2000.

[11] Roland Assaraf, Michel Caffarel, and Anatole Khelif. Diffusion Monte Carlo methods with a fixed number of walkers. *Physical Review E*, 61(4):4566, 2000.

[12] Baroni Stefano and Moroni Saverio. Reptation Quantum Monte Carlo: A method for unbiased ground-state averages and imaginary-time correlations. *Phys. Rev. Lett.*, 82:4745, 1999.

[13] Michel Mareschal. The early years of quantum Monte Carlo (1): the ground state. *The European Physical Journal H*, 46(1):11, 2021.
[14] Yao Ma, Lu Meng, Yan-Ke Chen, and Shi-Lin Zhu. Ground state baryons in the flux-tube three-body confinement model using diffusion Monte Carlo. *Physical Review D*, 107(5):054035, 2023.

[15] Matthew WMC Foulkes, Lubos Mitas, RJ Needs, and Guna Rajagopal. Quantum Monte Carlo simulations of solids. *Reviews of Modern Physics*, 73(1):33, 2001.

[16] Pierre Del Moral and Laurent Miclo. Particle approximations of Lyapunov exponents connected to Schrödinger operators and Feynman-Kac semigroups. *ESAIM: Probability and Statistics*, 7:171–208, 2003.

[17] Pierre Del Moral and Arnaud Doucet. Particle motions in absorbing medium with hard and soft obstacles. *Stochastic Analysis and Applications*, 22(5):1175–1207, 2004.

[18] Mathias Rousset. On the control of an interacting particle estimation of schrödinger ground states. *SIAM J. Math. Anal.*, 38:824–844, 2006.

[19] Eric Cancés, Benjamin Jourdain, and Tony Lelièvre. Quantum Monte Carlo simulations of fermions. A mathematical analysis of the fixed-node approximation. *Mathematical Models and Methods in Applied Sciences*, 16, 11 2011.

[20] Mohamed El Makrini, Benjamin Jourdain, and Tony Lelièvre. Diffusion monte carlo method: Numerical analysis in a simple case. *ESAIM: Mathematical Modelling and Numerical Analysis*, 41(2):189–213, 2007.

[21] Pierre Del Moral. *Feynman-Kac Formulae: Genealogical and Interacting Particle Systems With Applications*, volume 100. Springer, 2004.

[22] Pierre Del Moral. *Mean field simulation for Monte Carlo integration*. CRC press, 2013.

[23] Pierre Del Moral and Alice Guionnet. Large deviations for interacting particle systems: applications to non-linear filtering. *Stochastic processes and their applications*, 78(1):69–95, 1998.

[24] Donald A Dawson and Pierre Del Moral. Large deviations for interacting processes in the strong topology. In *Statistical modeling and analysis for complex data problems*, pages 179–208. Springer, 2005.

[25] Pierre Del Moral and Alice Guionnet. Central limit theorem for nonlinear filtering and interacting particle systems. *Annals of Applied Probability*, pages 275–297, 1999.

[26] Nicolas Chopin. Central limit theorem for sequential Monte Carlo methods and its application to Bayesian inference. *The Annals of Statistics*, 32:2385–2411, 2004.
[27] P. Del Moral and L. Miclo. *Branching and interacting particle systems approximations of Feynman-Kac formulae with applications to non-linear filtering*, pages 1–145. Springer Berlin Heidelberg, Berlin, Heidelberg, 2000.

[28] Pierre Del Moral and Alice Guionnet. On the stability of interacting processes with applications to filtering and genetic algorithms. In *Annales de l'Institut Henri Poincaré (B) Probability and Statistics*, volume 37, pages 155–194. Elsevier, 2001.

[29] Nick Whiteley. Sequential Monte Carlo samplers: error bounds and insensitivity to initial conditions. *Stochastic Analysis and Applications*, 30(5):774–798, 2012.

[30] P Del Moral and E Horton. Coupled quantum harmonic oscillators and Feynman-Kac path integrals for linear diffusive particles. *Communications in Mathematical Physics*, 402(2):2079–2127, 2023.

[31] Marc Arnaudon, Pierre Del Moral, and El Maati Ouhabaz. A lyapunov approach to stability of positive semigroups: An overview with illustrations. *Stochastic Analysis and Applications*, pages 1–80, 2023.

[32] Pierre Del Moral, Emma Horton, and Ajay Jasra. On the stability of positive semigroups. *The Annals of Applied Probability*, 33(6A):4424–4490, 2023.

[33] Malvin H. Kalos. Monte Carlo calculations of the ground state of three- and four-body nuclei. *Phys. Rev.*, 128:1791–1795, Nov 1962.

[34] Malvin H. Kalos, Dominique Levesque, and Loup Verlet. Helium at zero temperature with hard-sphere and other forces. *Phys. Rev. A*, 9:2178–2195, May 1974.

[35] Albert Messiah. *Quantum Mechanics*. Number vol. 2 in Quantum Mechanics. Elsevier Science, 1961.

[36] David J. Griffiths and Darrell F. Schroeter. *Introduction to quantum mechanics*. Cambridge University Press, Cambridge ; New York, NY, third edition edition, 2018.

[37] Gerhard Herzberg. *Molecular Spectra and Molecular Structure: Infrared and Raman of Polyatomic Molecules*. D. Van Nostrand Company, Inc., 1945.

[38] Dmitry N. Makarov. Quantum entanglement and reflection coefficient for coupled harmonic oscillators. *Phys. Rev. E*, 102:052213, Nov 2020.

[39] Artur K. Ekert. Quantum cryptography based on bell's theorem. *Phys. Rev. Lett.*, 67:661–663, Aug 1991.

[40] Elisabet Romero, Ramunas Augulis, Vladimir Novoderezhkin, Marco Ferretti, Jos Thieme, Donatas Zigmantas, and Rienk van Grondelle. Quantum coherence in photosynthesis for efficient solar energy conversion. *Nature Physics*, 10, 07 2014.
[41] Pierre Del Moral. Non linear filtering: Interacting particle solution. *Markov Processes and Related Fields*, 2:555–580, 03 1996.

[42] Rudolf E. Kálmán and Richard S. Bucy. New results in linear filtering and prediction theory. *Journal of Basic Engineering*, 83:95–108, 1961.

[43] Pierre Del Moral and Spiridon Penev. *Stochastic Processes: From Applications to Theory*. Chapman and Hall/CRC, 2017.

[44] Nick Whiteley. Stability properties of some particle filters. *The Annals of Applied Probability*, 23(6):2500 – 2537, 2013.

[45] Grégoire Ferré and Gabriel Stoltz. Error estimates on ergodic properties of discretized Feynman-Kac semigroups. 2019.

[46] P.Del Moral and L. Miclo. A Moran particle system approximation of Feynman-Kac formulae. *Stochastic Processes and their Applications*, 86(2):193–216, 2000.

[47] Marc Arnaudon and Pierre Del Moral. A duality formula and a particle Gibbs sampler for continuous time Feynman-Kac measures on path spaces. *Electronic Journal of Probability*, 25(none):1 – 54, 2020.

[48] Mathias Rousset. On the control of an interacting particle estimation of Schrödinger ground states. *SIAM Journal on Mathematical Analysis*, 38(3):824–844, 2006.

[49] Pierre Del Moral and Ajay Jasra. A note on random walks with absorbing barriers and sequential monte carlo methods. *Stochastic Analysis and Applications*, 36(3):413–442, 2018.

[50] Pierre Del Moral. Measure-valued processes and interacting particle systems. application to nonlinear filtering problems. *The Annals of Applied Probability*, 8(2):438–495, 1998.

[51] Tosio Kato. *Perturbation Theory for Linear Operators*. Springer Science & Business Media, 1966.

[52] Chi-Kwong Li and Fuzhen Zhang. Eigenvalue continuity and gersgorin's theorem. *Electronic Journal of Linear Algebra*, 35:619–625, 12 2019.

[53] Pierre Del Moral and Emma Horton. A note on riccati matrix difference equations. *SIAM Journal on Control and Optimization*, 60(3):1393–1409, 2022.

[54] Leandro Pardo. *Statistical Inference Based on Divergence Measures*. CRC Press, Abingdon, 2005.

[55] Nick Whiteley, Nikolas Kantas, and Ajay Jasra. Linear variance bounds for particle approximations of time-homogeneous Feynman-Kac formulae. *Stochastic Processes and their Applications*, 122(4):1840–1865, 2012.
[56] Pierre Del Moral and Jean Jacod. The Monte-Carlo method for filtering with discrete time observations. central limit theorems. In *Proceedings Workshop on Numerical Methods and Stochastics. The Fields Institute 1999*, editors T.J. Lyons and T.S Salisbury. American Mathematical Society, 2002.

[57] Pierre Del Moral and J. Jacod. Interacting particle filtering with discrete-time observations: asymptotic behaviour in the gaussian case. In T. Hida, R. Karandikar, H. Kunita, B. Rajput, S. Watanabe, and Xiong.J. editors, *Stochastics in finite/infinite dimensions*, pages 101–123. Trends in Mathematics, Birkhauser, 2001.