Adrien Laurent

A uniformly accurate scheme for the numerical integration of penalized Langevin dynamics

Keywords: constrained stochastic differential equations, penalized Langevin dynamics, manifolds, uniform accuracy, invariant measure AMS subject classification (2020): 60H35, 70H45, 37M25

In molecular dynamics, penalized overdamped Langevin dynamics are used to model the motion of a set of particles that follow constraints up to a parameter ε. The most used schemes for simulating these dynamics are the Euler integrator in R d and the constrained Euler integrator. Both have weak order one of accuracy, but work properly only in specific regimes depending on the size of the parameter ε. We propose in this paper a new consistent method with an accuracy independent of ε for solving penalized dynamics on a manifold of any dimension. Moreover, this method converges to the constrained Euler scheme when ε goes to zero. The numerical experiments confirm the theoretical findings, in the context of weak convergence and for the invariant measure, on a torus and on the orthogonal group in high dimension and high codimension.

Introduction

In molecular dynamics, the overdamped Langevin equation in R d is often used for modeling the behavior of a large set of particles in a high friction regime. It is given by dXptq f pXptqqdt σdW ptq, Xp0q X 0 , (1.1) where f is a smooth Lipschitz function (typically of the form f ¡∇V for V a smooth potential), σ ¡ 0 is a constant scalar, and W is a standard d-dimensional Brownian motion in R d on a probability space equipped with a filtration pF t q and fulfilling the usual assumptions. If the particles are subject to smooth constraints ζ : R d Ñ R q , such as strong covalent bonds between atoms or fixed angles in molecules, the dynamics follow the constrained overdamped Langevin equation dX 0 ptq Π M pX 0 ptqqfpX 0 ptqqdt σΠ M pX 0 ptqq ¥ dW ptq, X 0 p0q X 0 M, (1.2) where the solution lies on the manifold M tx R d , ζpxq 0u with codimension q thanks to Π M : R d Ñ R d¢d , the orthogonal projection on the tangent bundle of the manifold M.

In physical applications, constrained systems are often used as a limit model for stiff equations. For instance, in the dynamics of a diatomic molecule, the distance between the two atoms oscillates around an average length, called the bond length (see, for instance, [START_REF] Lelièvre | Free energy computations[END_REF]Sect. 1.2.1] on the interactions of particles). One can work with a simpler constrained dynamics where the distance between the atoms is fixed as a constraint, or with the original (possibly stiff) dynamics in R d . We refer the reader to [START_REF] Plecháč | Implicit mass-matrix penalization of Hamiltonian dynamics with application to exact sampling of stiff systems[END_REF][START_REF] Lelièvre | Free energy computations[END_REF], and references therein, for discussions on the use of constraints and penalizations in molecular dynamics. For overdamped Langevin dynamics (1.1), choosing a function r f ¡∇ r V with the potential r V V ¡ σ gives penalized Langevin dynamics [START_REF] Ciccotti | Projection of diffusions on submanifolds: application to mean force computation[END_REF] of the form

dY ε ptq f pY ε ptqqdt σdW ptq σ 2 4 ∇ lnpdetpGqqpY ε ptqqdt ¡ 1 ε pgζqpY ε ptqqdt, (1.3)
where we fix Y ε p0q X 0 , the parameter ε ¡ 0 is fixed with arbitrary size, f ¡∇V , g ∇ζ : R d Ñ R d¢q , and G g T g : R d Ñ R q¢q is the Gram matrix. It was shown in [16, Appx. C] that the solution Y ε of (1.3) converges strongly to the solution X 0 of the constrained dynamics (1.2) if X 0 M. The additional term σ 2 4 ∇ lnpdetpGqq is a correction term (called the Fixman correction) that is needed to obtain the convergence to the constrained dynamics (1.2) (see [START_REF] Lelièvre | Free energy computations[END_REF]Sect. 3.2.3.4] and references therein). Thus, for ε small, the trajectory of the solution of (1.3) lies in the vicinity of the manifold M. This penalization can also appear naturally when simulating Langevin dynamics with a stiff potential (see, for instance, [START_REF] Vilmart | Postprocessed integrators for the high order integration of ergodic SDEs[END_REF]Sect. 5.1]). One is then interested in numerical schemes that are robust with respect to the parameter ε and that lie on the manifold M in the limit ε Ñ 0. In this paper, we study the following similar penalized dynamics in R d to simulate trajectories in a vicinity of the manifold M:

dX ε ptq f pX ε ptqqdt σdW ptq σ 2 4 ∇ lnpdetpGqqpX ε ptqqdt ¡ 1
ε pgG ¡1 ζqpX ε ptqqdt, (1.4) where X ε p0q X 0 . It is a simpler version of (1.3) that also evolves in a vicinity of the manifold M in the limit ε Ñ 0. One result of this paper is the strong convergence of the solution X ε of (1.4) to the solution X 0 of (1.2) if X 0 M. We mention that in the deterministic setting, that is, when σ 0, equation (1.4) is a singular perturbation problem, and it converges to a differential algebraic equation (DAE) of index two in the limit ε Ñ 0 (see [22, Chaps. VI-VII]). We propose in this article a method that is robust with respect to the parameter ε for solving equations of the form (1.4), and we leave the creation of robust integrators for solving (1.3) for future work for the sake of clarity.

There are different ways to approximate the solution of the dynamics (1.4). A strong approximation focuses on approximating the realization of a single trajectory of (1.4) for a given realization of the Wiener process W . A weak approximation approximates the average of functionals of the solution at a fixed time T , that is, quantities of the form ErϕpX ε pTqqs for ϕ a smooth test function. In addition, under growth and smoothness assumptions on the vector fields in (1.4) (see, for instance, [START_REF] Hasminskii | Stochastic stability of differential equations[END_REF]), the dynamics (1.4) naturally satisfy an ergodicity property; that is, there exists a unique invariant measure dµ ε challenge when the dimension d is high, which is the case in the context of molecular dynamics where the dimension is proportional to the number of particles, as a standard quadrature formula becomes prohibitively expensive in high dimension. We emphasize that the invariant measure µ ε V becomes singular with respect to the Lebesgue measure on R d in the limit ε Ñ 0, and tends weakly as ε Ñ 0 to dµ 0 V , a measure that is absolutely continuous to dσ M , the canonical measure on M induced by the Euclidean metric of R d . In this paper, we propose weak convergence results for a new uniformly accurate integrator for solving (1.4), and numerical experiments in the weak context and for the invariant measure, as we recall that a scheme of weak order r automatically has order p ¥ r for the invariant measure (see, for instance, [START_REF] Mattingly | Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise[END_REF]).

The most used discretization for solving (1.4) is the explicit Euler integrator in R d (see [START_REF] Ciccotti | Blue moon sampling, vectorial reaction coordinates, and unbiased constrained dynamics[END_REF][START_REF] Lelièvre | Analyse de certains schémas de discrétisation pour des équations différentielles stochastiques contraintes[END_REF][START_REF] Lelièvre | Free energy computations[END_REF][START_REF] Lelièvre | Langevin dynamics with constraints and computation of free energy differences[END_REF], for instance),

X n 1 X n c hσξ hf pX n q h σ 2 4 ∇ lnpdetpGqqpX n q ¡ h ε pgG ¡1 ζqpX n q. (1.5)
This integrator has weak order one of accuracy, but it faces some severe stepsize restriction due to its instability, typically of the form h 3 ε, in order to be accurate in the regime ε Ñ 0. Since the solution X ε ptq of (1.4) converges to the solution X 0 ptq of (1.2) when ε Ñ 0, one can use integrators for the limit equation (1.2) and apply them to solve the original problem (1.4) when ε is close to zero. Indeed, one can prove that the solution X ε ptq of (1.4) stays at distance Op c εq of the constrained solution X 0 ptq of (1.2) (see Theorem 2.3). Thus, if the timestep of the integrator is small enough and satisfies ε 3 h, then this integrator is consistent for solving (1.4). The alternative for the discretization of (1.2) on the manifold of the explicit Euler scheme (1.5) is the constrained Euler scheme

X 0 n 1 X 0 n c hσξ hf pX 0 n q gpX 0 n qλ 0 n 1 , ζpX 0 n 1 q 0, (1.6)
where λ 0 n 1 R q acts as a Lagrange multiplier, and is entirely determined by the constraint ζpX 0 n 1 q 0. This integrator has weak order one for solving (1.2) (see [START_REF] Lelièvre | Free energy computations[END_REF]Sect. 3.2.4])

and it lies on the manifold M. It is a consistent approximation of (1.4) if ε is close to zero and ε 3 h. This integrator is, however, not appropriate for solving (1.4) if the size of ε is of the order of one, since the exact solution does not evolve in a neighborhood of the manifold in this regime. We mention a few other techniques to integrate numerically (1.4) or (1.2). In [START_REF] Laurent | Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs[END_REF][START_REF] Laurent | Order conditions for sampling the invariant measure of ergodic stochastic differential equations on manifolds[END_REF], high order Runge-Kutta methods are proposed for sampling the invariant measure in R d and on manifolds. The paper [START_REF] Lelièvre | Hybrid Monte Carlo methods for sampling probability measures on submanifolds[END_REF] presents a constrained integrator based on the RATTLE scheme (see [START_REF] Ryckaert | Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes[END_REF][START_REF] Andersen | Rattle: A "velocity" version of the Shake algorithm for molecular dynamics calculations[END_REF][START_REF] Hairer | Solving ordinary differential equations[END_REF]) in the context of the underdamped Langevin dynamics. Some of the previously cited discretizations can be combined with Metropolis-Hastings rejection procedures [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF][START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF], as done, for instance, in [START_REF] Girolami | Riemann manifold Langevin and Hamiltonian Monte Carlo methods[END_REF][START_REF] Brubaker | A family of MCMC methods on implicitly defined manifolds[END_REF][START_REF] Lelièvre | Langevin dynamics with constraints and computation of free energy differences[END_REF][START_REF] Zappa | Monte Carlo on manifolds: sampling densities and integrating functions[END_REF][START_REF] Lelièvre | Hybrid Monte Carlo methods for sampling probability measures on submanifolds[END_REF]. As for the Euler integrators (1.5)-(1.6), all of the previously mentioned methods are consistent provided that the parameter ε and the timestep satisfy h 3 ε in the context of methods in R d , and satisfy ε 3 h in the context of methods on the manifold M. When applied in the regime where ε and h share the same order of magnitude, the accuracy of these methods quickly deteriorates, and they may face stability issues.

In past decades, different solutions were proposed for treating the loss of accuracy in the intermediate regime h ε in the context of multiscale problems with the help of uniformly accurate (UA) methods. These methods are capable of solving dynamics indexed by a possibly stiff parameter ε with an accuracy and a cost both independent of ε. A uniformly accurate method is automatically asymptotic-preserving (AP), that is, it converges in the two regimes ε Ñ 0 and ε 1, but the converse is not true in general. We refer the reader to the review [START_REF] Jin | Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review[END_REF] and references therein for examples of AP integrators for solving multiscale problems. We mention in particular the paper [START_REF] Plecháč | Implicit mass-matrix penalization of Hamiltonian dynamics with application to exact sampling of stiff systems[END_REF], which proposes a penalized Hamiltonian dynamics, and AP discretizations for solving it, and the paper [START_REF] Bréhier | On asymptotic preserving schemes for a class of stochastic differential equations in averaging and diffusion approximation regimes[END_REF] that gives an AP scheme for the approximation of a class of multiscale SDEs. In [START_REF] Cohen | Convergence analysis of trigonometric methods for stiff secondorder stochastic differential equations[END_REF][START_REF] Vilmart | Weak second order multirevolution composition methods for highly oscillatory stochastic differential equations with additive or multiplicative noise[END_REF][START_REF] Laurent | Multirevolution integrators for differential equations with fast stochastic oscillations[END_REF], trigonometric and multirevolution integrators are considered for solving highly oscillatory SDEs (see also the deterministic works [START_REF] Melendo | A new approach to the construction of multirevolution methods and their implementation[END_REF][START_REF] Calvo | Approximate compositions of a near identity map by multi-revolution Runge-Kutta methods[END_REF][START_REF] Calvo | On explicit multi-revolution Runge-Kutta schemes[END_REF][START_REF] Chartier | Multi-revolution composition methods for highly oscillatory differential equations[END_REF]). We mention the recent papers [START_REF] Chartier | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF][START_REF] Chartier | A new class of uniformly accurate numerical schemes for highly oscillatory evolution equations[END_REF][START_REF] Chartier | Uniformly accurate numerical schemes for a class of dissipative systems[END_REF][START_REF] Almuslimani | Uniformly accurate schemes for oscillatory stochastic differential equations[END_REF] (see also the references therein) that introduce uniformly accurate methods for solving a variety of multiscale problems. There is a rich literature on AP and UA methods, but, to the best of our knowledge, the problem we study here and the techniques we consider are new. We propose in this paper a new consistent integrator with uniform accuracy and uniform cost for solving penalized Langevin dynamics, that is, a method for solving (1.4) whose accuracy and cost do not depend on the parameter ε.

The article is organized as follows. Section 2 is devoted to the presentation of the new integrator and to the main convergence results. In Section 3, we build a weak asymptotic expansion of the solution of (1.4) that is uniform in ε, and we use it for proving the uniform accuracy of our integrator. We compare in Section 4 the new integrator with the explicit Euler scheme in R d (1.5) and the constrained Euler scheme (1.6) on M in numerical experiments on a torus and on the orthogonal group to confirm its order of convergence in the weak context and for sampling the invariant measure. Finally, we present some possible extensions and future work in Section 5.

Uniformly accurate integrator for penalized Langevin dynamics

In this section, we present the new uniformly accurate integrator and the main convergence results of this paper. The proofs are postponed to Section 3. Let us first lay down a few notations and assumptions. We assume in the rest of the article that M tx R d , ζpxq 0u is a compact and smooth manifold of codimension q ¥ 1 embedded in R d , where the constraints are given by the smooth map ζ : R d Ñ R q . We write g ∇ζ : R d Ñ R d¢q and we assume that the Gram matrix Gpxq g T pxqgpxq R q¢q is invertible for all x in M. With these notations, the projection Π M on the tangent bundle is given by Π M pxq I d ¡G ¡1 pxqgpxqg T pxq.

f m φ fx i 1 . . . fx im pxq a 1 i 1 . . . a m im .
For the sake of clarity, we will often drop the coefficient x, and if m 1, we also use the notation φ I for the Jacobian matrix of φ. Moreover, for pe i q the canonical basis of R d , we write ∆φpxq d i1 φ P pxqpe i , e i q and pdiv φpxqq j d i1 φ I ij pxqpe i q.

In the rest of the paper, we make the following assumption in the spirit of the regularity assumptions made in [16, Under Assumption 2.1, the problems (1.2) and (1.4) are well posed, and we obtain the strong convergence of the penalized dynamics (1.4) to the constrained dynamics (1.2). Theorem 2.3. Under Assumption 2.1, the solution X ε of the penalized dynamics (1.4) converges strongly to X 0 , the solution of the constrained dynamics (1.2); that is, for all t ¤ T , there exists a constant C ¡ 0 such that, for all ε ¡ 0,

sup t¤T E § § X ε ptq ¡ X 0 ptq § § 2 % ¤ Cε.
Moreover, ζpX ε ptqq satisfies

sup t¤T E |ζpX ε ptqq| 2 % ¤ Cε.
This result was first introduced in [16, Appx. C] for slightly different penalized dynamics. The proof is almost identical, but we present it in Appx. A for the sake of completeness. In the deterministic setting with σ 0, the convergence to the manifold is of order 1 in ε instead of order 1{2. We introduce the following additional assumption. Assumption (2.4) is a stronger version of the inequality (2.1) and is in the spirit of the concept of admissible Lagrange multipliers [START_REF] Lelièvre | Hybrid Monte Carlo methods for sampling probability measures on submanifolds[END_REF]. It is always satisfied for x in a neighborhood of the manifold M and y in a ball centered on zero if we assume that the Gram matrix Gpxq G 0 pxq is invertible on M. We emphasize that we do not need this assumption for proving the uniform accuracy property of the new method, but we use it for obtaining uniform estimates in the regime ε Ñ 0 and on the numerical implementation of the uniformly accurate method.

We introduce the new integrator for approximating the penalized dynamics (1.4) with cost and accuracy independent of the parameter ε, and a cost comparable to that of the constrained Euler scheme (1.6) in terms of the number of evaluations of the functions f , ζ, g, and g I .

New Method (Uniform discretization of penalized overdamped Langevin dynamics)

X ε 0 X 0 M for n ¥ 0 do X ε n 1 X ε n c hσξ n hf pX ε n q p1 ¡ e ¡h{ε q 2 2 pg I pgG ¡1 ζqG ¡1 ζqpX ε n q σ 2 ε 8 p1 ¡ e ¡2h{ε q∇ lnpdetpGqqpX ε n q gpX ε n qλ ε n 1 , (2.2) ζpX ε n 1 q e ¡h{ε ζpX ε n q σ ε 2 p1 ¡ e ¡2h{ε qg T pX ε n qξ n εp1 ¡ e ¡h{ε q g T f σ 2 4 g T ∇ lnpdetpGqq σ 2 2 divpgq ¨pX ε n q σ 2 ¡ εp1 ¡ e ¡h{ε q ¡ εh 2 p1 ¡ e ¡2h{ε q © ¢ ¡ d i1 pg I pgG ¡1 g T e i qq T gG ¡1 g T e i ¡ d i1 pg I pe i qq T gG ¡1 g T e i © pX ε n q.
end for

The new method works in a way similar to the constrained Euler integrator (1.6). Knowing the approximation X ε n of X ε pnhq, we project a modified Euler step on a modified manifold defined by the constraint given in (2.2) (in place of ζpX ε n 1 q 0 for the constrained Euler integrator (1.6)). We project the modified step in the direction gpX ε n q with the help of a Lagrange multiplier λ ε n 1 . For the implementation of the method, one can use, for instance, a fixed point iteration or a Newton method at each step to find the solution pX ε n 1 , λ ε n 1 q of the implicit system of equations (2.2).

In order for the discretization (2.2) to be well-defined, we use bounded random variables. The ξ n are independent and bounded discrete random vectors that have the same moments as standard Gaussian random vectors up to order four, in the spirit of [START_REF] Milstein | Stochastic numerics for mathematical physics[END_REF]Chap. 2], that is, for instance, that their components satisfy Ppξ i 0q 2 3 and Ppξ i ¨c3q 1 6 , i 1, . . . , d.

(2.3) Note that using truncated Gaussian random variables would also work.

Remark 2.5. The new integrator is related to the popular idea of backward error analysis and modified equations for SDEs (see, for instance, [START_REF] Zygalakis | On the existence and the applications of modified equations for stochastic differential equations[END_REF][START_REF] Abdulle | High weak order methods for stochastic differential equations based on modified equations[END_REF][START_REF] Debussche | Weak backward error analysis for SDEs[END_REF][START_REF] Kopec | Weak backward error analysis for Langevin process[END_REF][START_REF] Kopec | Weak backward error analysis for overdamped Langevin processes[END_REF]). The idea is to define a projection method (see [START_REF] Hairer | Geometric numerical integration[END_REF]Sect. IV.4]) with a modified constraint in place of ζpX n q 0. Instead of evaluating the stiff term h ε gG ¡1 ζ as in the Euler scheme (1.5), we project a modified step of the explicit Euler scheme in R d on a manifold that is close to M when ε 3 1 and whose constraint is given by a truncation of a uniform expansion of ζpX ε q. When ε Ñ V, the expression of the constraint ζpX ε n 1 q in (2.2) tends to a truncated Taylor expansion in h around X ε n , while for ε Ñ 0, ζpX ε n 1 q tends to zero, which enforces that the integrator lies on M. These intuitions will be made rigorous in Section 3.

Remark 2.6. In the context of a manifold M of codimension q 1, the Gram matrix Gpxq and gpxq T e i g i pxq are real numbers, so that ∇ lnpdetpGqq 2G ¡1 g I pgq and

d i1 pg I pgG ¡1 g T e i qq T gG ¡1 g T e i G ¡1 pg I pgqq T g d i1 pg I pe i qq T gG ¡1 g T e i .
The discretization (2.2) thus reduces to

X ε n 1 X ε n c hσξ n hf pX ε n q p1 ¡ e ¡h{ε q 2 2 pζ 2 G ¡2 g I pgqqpX ε n q σ 2 ε 4 p1 ¡ e ¡2h{ε qpG ¡1 g I pgqqpX ε n q gpX ε n qλ ε n 1 , (2.4
)

ζpX ε n 1 q e ¡h{ε ζpX ε n q σ ε 2 p1 ¡ e ¡2h{ε qg T pX ε n qξ n εp1 ¡ e ¡h{ε qpg T f σ 2 2 G ¡1 g T g I pgq σ 2 2 divpgqqpX ε n q.
We present in the rest of the section the uniform accuracy property of the discretization (2.2), and we show that the integrator converges to the constrained Euler scheme (1.6) when ε Ñ 0. The different convergence results are summarized by the following commutative diagram, where T N h is fixed. Note that, as we present a convergence result in h that is uniform in ε, the two arrows for the convergence in h rely on the same Theorem 2.7.

integrator X ε N in R d integrator X 0 N on M solution X ε pTq of (1.4) solution X 0 pTq of (1.2) hÑ0 pThm. 2.7q εÑ0 pThm. 2.9q hÑ0 pThm. 2.7q εÑ0 pThm. 2.3q
We now state the main result of this work, that is, the uniform accuracy of the discretization (2.2) for approximating the solution of the penalized Langevin dynamics (1.4).

Theorem 2.7. Under Assumption 2.1, the integrator pX ε n q given by (2.2) is a consistent uniformly accurate approximation of the solution X ε ptq of the penalized Langevin dynamics (1.4); that is, for a given test function ϕ C 5 P , there exist h 0 ¡ 0, C ¡ 0 such that for all ε ¡ 0, h ¤ h 0 , the following estimate holds:

|ErϕpX ε n qs ¡ ErϕpX ε pnhqqs| ¤ C c h, n 0, 1, . . . , N, N h T.
Remark 2.8. Note that Theorem 2.7 states the uniform consistency, but not the uniform weak order one, as one could expect. The discretization (2.2) has weak order one if ε ε 0 is fixed or in the limit ε Ñ 0, but an order reduction occurs in the intermediate regime 0 ε ε 0 , and the integrator only has weak order 1{2 with respect to h in general. For the sake of simplicity, we leave the creation of uniformly accurate integrators of higher weak order for future works.

We present the proof of Theorem 2.7 in Section 3. It relies on a weak expansion in h of the solution of (1.4) that is uniform in ε. One could directly use this uniform expansion as an explicit numerical integrator for solving (1.4). It would also yield a uniformly accurate scheme and would not require one to solve a fixed point problem. However, in the limit ε Ñ 0, this integrator would almost surely not stay on the manifold. The crucial geometric property that the integrator lies on the manifold when ε Ñ 0 is satisfied for the new method, as stated in the following result. Theorem 2.9. Under Assumption 2.1, the integrator pX ε n q in (2.2) converges to the Euler scheme on the manifold (1.6) when ε Ñ 0; that is, for h and N fixed such that T N h, there exists a constant C h ¡ 0 that depends on h but not on ε such that

|ζpX ε n q| ¤ C h c ε, n 0, 1, . . . , N. In addition, if Assumption 2.4 is satisfied, then, for h 0 small enough, for h ¤ h 0 fixed, there exists a constant C h ¡ 0 such that § § X ε n ¡ X 0 n § § ¤ C h c ε, n 0, 1, . . . , N, N h T.
(

Remark 2.10. In the deterministic context (i.e. when σ 0), the uniform accuracy of the discretization (2.2) still holds, and the speed of convergence to the manifold M of both the exact solution and the integrator are in Opεq. To the best of our knowledge, the integrator given by (2.2) is the first integrator with the uniform accuracy property for solving the singular perturbation problem (1.4) with σ 0. However, similar expansions that are uniform with respect to ε are presented in [22, Chaps. VI-VII] and references therein.

Remark 2.11. Another widely used scheme on the manifold is the Euler scheme with implicit projection direction,

X 0 n 1 X 0 n hf pX 0 n q c hσξ n gpX 0 n 1 qλ 0 n 1 , ζpX 0 n 1 q 0, (2.7)
where the Lagrange multiplier λ 0 n 1 is determined by the constraint ζpX 0 n 1 q 0. The uniformly accurate discretization given in (2.2) can be modified so that it converges to the integrator (2.7) when ε Ñ 0. It suffices to replace the first line of (2.2)

X ε n 1 X ε n c hσξ n hf pX ε n q ¡ p1 ¡ e ¡h{ε q 2 2 pg I pgG ¡1 ζqG ¡1 ζqpX ε n q σ 2 4 ¡ 2εhp1 ¡ e ¡2h{ε q ¡ ε 2 p1 ¡ e ¡2h{ε q © ∇ lnpdetpGqqpX ε n q gpX ε n 1 qλ ε n 1 ,
and to keep the same expansion for the constraint ζpX ε n 1 q. The methodology for the uniform expansion of the integrator that we present in Section 3.2 extends to this context, so that the convergence results persist. Similarly, one could change the direction of projection gpX ε n q into gpY ε n q, where Y ε n is any consistent one-step approximation of X ε n , in the spirit of the class of projected Runge-Kutta methods presented in [START_REF] Laurent | Order conditions for sampling the invariant measure of ergodic stochastic differential equations on manifolds[END_REF]. Finding a class of uniformly accurate discretizations that converge to a more general class of Runge-Kutta methods on the manifold M is matter for future works.

The uniform discretization (2.2) is implicit and requires one to solve a fixed point problem at each step with, for instance, a fixed point iteration or a Newton method. The following result, in the spirit of [START_REF] Hairer | Solving ordinary differential equations[END_REF]Chap. VII] for deterministic DAEs and [START_REF] Laurent | Order conditions for sampling the invariant measure of ergodic stochastic differential equations on manifolds[END_REF]Lemma 3.3] for the constrained dynamics (1.2), confirms that the associated implicit system is not stiff, that is, that its complexity does not depend on the stiff parameter ε. Theorem 2.12. Under Assumption 2.1, each step of the integrator pX ε n q given by (2.2) can be rewritten as a solution of a fixed point problem of the form

X ε n 1 F ε h pX ε n 1 q, where F ε h : R d Ñ R d depends on X ε n , ξ n , h, and ε. Moreover, if Assumption 2.4 is satisfied, then there exists h 0 ¡ 0 independent of ε such that for all h ¤ h 0 , F ε h is a uniform contraction,
that is, there exists a positive constant L 1 independent of h and ε such that, for all y 1 , y 2

R d , |F ε h py 2 q ¡ F ε h py 1 q| ¤ L |y 2 ¡ y 1 | .

Weak convergence analysis

In this section, we present the uniform weak expansion and the stability properties of the solution X ε ptq to the penalized dynamic (1.4) and of the uniform integrator pX ε n q given in (2.2).

We then use these results to prove Theorem 2.7, Theorem 2.9 and Theorem 2.12. Let us begin the analysis with a few technical lemmas and notations.

Lemma 3.1. Let pX ε n q be given by (2.2), then, under Assumption 2.1, there exists a constant C 0 ¡ 0 independent of X 0 M, ε, and h such that, for all n ¥ 0, p1 ¡ e ¡h{ε q |ζpX ε n q|

¤ C 0 c h. (3.1)
Proof. As f , ξ n , g, and g I are bounded, we obtain from the definition of the integrator given by (2.2) that

|ζpX ε n 1 q| ¤ e ¡h{ε |ζpX ε n q| C εp1 ¡ e ¡2h{ε
q, where we used that the function p1 ¡ e ¡x q{x is bounded for x ¡ 0. Thus, as ζpX 0 q 0,

|ζpX ε n q| ¤ e ¡nh{ε |ζpX 0 q| C εp1 ¡ e ¡2h{ε q n¡1 ķ0 e ¡kh{ε ¤ C εp1 ¡ e ¡2h{ε q 1 ¡ e ¡h{ε ¤ C 0 c h 1 ¡ e ¡h{ε ,
where the constant C 0 does not depend on X 0 , ε, n, and h. This yields the estimate (3.1).

The estimate (3.1) is a direct consequence of our choice of using bounded random variables in (2.2). We shall use this estimate extensively in the rest of this section. Thus, we denote by M ε h the set of vectors x R d that satisfy the estimate p1 ¡ e ¡h{ε q |ζpxq| ¤ C 0 c h, (

where C 0 is the constant given in Lemma 3.1. The set M ε h is a closed subset of R d that contains M. The numerical scheme given by (2.2) takes values in M ε h . We mention that the convergence results are still valid if the initial condition X 0 of (2.2) is chosen in M ε h instead of M.

As we aim at writing uniform expansions, we introduce the convenient notation R ε h pxq for any remainder that satisfies at least |ErR ε h pxqs| ¤ Ch 3{2 , where C is independent of ε, h, and x.

The following result serves as a technical tool for simplifying the calculations in the uniform expansions of Subsection 3.1 and Subsection 3.2. It is proved with elementary computations. Lemma 3.2. The Fixman correction can be rewritten in the following way:

σ 2 4 ∇ lnpdetpGqq σ 2 2 d i1 g I pe i qG ¡1 g T e i σ 2 2 d i1 g I pgG ¡1 g T e i qG ¡1 g T e i ,
where pe i q is the canonical basis of R d . Proof of Lemma 3.2. Using that G g T g is a symmetric matrix, that g ∇ζ is a gradient (which implies x T g I pyq y T g I pxq) and the standard properties of the trace operator Tr, we deduce that f j plnpdetpGqqq TrpG ¡1 f j Gq 2 TrpG ¡1 pg I pe j qq T gq 2 Trpg I pe j qG ¡1 g T q

2 d i1 e T i g I pe j qG ¡1 g T e i 2 d i1 e T j g I pe i qG ¡1 g T e i ,
that is, ∇ lnpdetpGqq 2 °d i1 g I pe i qG ¡1 g T e i . For the second equality, we have

d i1 e T j g I pgG ¡1 g T e i qG ¡1 g T e i d i1
e T i gG ¡1 g T g I pe j qG ¡1 g T e i TrpgG ¡1 g T g I pe j qG ¡1 g T q Trpg I pe j qG ¡1 g T q d i1 e T i g I pe j qG ¡1 g T e i d i1 e T j g I pe i qG ¡1 g T e i .

Hence we get the result.

Uniform expansion of the exact solution

We consider the exact solution X ε ptq of the penalized Langevin dynamics (1.4), with the initial condition X 0 x, that we assume is deterministic for simplicity. Then, X ε satisfies the following expansion in h that is uniform with respect to ε. Proposition 3.3. Under Assumption 2.1, there exists h 0 ¡ 0 such that for all h ¤ h 0 , if X ε is the solution of the penalized Langevin dynamics (1.4) starting at x M ε h , then, for all ϕ C 3 P , the following estimate holds: § § §ErϕpX

ε phqqs ¡ Erϕpx c hA ε h pxq hB ε h pxqqs § § § ¤ Cp1 |x| K qh 3{2 , (3.3)
where C is independent of h and ε and where the functions A ε h and B ε h are given by

A ε h σξ e ¡h{ε ¡ 1 c h gG ¡1 ζ σ ¡ ε 2h p1 ¡ e ¡2h{ε q ¡ 1 © gG ¡1 g T ξ, B ε h f ¡ ε h p1 ¡ e ¡h{ε q ¡ 1 © gG ¡1 g T f σ 2 2 ¡ ε h p1 ¡ e ¡h{ε q ¡ 1 © gG ¡1 divpgq σ 2 ε 8h p1 ¡ e ¡2h{ε q∇ lnpdetpGqq σ 2 ε 8h p1 ¡ e ¡h{ε q 2 gG ¡1 g T ∇ lnpdetpGqq 1 2h pe ¡h{ε ¡ 1q 2 ¡ g I pgG ¡1 ζqG ¡1 ζ ¡ gG ¡1 g T g I pgG ¡1 ζqG ¡1 ζ ¡ gG ¡1 pg I pgG ¡1 ζqq T gG ¡1 ζ © σ 2 ¡ 1 ε h pe ¡h{ε ¡ 1q © d i1 gG ¡1 pg I pe i qq T gG ¡1 g T e i σ 2 4 ¡ ε h pe ¡2h{ε ¡ 4e ¡h{ε 3q ¡ 2 © d i1 gG ¡1 pg I pgG ¡1 g T e i qq T gG ¡1 g T e i ,
with ξ a discrete bounded random vector that satisfies (2.3), and where the functions A ε h and B ε h are bounded uniformly in ε and h on M ε h .

Note that for q 1, we obtain the first step of the uniformly accurate discretization (2.4) by gathering all the terms of the form gM with M R q of the weak approximation given in Proposition 3.3 in a Lagrange multiplier gλ ε 1 R q and by adding the truncated expansion of ζpX ε phqq as a constraint.

The proof of Proposition 3.3 relies on the change of coordinate ψ given in Assumption 2.1. Instead of discretizing directly the penalized dynamics (1.4), we first apply the change of coordinate ψ and we derive an expansion in time of X ε ptq that is uniform in the parameter ε.

The following result is used for proving Proposition 3.3.

Lemma 3.4. With the same notations and assumptions as in Proposition 3.3, the following

estimates hold for all h ¤ h 0 and all x M ε h :

Er|X ε phq ¡ x| 2 s 1{2 ¤ C c h, ( 3.4
)

Er|X ε phq ¡ px c h p A ε h pxqq| 2 s 1{2 ¤ Ch, (3.5) § § §ErψpX ε phqqs ¡ Erψpx c h p A ε h pxq hB ε h pxqqs § § § ¤ Ch 3{2 , (3.6)
where C is independent of h and ε, B ε h is defined in Proposition 3.3, and p A ε h is given by

p A ε h σ W phq c h e ¡h{ε ¡ 1 c h gG ¡1 ζ σ c h gG ¡1 g T » h 0 pe ps¡hq{ε ¡ 1qdW psq.
Proof of Proposition 3.3. The uniform bounds on A ε h pX ε ptqq and B ε h pX ε ptqq are obtained straightforwardly by using Assumption 2.1 and the fact that x M ε h . We prove the local weak order one of the approximation Y ε phq x c

h p A ε h pxq hB ε h pxq given in Lemma 3.4. Let ϕ C 3 P and r ϕ ϕ ¥ ψ ¡1 ; then a Taylor expansion around x yields § § §ErϕpX ε phqq ¡ ϕpY ε phqqs § § § § § §Er r ϕ ¥ ψpX ε phqqs ¡ Er r ϕ ¥ ψpY ε phqqs § § § ¤ § § §Er r ϕ I pψpxqqpψpX ε phqq ¡ ψpY ε phqqqs § § § 1 2 § § §Er r ϕ P pψpxqqpψpX ε phqq ¡ ψpxq, ψpX ε phqq ¡ ψpxqq ¡ r ϕ P pψpxqqpψpY ε phqq ¡ ψpxq, ψpY ε phqq ¡ ψpxqqs § § § Cp1 |x| K qh 3{2 ¤ § § § r ϕ I pψpxqq § § § § § §ErψpX ε phqq ¡ ψpY ε phqqs § § § 1 2 § § § r ϕ P pψpxqq § § §Er|ψpX ε phqq ψpY ε phqq ¡ 2ψpxq| 2 s 1{2 ¤ Er|ψpX ε phqq ¡ ψpY ε phqq| 2 s 1{2 Cp1 |x| K qh 3{2
where we used (3.4), Assumption 2.1, and the bilinearity of r ϕ P pψpxqq. With Lemma 3.4 and the regularity properties of r ϕ and ψ, we get § § §ErϕpX

ε phqqs ¡ Erϕpx c h p A ε h pxq hB ε h pxqqs § § § ¤ Cp1 |x| K qh 3{2 . ( 3
A ε h,j pxqq σ 2 δ ij 2σ 2 h pgG ¡1 g T q ij » h 0 pe ps¡hq{ε ¡ 1qds σ 2 h d ķ1 pgG ¡1 g T q ik pgG ¡1 g T q jk » h 0 pe ps¡hq{ε ¡ 1q 2 ds σ 2 δ ij σ 2 pgG ¡1 g T q ij ¡ ε 2h p1 ¡ e ¡2h{ε q ¡ 1 © ,
where we used that g T g G and the Itô isometry. On the other hand, a similar calculation yields

CovpA ε h,i pxq, A ε h,j pxqq σ 2 δ ij 2σ 2 ¡ ε 2h p1 ¡ e ¡2h{ε q ¡ 1 © pgG ¡1 g T q ij σ 2 ¡ ε 2h p1 ¡ e ¡2h{ε q ¡ 1 © 2 d ķ1 pgG ¡1 g T q ik pgG ¡1 g T q jk σ 2 δ ij σ 2 pgG ¡1 g T q ij ¡ ε 2h p1 ¡ e ¡2h{ε q ¡ 1 © .
Replacing p A ε h pxq by A ε h pxq in the weak expansion (3.7) gives the estimate (3.3).

The main ingredient of the proof of Lemma 3.4 is the decomposition of the terms of the expansion in a part that stays on the tangent space and a part of the form gM with M R q that is orthogonal to the tangent space.

Proof of Lemma 3.4. As ψ ¡1 is Lipschitz, we have

Er|X ε phq ¡ x| 2 s 1{2 ¤ CEr|ψpX ε phqq ¡ ψpxq| 2 s 1{2 ¤ CEr|φpX ε phqq ¡ φpxq| 2 s 1{2 CEr|ζpX ε phqq ¡ ζpxq| 2 s 1{2 .
On the one hand, applying the Itô formula to φpX ε q yields φpX ε phqq φpxq σ

» h 0 φ I pX ε psqqdWpsq (3.8) » h 0 rφ I f σ 2 4 φ I ∇ lnpdetpGqq σ 2 2 ∆φspX ε psqqds,
where the term in ε vanishes as φ I g 0 (see Assumption 2.1). Assumption 2.1 allows us to write the uniform strong expansion

Er|φpX ε phqq ¡ φpxq| 2 s 1{2 ¤ C c h.
On the other hand, for ζpX ε q, we have

dζpX ε q σg T pX ε qdW g T f σ 2 4 g T ∇ lnpdetpGqq σ 2 2 divpgq ¡ 1 ε ζ % pX ε qdt.
With the variation of constants formula, it rewrites into

ζpX ε phqq e ¡h{ε ζpxq σ » h 0 e ps¡hq{ε g T pX ε psqqdWpsq (3.9) » h 0 e ps¡hq{ε g T f σ 2 4 g T ∇ lnpdetpGqq σ 2 2 divpgq
% pX ε psqqds.

As the integrands in (3.9) are bounded (using Assumption 2.1), we get

Er|ζpX ε phqq ¡ ζpxq| 2 s 1{2 ¤ Cppe ¡h{ε ¡ 1q 2 |ζpxq| 2 hq 1{2 ¤ C c h,
where we used that x M ε h . We thus get the desired estimate (3.4). The estimate (3.5) is obtained with the same arguments by keeping track of the terms of size Op c hq in the expansions.

We now prove the weak estimate (3.6). We denote for simplicity Y ε phq x c

h p A ε h pxq hB ε h pxq.
Let us first look at the approximation of φpX ε phqq. On the one hand, applying the Itô formula to φpX ε ptqq gives

φpX ε phqq φpxq hφ I f pxq h σ 2 4 φ I ∇ lnpdetpGqqpxq h σ 2 2 ∆φpxq R ε h pxq,
where we used (3.4) and we put in R ε h pxq all the terms that are zero in average. On the other hand, an expansion in h of φpY ε phqq yields

φpY ε phqq φ c hφ I p A ε h h φ I B ε h 1 2 φ P p p A ε h , p A ε h q % R ε h φ hφ I f σ 2 ε 8 p1 ¡ e ¡2h{ε qφ I ∇ lnpdetpGqq 1 2 pe ¡h{ε ¡ 1q 2 φ I g I pgG ¡1 ζqG ¡1 ζ σ 2 2 φ P pWphq, W phqq 1 2 pe ¡h{ε ¡ 1q 2 φ P pgG ¡1 ζ, gG ¡1 ζq σ 2 2 » h 0 » h 0 pe ps¡hq{ε ¡ 1qpe pu¡hq{ε ¡ 1qφ P pgG ¡1 g T dW psq, gG ¡1 g T dW puqq σ 2 » h 0 pe ps¡hq{ε ¡ 1qφ P pWphq, gG ¡1 g T dW psqq R ε h ,
where we use that φ I g 0, we omit the dependence in x for conciseness, and we put in R ε h pxq all the terms that are zero in average. We now replace the random terms by their expectation,

φpY ε phqq φ hφ I f σ 2 ε 8 p1 ¡ e ¡2h{ε qφ I ∇ lnpdetpGqq h σ 2 2 ∆φ 1 2 pe ¡h{ε ¡ 1q 2 pφ I pg I pgG ¡1 ζqG ¡1 ζq φ P pgG ¡1 ζ, gG ¡1 ζqq σ 2 2 » h 0 pe ps¡hq{ε ¡ 1q 2 ds d i1 φ P pgG ¡1 g T e i , gG ¡1 g T e i q σ 2 » h 0 pe ps¡hq{ε ¡ 1qds d i1 φ P pe i , gG ¡1 g T e i q R ε h .
Letting M R q , we differentiate the equality φ I pgMq 0. We obtain that for any M R q and v R d we have φ I pg I pvqMq φ P pgM, vq 0. We deduce that φ I pg I pgG ¡1 ζqG ¡1 ζq φ P pgG ¡1

ε phqqs ¡ Erφpx c h p A ε h pxq hB ε h pxqqs § § § ¤ Ch 3{2 .
(3.10)

For the one-step approximation of ζpX ε phqq, the Itô formula and the variation of constants formula yield

ζpX ε phqq e ¡h{ε ζpxq εp1 ¡ e ¡h{ε qpg T f σ 2 4 g T ∇ lnpdetpGqq σ 2 2 divpgqqpxq R ε h pxq.
For ζpY ε phqq, using g T g G, we get with the same arguments as for φpY ε phqq that

ζpY ε phqq ζ c hg T p A ε h h g T B ε h 1 2 pg I p p A ε h qq T p A ε h % R ε h e ¡h{ε ζ εp1 ¡ e ¡h{ε qg T f σ 2 2 ¡ εp1 ¡ e ¡h{ε q ¡ h © divpgq σ 2 2 pg I pWphqqq T W phq σ 2 ε 4 p1 ¡ e ¡h{ε q 2 g T ∇ lnpdetpGqq σ 2 ¡ h εpe ¡h{ε ¡ 1q © d i1 pg I pe i qq T gG ¡1 g T e i σ 2 4 ¡ εpe ¡2h{ε ¡ 4e ¡h{ε 3q ¡ 2h © d i1 pg I pgG ¡1 g T e i qq T gG ¡1 g T e i σ 2 2 » h 0 » h 0 pe ps¡hq{ε ¡ 1qpe pu¡hq{ε ¡ 1qpg I pgG ¡1 g T dW puqqq T gG ¡1 g T dW psq σ 2 » h 0 » h 0 pe ps¡hq{ε ¡ 1qpg I pWphqqq T gG ¡1 g T dW psq R ε h .
We now replace the stochastic integrals by their expectations (putting the remainders in R ε h ), and we use Lemma 3.2 to simplify the expansion. It yields

ζpY ε phqq e ¡h{ε ζ εp1 ¡ e ¡h{ε qpg T f σ 2 4 g T ∇ lnpdetpGqq σ 2 2 divpgqq R ε h , which implies § § §ErζpX ε phqqs ¡ Erζpx c h p A ε h pxq hB ε h pxqqs § § § ¤ Ch 3{2 . (3.11)
Combining the inequalities (3.10) and (3.11) gives the desired weak estimate (3.6).

To end this subsection, we recall the growth properties of X ε that will be of use in Subsection 3.3. For this particular result, we add the dependency in the initial condition x of the exact solution of (1.4) with the notation X ε pt, xq. Lemma 3.5. Under Assumption 2.1, for ϕ C 5 P and t ¤ T fixed, the function r ϕ ε pxq ErϕpX ε pt, xqqs lies in C 3 P with constants independent of t and ε.

Proof. The standard theory (see, for instance, the textbook [START_REF] Milstein | Stochastic numerics for mathematical physics[END_REF]) gives r ϕ ε C 3 . With the regularity assumptions on ψ and its derivatives, it is sufficient to prove that ErϕpY ε pt, yqqs is in C 3 P , where Y ε pt, yq ψpX ε pt, ψ ¡1 pxqqq (replacing ϕ by ϕ ¥ ψ ¡1 ). We recall from the proof of Lemma 3.4 that φpX ε pt, xqq satisfies the integral formulation (3.8) and that ζpX ε pt, xqq satisfies (3.9). Putting together these equations, we deduce that Y ε pt, yq satisfies an equation of the form

Y ε pt, yq A ε ptqy » t 0 A ε pt ¡ sqF pY ε ps, yqqds » t 0 A ε pt ¡ sqGpY ε ps, yqqdW psq,
where F and G are in C 3 P and do not depend on ε, and

A ε ptq ¢ I d¡q 0 0 e ¡t{ε I q . As |A ε ptq| ¤ C, the process Y ε pt, yq satisfies Er|Y ε pt, yq| 2p s ¤ Cp1 |y| K q, with C and K independent of ε. Thus we have ErϕpY ε pt, yqqs ¤ Cp1 Er|Y ε pt, yq| K sq ¤ Cp1 Er|y| K sq.
For the derivatives, we recall from [START_REF] Gikhman | Stochastic differential equations[END_REF] that Z ε pt, yq f y Y ε pt, yq satisfies the equation

Z ε pt, yq A ε ptqI d » t 0 A ε pt ¡ sqF I pY ε ps, yqqZ ε pt, yqds » t 0 A ε pt ¡ sqG I pY ε ps, yqqZ ε pt, yqdW psq.
Applying the same arguments as for Y ε pt, yq yields that Z ε pt, yq has bounded moments of all order and that ErϕpZ ε pt, yqqs ¤ Cp1 Er|y| K sq. The same methodology extends to f 2 y Y ε pt, yq and f 3 y Y ε pt, yq.

Uniform expansion and bounded moments of the numerical solution

In this subsection, we show that the integrator given by (2.2) has bounded moments of all order, that it lies on the manifold M in the limit ε Ñ 0, and that it satisfies the same local weak uniform expansion as the exact solution of (1.4) (see Proposition 3.3). First, the integrator given by (2.2) satisfies the following bounded moments property.

Proposition 3.6. Under Assumption 2.1, pX ε n q has bounded moments of any order along time; i.e., for all timestep h ¤ h 0 small enough such that N h T is fixed, for all integer m ¥ 0, sup

n¤N Er|X ε n | 2m s ¤ C m , where the constant C ¡ 0 is independent of ε and h.
The integrator pX ε n q given in (2.2) also satisfies a uniform local expansion that is similar to its continuous counterpart presented in Proposition 3.3.

Proposition 3.7. Under Assumption 2.1, there exists h 0 ¡ 0 such that for all h ¤ h 0 , if pX ε n q is the numerical discretization given by (2.2) beginning at x M ε h (assumed deterministic for simplicity), then, for all test functions ϕ C 3 P , the following estimate holds:

§ § §ErϕpX ε 1 qs ¡ Erϕpx c hA ε h pxq hB ε h pxqqs § § § ¤ Cp1 |x| K qh 3{2 , (3.12)
where C is independent of h and ε, and A ε h and B ε h are the functions given in Proposition 3.3.

To prove Proposition 3.6 and Proposition 3.7, we rely on the following lemma, whose proof is postponed to the end of this subsection. We emphasize that an inequality of the form |ζpxq| ¤ C does not imply in general that x stays close to M. That is why we rely in Lemma 3.8 on an estimate of the Lagrange multipliers (using the inequality (2.1)). This estimate ensures that the method evolves in a neighborhood of the manifold.

Lemma 3.8. Under Assumption 2.1 and if x M ε

h , there exists h 0 ¡ 0 such that, for all timestep h ¤ h 0 , the one-step approximation X ε 1 and the Lagrange multiplier

λ ε 1 in the discretization (2.2) satisfy |X ε 1 ¡ x| ¤ C c h, λ ε 1 c hG ¡1 pxqλ ε 1,p1{2q hG ¡1 pxqλ ε 1,p1q R ε h pxq, where |λ ε 1,p1{2q | ¤ C, |λ ε 1,p1q | ¤ C, and |R ε h pxq| ¤ Ch 3{2 with C independent of ε, h and x.
For proving Proposition 3.6, we apply the change of variable ψ and we adapt the standard methodology presented in [START_REF] Milstein | Stochastic numerics for mathematical physics[END_REF]

, Lemma 1.1.6 & Lemma 2.2.2].
Proof of Proposition 3.6. We derive from (2.2) that § § §ErζpX

ε n 1 q ¡ e ¡h{ε ζpX ε n q|X ε n s § § § ¤ Ch, § § §ζpX ε n 1 q ¡ e ¡h{ε ζpX ε n q § § § ¤ C c h.
We prove that ζpX n q has bounded moments by induction on n. The binomial formula yields

Er § § ζpX ε n 1 q § § 2m s Er § § §e ¡h{ε ζpX ε n q ζpX ε n 1 q ¡ e ¡h{ε ζpX ε n q § § § 2m s ¤ e ¡2mh{ε Er|ζpX ε n q| 2m s CEr|ζpX ε n q| 2m¡1 § § §ErζpX ε n 1 q ¡ e ¡h{ε ζpX ε n q|X ε n s § § §s C 2m ķ2 Er|ζpX ε n q| 2m¡k § § §ζpX ε n 1 q ¡ e ¡h{ε ζpX ε n q § § § k s ¤ Er|ζpX ε n q| 2m s Cp1 Er|ζpX ε n q| 2m sqh. Following [40, Lemma 1.1.6],
as X 0 M is bounded, it implies that Er|ζpX n q| 2m s is bounded uniformly in n 0, . . . , N and ε.

Using Lemma 3.8 and the equality φ I g 0, a direct calculation gives

§ § ErφpX ε n 1 q ¡ φpX ε n q|X ε n s § § ¤ Ch, (3.13) § § φpX ε n 1 q ¡ φpX ε n q § § ¤ C c h. (3.14)
Following the same methodology as for Er|ζpX n q| 2m s, the estimates (3.13)-(3.14) imply that the quantity Er|φpX n q| 2m s is bounded uniformly in n 0, . . . , N and ε. Then, as ψ ¡1 is Lipschitz, we have

Er|X ε n | 2m s ¤ Cp1 Er|ψpX ε n q| 2m sq ¤ Cp1 Er|φpX ε n q| 2m s Er|ζpX ε n q| 2m sq ¤ C.
Hence we get the result.

We obtain the uniform expansion of the numerical solution by writing explicitly a uniform expansion of the Lagrange multiplier λ ε 1 , in the spirit of [START_REF] Lelièvre | Free energy computations[END_REF]Lemma 3.25].

Proof of Proposition 3.7. Using Lemma 3.8 and Assumption 2.1, we obtain

X ε 1 x c h σξ pgG ¡1 qpxqλ ε 1,p1{2q % R ε
h pxq, where the remainder satisfies |R ε h pxq| ¤ Ch. The constraint is then given by

ζpX ε 1 q ζpxq c h σg T pxqξ λ ε 1,p1{2q % R ε h pxq.
On the other hand, we get from the definition of the integrator (2.2) that

ζpX ε 1 q ζpxq c h e ¡h{ε ¡ 1 c h ζpxq σ ε 2h p1 ¡ e ¡2h{ε qg T pxqξ % R ε h pxq.
By identifying the two terms in c h in the expansions of ζpX ε 1 q, we deduce the value of λ ε 1,p1{2q , that is,

λ ε 1,p1{2q e ¡h{ε ¡ 1 c h ζ σ ¡ ε 2h p1 ¡ e ¡2h{ε q ¡ 1 © g T ξ.
The expression of λ ε 1,p1{2q and Lemma 3.8 give

X ε 1 x c hA ε h pxq h f pxq p1 ¡ e ¡h{ε q 2 2h pg I pgG ¡1 ζqG ¡1 ζqpxq σ 2 ε 8h p1 ¡ e ¡2h{ε q∇ lnpdetpGqqpxq pgG ¡1 qpxqλ ε 1,p1q % R ε h pxq, (3.15)
where |R ε h pxq| ¤ Ch 3{2 . We then compute the expansion of ζpX ε 1 q, and we compare it with the definition of the integrator (2.2) to obtain the expression of λ ε 1,p1q . Inserting this expression in (3.15) gives

X ε 1 x c hA ε h pxq hB ε h pxq R ε h pxq, where the remainder satisfies |ErR ε h pxqs| ¤ Ch 3{2 and |R ε h pxq| ¤ Ch. A Taylor expansion of ϕpX ε 1 q around x c hA ε h pxq hB ε h pxq yields the estimate (3.12).
The proof of Lemma 3.8 mainly relies on the estimates (2.1) and (3.2). We refer the reader to [START_REF] Laurent | Order conditions for sampling the invariant measure of ergodic stochastic differential equations on manifolds[END_REF]Lemma 3.3] and [22, Chap. VII] for similar proofs where explicit expressions of Lagrange multipliers are derived.

Proof of Lemma 3.8. Let x M ε h . For brevity, we rewrite the discretization (2.2) as

X ε 1 x c hY ε p1{2q pxq hY ε p1q pxq gpxqλ ε 1 , ζpX ε 1 q ζpxq c hζ ε p1{2q pxq hζ ε p1q pxq,
where the functions Y ε p1{2q , Y ε p1q , ζ ε p1{2q , and ζ ε p1q are given by

Y ε p1{2q σξ, Y ε p1q f p1 ¡ e ¡h{ε q 2 2h pg I pgG ¡1 ζqG ¡1 ζq σ 2 ε 8h p1 ¡ e ¡2h{ε q∇ lnpdetpGqq, ζ ε p1{2q e ¡h{ε ¡ 1 c h ζ σ ε 2h p1 ¡ e ¡2h{ε qg T ξ, ζ ε p1q ε h p1 ¡ e ¡h{ε q g T f σ 2 4 g T ∇ lnpdetpGqq σ 2 2 divpgq ¨ σ 2 ¡ ε h p1 ¡ e ¡h{ε q ¡ ε 2h p1 ¡ e ¡2h{ε q © d i1 ppg I pgG ¡1 g T e i qq T gG ¡1 g T e i q ¡ σ 2 ¡ ε h p1 ¡ e ¡h{ε q ¡ ε 2h p1 ¡ e ¡2h{ε q © d i1 ppg I pe i qq T gG ¡1 g T e i q.
Using Assumption 2.1 and the estimate (3.2), the following uniform estimates hold:

§ § §Y ε piq pxq § § § ¤ C, § § §ζ ε piq pxq § § § ¤ C, i t1{2, 1u.
The fundamental theorem of calculus yields

ζpX ε 1 q ¡ ζpxq » 1 0 g T px τ pX ε 1 ¡ xqqdτ pX ε 1 ¡ xq c hζ ε p1{2q pxq hζ ε p1q pxq. Substituting X ε 1 ¡ x then gives » 1 0 g T px τ pX ε 1 ¡ xqqdτ p c hY ε p1{2q pxq hY ε p1q pxq gpxqλ ε 1 q c hζ ε
p1{2q pxq hζ ε p1q pxq.

Using Assumption 2.1, we get the following explicit expression of λ ε 1 :

λ ε 1 c hG ¡1 X ε 1 ¡x pxq ¡ ζ ε p1{2q pxq ¡ » 1 0 g T px τ pX ε 1 ¡ xqqdτ Y ε p1{2q pxq © (3.16) hG ¡1 X ε 1 ¡x pxq ¡ ζ ε p1q pxq ¡ » 1 0 g T px τ pX ε 1 ¡ xqqdτ Y ε p1q pxq © .
Then, the growth assumption (2.1) on G ¡1 y pxq allows us to write

|λ ε 1 | ¤ C c hp1 |X ε 1 ¡ x|q and |X ε 1 ¡ x| ¤ C c hp1 |X ε 1 ¡ x|q. Hence, for h ¤ h 0 small enough, we deduce that |X ε 1 ¡ x| ¤ C c h and λ ε 1 c hG ¡1 pxqpζ ε p1{2q ¡ g T Y ε p1{2q qpxq R ε h , ( 3.17)
where |R ε h | ¤ Ch. For the term of size Ophq, we first deduce from (3.17) that § § §X

ε 1 ¡ x ¡ c hpY ε p1{2q gG ¡1 ζ ε p1{2q ¡ gG ¡1 g T Y ε p1{2q qpxq § § § ¤ Ch.
By using this estimate in (3.16), a Taylor expansion yields the desired expansion of λ ε 1 .

Proofs of the convergence theorems

Now that we have the local uniform expansion of the exact solution and the numerical scheme, as well as the stability property of Proposition 3.6, we are able to prove the main convergence theorems.

Proof of Theorem 2.7. We derive the global weak consistency (2.5)

E ε h ¤ N ņ1 § § ErϕpX ε pnh, X ε N ¡n pX 0 qqq ¡ ϕpX ε ppn ¡ 1qh, X ε N ¡n 1 pX 0 qqqs § § ¤ N ņ1 § § Erϕ n pX ε ph, X ε N ¡n pX 0 qqq ¡ ϕ n pX ε 1 pX ε N ¡n pX 0 qqqs § § ¤ N ņ1 Er § § Erϕ n pX ε ph, xqq ¡ ϕ n pX ε 1 pxqq|x X ε N ¡n pX 0 qs § § s ¤ N ņ1 Cp1 Er § § X ε N ¡n pX 0 q § § K sqh 3{2 ¤ Ch 1{2 ,
where we used the bounded moments property of Proposition 3.6 and X ε n M ε h (Lemma 3.1).

Proof of Theorem 2.9. We obtain straightforwardly from the expression of the integrator (2.2)

that |ζpX ε n q| ¤ C h c ε.
Using this estimate and the notation introduced in the proof of Lemma 3.8, we observe that § § §Y

ε piq pX ε n q ¡ Y 0 piq pX ε n q § § § ¤ C h c ε, § § §ζ ε piq pX ε n q § § § ¤ C h c ε, i t1{2, 1u,
where Y 0 p1{2q pxq σξ n and Y 0 p1q pxq f pxq. The Lagrange multiplier given by (3.16) therefore

satisfies § § §λ ε n 1 ¡ λ ε n 1 § § § ¤ C h c ε, where λ ε n 1 ¡G ¡1 X ε n 1 ¡X ε n pX ε n q » 1 0 g T pX ε n τ pX ε n 1 ¡ X ε n qqdτpσ c hξ n hf pX ε n qq.
Similarly to (3.16), the Lagrange multiplier λ 0 n 1 of the constrained Euler integrator (1.6) is given by

λ 0 n 1 ¡G ¡1 X 0 n 1 ¡X 0 n pX 0 n q » 1 0 g T pX 0 n τ pX 0 n 1 ¡ X 0 n qqdτpσ c hξ n hf pX 0 n qq.
Using that G ¡1 y pxq is bounded, a straightforward calculation shows that G ¡1 y pxq is Lipschitz in x, y R d , that is, there exists a constant L ¡ 0 such that § § G ¡1

y 1 px 1 q ¡ G ¡1 y 2 px 2 q § § ¤ Lp|x 1 ¡ x 2 | |y 1 ¡ y 2 |q, x 1 , x 2 , y 1 , y 2 R d .
Thus, we get the estimate

§ § λ ε n 1 ¡ λ 0 n 1 § § ¤ C h p c ε § § X ε n ¡ X 0 n § § q C c h § § X ε n 1 ¡ X 0 n 1 § § ,
and, as λ ε n 1 and λ 0 n 1 are bounded uniformly in ε, we have

§ § gpX ε n qλ ε n 1 ¡ gpX 0 n qλ 0 n 1 § § ¤ C h p c ε § § X ε n ¡ X 0 n § § q C c h § § X ε n 1 ¡ X 0 n 1 § § .
From the definitions of X ε n 1 in (2.2) and X 0 n 1 in (1.6), we deduce that

§ § X ε n 1 ¡ X 0 n 1 § § ¤ C h p c ε § § X ε n ¡ X 0 n § § q C c h § § X ε n 1 ¡ X 0 n 1 § § . If h ¤ h 0 is small enough, we obtain § § X ε n 1 ¡ X 0 n 1 § § ¤ C h p c ε § § X ε n ¡ X 0 n § § q.
We deduce the estimate (2.6) by iterating this inequality for n ¤ N T {h.

Proof of Theorem 2.12. We take over the notations and the expression (3.16) of the Lagrange multiplier λ ε 1 that we used in the proof of Lemma 3.8. Without loss of generality, we concentrate on the first step of the algorithm with X 0 x. Replacing λ ε 1 by the explicit formula (3.16) in (2.2) yields that X ε 1 is a fixed point of the following map:

F ε h pyq x c hY ε p1{2q pxq hY ε p1q pxq gpxqG ¡1 y¡x pxq c hζ ε p1{2q pxq hζ ε p1q pxq ¡ » 1 0 g T px τ py ¡ xqqdτ p c hY ε p1{2q pxq hY ε p1q pxqq % . For y 1 , y 2 R d , Assumption 2.1 gives |F ε h py 2 q ¡ F ε h py 1 q| ¤ C c h |y 2 ¡ y 1 | ,
where we used that G ¡1 y pxq is Lipschitz in x, y R d . We deduce that F ε h is a uniform contraction for h ¤ h 0 small enough.

Numerical experiments

In this section, we perform numerical experiments to confirm the theoretical findings, on a torus in R 3 and on the orthogonal group in high dimension and codimension, in the spirit of the experiments in [START_REF] Zappa | Monte Carlo on manifolds: sampling densities and integrating functions[END_REF][START_REF] Zhang | Ergodic SDEs on submanifolds and related numerical sampling schemes[END_REF][START_REF] Laurent | Order conditions for sampling the invariant measure of ergodic stochastic differential equations on manifolds[END_REF].

Uniform approximation for the invariant measure on a torus

We consider the example in codimension one of a torus in R 3 . We apply the new method given by the discretization (2.4) for sampling the invariant measure of (1.4) for different steps h and parameters ε, and we compare it with the Euler integrators (1.5) in R d and (1.6) on the manifold. We recall that for dynamics of the form (1.4), the weak convergence implies the convergence for the invariant measure (see [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF]). The numerical experiments of this subsection hint that the uniform accuracy property also extends to the approximation of the invariant measure, but we leave the mathematical analysis for the invariant measure for future work. We consider the constraint ζpxq px 2

1 x 2 2 x 2 3 R 2 ¡r 2 q 2 ¡4R 2 px 2 1 x 2
2 q, with R 3 and r 1, and we choose the map f pxq ¡25px 1 ¡R r, x 2 , x 3 q, with σ c 2, the test function ϕpxq |x| 2 , M 10 7 trajectories, the final time T 10, and the initial condition X 0 pR ¡ r, 0, 0q.

Increasing the value of T does not modify the computed averages, which hints that we reached the equilibrium. The factor 25 in f confines the solution in a reasonably small neighborhood of the torus, which allows a faster convergence to equilibrium and to take fewer trajectories. We compute the Monte-Carlo estimator s J 1 M °M m1 ϕpX pmq N q ErϕpX N qs, where X pmq n is the mth realization of the integrator at time t n nh, and N is an integer satisfying N h T . We compare this approximation with a reference value of V . We observe in Figure 1 that the accuracy of the constrained integrator (1.6) for solving the unconstrained problem (1.4) deteriorates when ε grows larger, as µ ε V deviates from µ 0 V . The explicit Euler scheme (1.5) faces stability issues when ε Ñ 0. The accuracy of the new method for solving (1.4) does not deteriorate depending on ε, and it shares a behavior similar to the constrained Euler scheme (1.6) when ε Ñ 0, which is in agreement with Theorem 2.9. The right graph of Figure 1 shows that the behavior of the error in ε is the same for any fixed value of h. This is a numerical confirmation of the uniform accuracy property of the discretization (2.2) (in the spirit of the numerical experiments in [START_REF] Chartier | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF]), as stated in Theorem 2.7. For any fixed ε, the right graph of Figure 1 also shows that the error decreases when h Ñ 0. A plot of the error against h for a fixed ε (not included for conciseness) shows a slope of order one (see Remark 2.8).

Figure 1: Error for sampling the invariant measure of penalized Langevin dynamics on a torus in R 3 of the uniform discretization (2.2) and the Euler integrators (1.5) and (1.6) for different values of ε with h 2 ¡9 T (left), and error curves versus ε of the uniformly accurate method for different timesteps h 2 ¡i T and i 6, . . . , 10 (right), with the final time T 10, the maps f pxq ¡25px 1 ¡ R r, x 2 , x 3 q, ϕpxq |x| 2 , and M 10 7 trajectories.

Weak approximation on the orthogonal group

We apply the uniformly accurate method on a compact Lie group (in the spirit of the numerical experiments in [START_REF] Zappa | Monte Carlo on manifolds: sampling densities and integrating functions[END_REF][START_REF] Zhang | Ergodic SDEs on submanifolds and related numerical sampling schemes[END_REF]) to see how it performs in high dimension and codimension. We choose the orthogonal group Opmq tM R m¢m , M T M I m u, seen as a submanifold of R m 2 of codimension q mpm 1q{2. We compare the explicit Euler scheme (1.5), the constrained Euler scheme (1.6), and the new method on M Opmq for m 2, . . . , 5 with the parameters ε 0.005, T 1, and h 2 ¡7 . Note that, as h and ε share the same order of magnitude, the explicit Euler scheme (1.5) can face stability issues, and the solution does not lie on the manifold M. Thus, we are in the regime where the convergence results for both Euler schemes (1.5)-(1.6) do not apply. We choose f ¡∇V , where V is given by V pxq 50 Trppx ¡ I m 2 q T px ¡ I m 2 qq with the parameters σ c 2, X 0 I d , ϕpxq Trpxq, and M 10 6 trajectories. The reference solution for Jpmq ErϕpX ε pTqqs is computed with the uniformly accurate integrator with h ref 2 ¡9 . With the factor 50 in the potential (4.1), the trajectories stay close to I m 2 , and Jpmq is close to ϕpI m 2 q m. This choice of factor permits one to explore a reasonably small area of Opmq, in order to avoid zones close to M where the Gram matrix G has a bad condition number or is singular, and to reduce the number of trajectories needed. We observe numerically that replacing the factor 50 by 1 in (4.1) induces a severe timestep restriction. We present the results of the experiment in Table 1. We omit the results for the explicit Euler scheme (1.5) as the method is inaccurate in this regime (error of size 1). We observe that, in the regime where h and ε share the same order of magnitude, the uniformly accurate integrator performs significantly better than the Euler schemes (1.6) and (1.5) for solving the problem (1.4). In this regime, the Euler method (1.5) faces stability issues, and it is inappropriate to use the constrained Euler scheme (1.6) as the solution X ε ptq of (1.4) is not close to the solution X 0 ptq of (1.2). Moreover, the cost in time of the new method stays the same in average for any value of ε (results not included for conciseness). This confirms numerically the uniform cost of solving the fixed point problem (2.2), as stated in Theorem 2.12. J M ¡1 °M k1 ϕpX pkq N q, where pX n q is given by the uniform discretization (2.2) for s J UA and by the constrained Euler scheme (1.6) for s J EC with their respective errors. The average is taken over M 10 6 trajectories with the potential (4.1), ϕpxq Trpxq, the final time T 1, the stiff parameter ε 0.005, and the timestep h 2 ¡7 .

m

Conclusion and future work

In this work, we presented a new method for the weak numerical integration of penalized Langevin dynamics evolving in the vicinity of manifolds of any dimension and codimension.

On the contrary of the other existing discretizations, the accuracy of the proposed integrator is independent of the size of the stiff parameter ε. Moreover, its cost does not depend on ε, and it converges to the Euler scheme on the manifold when ε Ñ 0. Throughout the analysis, we gave an expansion in time of the solution to the penalized Langevin dynamics that is uniform in ε, as well as new tools for the study of stochastic projection methods for solving stiff SDEs.

Multiple questions arise from the work presented in this paper, with many of great interest for physical applications. First, it would be interesting to get convergence results with weaker assumptions or to develop uniformly accurate integrators for different penalized dynamics with the same limit when ε Ñ 0 such as the original penalized dynamics (1.3) (see [START_REF] Ciccotti | Projection of diffusions on submanifolds: application to mean force computation[END_REF]). One could build integrators for penalized dynamics of the form

dX ε f pX ε qdt σdW σ 2 4 ∇ lnpdetpGqq ¡ 1 ε pgG ¡1 ζ 1 qpX ε qdt ¡ 1 ν pgG ¡1 ζ 2 qpX ε qdt,
where ε and ν do not share the same order of magnitude, or for constrained dynamics with a penalized term. One could also build a uniformly accurate numerical scheme with high order in the weak context, or just in the context of the invariant measure (in the spirit of the works [START_REF] Bou-Rabee | Long-run accuracy of variational integrators in the stochastic context[END_REF][START_REF] Leimkuhler | Rational construction of stochastic numerical methods for molecular sampling[END_REF][START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF][START_REF] Abdulle | Long time accuracy of Lie-Trotter splitting methods for Langevin dynamics[END_REF][START_REF] Leimkuhler | The computation of averages from equilibrium and nonequilibrium Langevin molecular dynamics[END_REF][START_REF] Laurent | Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs[END_REF][START_REF] Laurent | Order conditions for sampling the invariant measure of ergodic stochastic differential equations on manifolds[END_REF] where numerical schemes of high order for the invariant measure and weak order one were introduced). Postprocessors [START_REF] Vilmart | Postprocessed integrators for the high order integration of ergodic SDEs[END_REF] proved to be an efficient tool for reaching high order for the invariant measure without increasing the cost of the method and could be used in this context. Moreover, the order conditions presented in [START_REF] Laurent | Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs[END_REF] for Runge-Kutta methods for solving Langevin dynamics in R d both in the weak sense and for the invariant measure do not match with the order conditions for solving Langevin dynamics constrained on the manifold M, as presented in [START_REF] Laurent | Order conditions for sampling the invariant measure of ergodic stochastic differential equations on manifolds[END_REF]. It would be interesting to create a unified class of high order Runge-Kutta methods with the same order conditions in R d , on the manifold M and in the vicinity of M. The discretizations presented in this paper could also be combined with Metropolis-Hastings rejection procedures [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF][START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF], in the spirit of the works [START_REF] Girolami | Riemann manifold Langevin and Hamiltonian Monte Carlo methods[END_REF][START_REF] Brubaker | A family of MCMC methods on implicitly defined manifolds[END_REF][START_REF] Lelièvre | Langevin dynamics with constraints and computation of free energy differences[END_REF][START_REF] Zappa | Monte Carlo on manifolds: sampling densities and integrating functions[END_REF][START_REF] Lelièvre | Hybrid Monte Carlo methods for sampling probability measures on submanifolds[END_REF], in order to get an exact approximation for the invariant measure with a rejection rate that does not deteriorate in the regime ε Ñ 0.

2

 2

Assumption 2 . 4 . 1 0g

 241 There exists c ¡ 0 such that, for x, y R d , |G y pxq| ¥ c, where G y pxq » T px τ yqdτ gpxq.

³R

  scheme(1.6), it amounts to comparing an approximation of³ M ϕdµ 0V with the reference value

  σ

  Rq, the vector space of C V functions ϕpxq such that all partial derivatives up to any order have a polynomial growth of the form § § ϕ pkq pxq § § ¤ Cp1 |x| K q, where the constants C and K are independent of x R d (but can depend on k), and where we denote by |x| px T xq 1{2 the Euclidean norm in R d . Similarly, we denote by C p P pR d , Rq the space of C p functions whose partial derivatives up to order p have polynomial growth.

	functions typically belong to a subspace of C V	The test
	Letting φ : R d Ñ R d¢k , we use the following notations for differentials: for all vectors x, a 1 , . . . , a m R d , we denote
	φ pmq pxqpa 1 , . . . , a m q	d i1
		,...,im1

P pR d ,

  R d Ñ R d¡q satisfies φ I g 0. The map ψ lies in C 5 P and is invertible, ψ I and ψ P are Lipschitz and there exist two constants c, C ¡ 0 such that c ¤ |ψ I pxq| ¤ C.

	Appx. C] and [40, Chap. 2].
	Assumption 2.1. The map f is bounded, is Lipschitz and lies in C 3 P . The maps g and g I are bounded in R d , and there exists c ¡ 0 such that, for x, y R d ,
	|G y pxq| ¥ cp1 |y|q ¡1 , where G y pxq	» 1	g T px τ yqdτ gpxq.	(2.1)
		0
	In addition, there exists a smooth change of coordinate
	ψ : R d Ñ R d ¢ φpxq x Þ Ñ ζpxq
	where φ : Remark 2.2. Assumption 2.1 is almost the same as the one given in [16, Appx. C] to prove
	the strong convergence of the dynamics (1.3) to the constrained dynamics (1.2). The differ-
	ence lies in the additional estimate (2.1) that replaces the weaker assumption that the Gram matrix Gpxq G 0 pxq is invertible on M. We use the estimate (2.1) for obtaining a uniform
	expansion of the Lagrange multipliers in the new method and for proving that the new method
	evolves in a neighborhood of the manifold (see Lemma 3.8). Note that the existence of the
	change of coordinate ψ is always valid in a neighborhood of the smooth manifold M. The same
	goes for the estimate (2.1) for x in a neighborhood of M and y in a ball centered on zero.
	Assumption 2.1 is valid in particular if M is a vector subspace of R d , but it is quite restrictive.
	It would be interesting to extend the results of this paper under simpler regularity assumptions
	made only on the manifold, as numerical experiments hint that the results presented in this
	paper still stand without global assumptions. This is matter for future work.

  ζ, gG ¡1 ζq 0.

	Applying Lemma 3.2, we get by a direct calculation that
	φpY ε phqq φpxq hφ I f pxq h σ 2 4	φ I ∇ lnpdetpGqqpxq h	σ 2 2	∆φpxq R ε h pxq,
	which gives the desired estimate
	§ § §ErφpX

  with techniques similar to the ones presented in [40, Chap. 2]. We denote by X ε pt, xq the solution of the penalized dynamics with initial condition x and X ε n pxq the numerical solution with initial condition x. For x M ε h , Proposition 3.3 and Proposition 3.7 yield |ErϕpX ε ph, xqq ¡ ϕpX ε 1 pxqq|xs| ¤ Cp1 |x| K qh 3{2 , where ϕ C 3 P . Lemma 3.5 gives that ϕ n pxq ErϕpX ε ppn ¡ 1qh, xqq|xs is in C 3 P . We rewrite the global error, given by E ε h |ErϕpX ε pT, X 0 qq ¡ ϕpX ε N pX 0 qqs|, with a telescopic sum,

Table 1 :

 1 ¤ 10 ¡2 4.91298 1.1 ¤ 10 ¡1 Numerical approximation of Jpmq ErϕpX ε pTqqs for 2 ¤ m ¤ 5 with the estima-

		dimpMq q	Jpmq	s J UA	error of s J UA	s J EC	error of s J EC
	2 3 4	1 3 6	3 6 10	2.00934 3.01458 4.02050	2.00619 3.00821 4.00972	3.1 ¤ 10 ¡3 1.99165 1.8 ¤ 10 ¡2 6.4 ¤ 10 ¡3 2.97460 4.0 ¤ 10 ¡2 1.1 ¤ 10 ¡2 3.94846 7.2 ¤ 10 ¡2
	5	10	15	5.02669	5.00842	1.8

tor s

  2 4 pφ I ∇ lnpdetpGqqqpX ε ptqqdt σ 2 2 ∆φpX ε ptqqdt,dζpX ε ptqq pg T f qpX ε ptqqdt σg T pX ε ptqqdWptq σ 2 4 pg T ∇ lnpdetpGqqqpX ε ptqqdtwhere we used that φ I g 0 and that ζ I g T . Therefore, ζpX ε ptqq satisfiesζpX ε ptqq σ pX ε psqqdt,and a bound on ζpX ε ptqq is given byEr|ζpX ε ptqq| 2 s ¤ COn the other hand, we rewrite (1.2) in the orthogonal coordinates asdφpX 0 ptqq pφ I f qpX 0 ptqqdt σφ I pX 0 ptqqdW σ 2 4 pφ I ∇ lnpdetpGqqqpX 0 ptqqdt σ 2 Er § § φpX ε ptqq ¡ φpX 0 ptqq § § 2 s ¤ C Er § § X ε psq ¡ X 0 psq § § 2 sds,and from (A.1) we deduce thatEr § § ζpX ε ptqq ¡ ζpX 0 ptqq § § 2 s ¤ Cε.Thus, ψpX ε ptqq satisfiesEr § § ψpX ε ptqq ¡ ψpX 0 ptqq § § 2 s ¤ C Er § § X ε psq ¡ X 0 psq § § 2 sds Cε,and, as ψ ¡1 is Lipschitz, we haveEr § § X ε ptq ¡ X 0 ptq § § 2 s ¤ C Er § § X ε psq ¡ X 0 psq § § 2 sds Cε,which, with the use of the Gronwall lemma, gives the desired estimate.

		σ 2 2	divpgqpX ε ptqqdt ¡ 1 ε	ζpX ε ptqqdt,
			» t	e ps¡tq{ε g T pX ε psqqdW
			» t 0	0 e ps¡tq{ε	g T f		σ 2 4	g T ∇ lnpdetpGqq	σ 2 2	divpgq
						» t	e 2ps¡tq{ε dt C	¡ » t	e ps¡tq{ε dt	© 2 ¤ Cε.	(A.1)
							0	0
		dζpX 0 ptqq 0,			2	∆φpX 0 ptqqdt,
	where we used that
	1 2	φ I ∇ lnpdetpGqq ∆φ	d i1	φ I pΠ I
								» t
								0
								» t
								0
								» t
								0

% M pΠ M e i qe i q d i1 φ P pΠ M e i , Π M e i q.

With Assumption 2.1, we obtain

Department of Mathematics, University of Bergen, Norway. Adrien.Laurent@uib.no.

Acknowledgements. The author would like to thank Gilles Vilmart for helpful discussions. This work was partially supported by the Swiss National Science Foundation, grants No. 200020_184614, No. 200021_162404 and No. 200020_178752. The computations were performed at the University of Geneva on the Baobab cluster using the Julia programming language.

Appendices A Proof of Theorem 2.3

In this section, we prove the convergence of the penalized Langevin dynamics (1.4) to the constrained dynamics (1.2) when ε Ñ 0, as stated in Theorem 2.3. The proof uses techniques and arguments similar to those in [16, Appx. C]. However, since we rescaled the stiff term in (1.4), there is no need for a change of time to prove the convergence to the constrained dynamics.

Proof of Theorem 2.3. In the orthogonal coordinates system given by Assumption 2.1, equation (1.4) becomes dφpX ε ptqq pφ I f qpX ε ptqqdt σφ I pX ε ptqqdWptq