Grégoire Ferré
email: gregoire.ferre@enpc.fr

Hugo Touchette
email: htouchette@sun.ac.za

Adaptive sampling of large deviations

We introduce and test an algorithm that adaptively estimates large deviation functions characterizing the fluctuations of additive functionals of Markov processes in the long-time limit. These functions play an important role for predicting the probability and pathways of rare events in stochastic processes, as well as for understanding the physics of nonequilibrium systems driven in steady states by external forces and reservoirs. The algorithm uses methods from risk-sensitive and feedback control to estimate from a single trajectory a new process, called the driven process, known to be efficient for importance sampling. Its advantages compared to other simulation techniques, such as splitting or cloning, are discussed and illustrated with simple equilibrium and nonequilibrium diffusion models.

I. INTRODUCTION

.

.

, which we adapt to continuous-time diffusion processes and to the problem of estimating large deviation functions. The algorithm works by estimating or learning "on the fly" a modified process, called the auxiliary or driven process, which corresponds to the process that is asymptotically

equivalent to the original process conditioned on the rare event of interest [START_REF] Chetrite | Nonequilibrium Markov processes conditioned on large deviations[END_REF] or, alternatively, to the exponential tilting of the original process, known to be efficient for importance sampling [START_REF] Chetrite | Variational and optimal control representations of conditioned and driven processes[END_REF]. This modified process is given by a principal eigenvalue problem related to the Feynman-Kac equation or, equivalently, by a stochastic optimal control problem [START_REF] Chetrite | Variational and optimal control representations of conditioned and driven processes[END_REF] that we solve iteratively using stochastic approximation and feedback control methods.

The main advantage of this algorithm, compared to splitting or cloning, is that it does not require the simulation of many copies of the process considered -it runs on one long trajectory of that process, modified with a feedback-reinforcement rule, to adaptively learn the driven process, thereby reducing significantly the complexity of estimating large deviation functions. The calculation of error bars for the estimated quantities is also simplified compared to other techniques, as the algorithm is based on simple time averages and stochastic approximations [START_REF] Chauveau | Estimation of the asymptotic variance in the CLT for Markov chains[END_REF][START_REF] Roberts | General state space Markov chains and MCMC algorithms[END_REF][START_REF] Benveniste | Adaptive Algorithms and Stochastic Approximations[END_REF]. Finally, the errors incurred by discretizing continuous-time processes and functionals can be analysed in a precise way, in principle, via Feynman-Kac semi-groups [START_REF] Ferré | Error estimates on ergodic properties of discretized Feynman-Kac semigroups[END_REF].

We discuss these advantages and test the algorithm in Sec. IV with simple equilibrium and nonequilibrium diffusions, after introducing the general model and notations in Sec. II and the algorithm in Sec. III. The results at this point are preliminary and are presented as a proof of concept of the algorithm. More detailed results about the time-discretization and sampling errors will be addressed in future works, together with more complex applications involving interacting particle systems and higher-dimensional diffusions.

II. FRAMEWORK A. Model and notations

We consider an ergodic diffusion (X t ) t≥0 evolving in a state space X ⊂ R d according to the following stochastic differential equation (SDE):

dX t = b(X t ) dt + σ dW t , (1)
where b : X → R d is a smooth function, called the drift, W t is an m-dimensional Brownian motion, and σ is a d × m matrix, assumed to be constant for simplicity. 1 The generator of this diffusion reads

L = b • ∇ + 1 2 ∇ • D∇, (2)
where • denotes the scalar product and D = σσ T , with T as the transpose, is the diffusion matrix, assumed to be positive definite. This is the generator of the evolution semi-group P T , defined by

P T ϕ(x) = E [ϕ(X T ) | X 0 = x] , (3)
for all time T ≥ 0 and any smooth test function ϕ. The dual L † of L in the space L 2 (dx) of squareintegrable functions with respect to the Lebesgue measure is the generator of the Fokker-Planck equation

∂ t ρ(x, t) = L † ρ(x, t), (4)
which gives the evolution of the probability density ρ(•, t) of X t starting from some initial density ρ(•, 0) for X 0 .

Our goal here is to study the fluctuations of time-integrated functionals of X t , called observables, having the general form

A T = 1 T T 0 f (X t ) dt + 1 T T 0 g(X t ) • dX t , (5)
where f : X → R and g : X → R d are reasonably smooth functions (e.g. continuous) and • denotes the Stratonovich product [START_REF] Pavliotis | Stochastic Processes and Applications[END_REF]. Such a functional defined over the time horizon [0, T ] can represent, for example, a control cost associated with the state X t and its increments [START_REF] Chernyak | Stochastic optimal control as nonequilibrium statistical mechanics: Calculus of variations over density and current[END_REF] or a physical quantity integrated in time, such as the work performed on a particle by external forces or the heat exchanged by a particle with its environment [START_REF] Sekimoto | Stochastic Energetics[END_REF].

Assuming that the process is ergodic with respect to an invariant measure µ(dx) = ρ * (x) dx with smooth density ρ * , we have

A T T →∞ ----→ X f (x)ρ * (x) dx + X g(x) • J * (x) dx = a * , (6)
almost surely, where

J * (x) = b(x)ρ * (x) - D 2 ∇ρ * (x) (7)
is the stationary current field associated with ρ * [START_REF] Chetrite | Variational and optimal control representations of conditioned and driven processes[END_REF]. The theory of large deviations [START_REF] Dembo | Large Deviations Techniques and Applications[END_REF] refines this ergodic theorem, generalized here with the additional g term, by providing estimates for the rate at which the probability distribution of A T concentrates on its ergodic value a * . Such estimates can be derived under general conditions (see [START_REF] Dembo | Large Deviations Techniques and Applications[END_REF]) and take, in the simplest case, the form lim

T →∞ - 1 T log P(A T ∈ B) = min a∈B I(a), (8)
for any Borel subset B of R, where I : R → R + is a positive function such that I(a * ) = 0. When this limit exists, A T is said to satisfy the large deviation principle (LDP) with rate function I. Formally, this means that

P(A T ∈ da) = e -T I(a)+o(T ) da, (9)
where o(T ) denotes corrections in the exponential that grow slower than linearly in T . Thus, we see that the rate function provides useful information about the fluctuations of A T : the likelihood that A T = a decays exponentially with time for all a = a * , since I(a) > 0 in this case, and converges otherwise to 1 as T → ∞, since I(a * ) = 0. Moreover, the rate function is in general not a parabola, meaning that it describes fluctuations that are generally not Gaussian. In this sense, large deviation theory is often seen as an extension of both the ergodic theorem, which describes the concentration of A T towards its mean, and the central limit theorem, which describes the local Gaussian fluctuations of A T around its mean [START_REF] Dembo | Large Deviations Techniques and Applications[END_REF].

B. Large deviation functions and driven process

In practice, the rate function I(a) can be calculated in many different ways other than by direct sampling, which requires exponentially large samples with T [START_REF] Bucklew | Introduction to Rare Event Simulation[END_REF]. The most common method proceeds from the scaled cumulant generating function (SCGF), defined for k ∈ R by

λ(k) = lim T →∞ 1 T log E e kT A T . ( 10
)
By the Gärtner-Ellis Theorem [START_REF] Dembo | Large Deviations Techniques and Applications[END_REF], the Legendre-Fenchel transform of this function yields the rate function:

I(a) = sup k∈R {ka -λ(k)}, (11)
provided, essentially, that λ(k) is differentiable; see [START_REF] Dembo | Large Deviations Techniques and Applications[END_REF] for more precise conditions. For the SDE (1) and the additive functional (5), the SCGF is known to be given by the principal eigenvalue of the operator

L k = b • (∇ + kg) + 1 2 (∇ + kg) • D(∇ + kg) + kf, (12)
which is the generator of the Feynman-Kac semi-group P k T , defined by

P k T ϕ(x) = E ϕ(X T ) e T kA T X 0 = x , (13)
for any smooth function ϕ; see Appendix A.2 of [START_REF] Chetrite | Nonequilibrium Markov processes conditioned on large deviations[END_REF]. In the end, the rate function can thus be computed by solving the spectral problem

L k h k = λ(k)h k , (14)
where λ(k) is the principal eigenvalue of L k and h k its corresponding eigenvector. This holds provided that this operator has reasonable spectral properties, made precise in the following assumption.

Assumption 1. The operator L k defined in [START_REF] Hedges | Dynamic order-disorder in atomistic models of structural glass formers[END_REF] acting on L 2 (µ) has an isolated largest eigenvalue λ(k). Its multiplicity is one and it is associated with a regular eigenvector

h k ∈ L 2 (µ) such that for all x ∈ X , h k (x) > 0.
This assumption holds for many systems, in particular when X is bounded [START_REF] Ferré | Error estimates on ergodic properties of discretized Feynman-Kac semigroups[END_REF][START_REF] Bierkens | Linear PDEs and eigenvalue problems corresponding to ergodic stochastic optimization problems on compact manifolds[END_REF] or when b and g are gradient fields with appropriate growth conditions; see [START_REF] Lelièvre | Partial differential equations and stochastic methods in molecular dynamics[END_REF]Sec. 2.5]. In practice, the spectral problem ( 14) can be solved numerically using standard projection or discretization (Galerkin) methods, which work well for low-dimensional systems [START_REF] Chatelin | Spectral Approximation of Linear Operators[END_REF], or more involved real space renormalization methods when dealing with higher-dimensional systems [START_REF] Gorissen | Finite size scaling of current fluctuations in the totally asymmetric exclusion process[END_REF]. Note that, for g = 0, L k is the usual Feynman-Kac generator with source term kf . Moreover, for k = 0, L 0 = L is the generator of the SDE, so that λ(0) = 0 and h 0 = 1, the constant unit function.

The numerical method that we propose in the next section attempts to estimate the spectral elements λ(k) and h k in a different way using the fact that they are solutions of the family of eigenproblems

P k T h k = e T λ(k) h k , ∀ T > 0, (15)
which can be approximated stochastically. The method also exploits a connection between large deviations and control theory showing that λ(k) is the ergodic limit of an optimal control cost satisfying a stochastic Hamilton-Jacobi-Bellman equation [START_REF] Fleming | Controlled Markov Processes and Viscosity Solutions[END_REF], and that h k determines the controlled diffusion achieving the optimal cost, which is the driven process mentioned in the introduction. These results are explained in detail in [START_REF] Chetrite | Variational and optimal control representations of conditioned and driven processes[END_REF] (see also references therein); for the purpose of this paper, we only state them without proofs.

The controlled diffusion, denoted by ( Xt ) t≥0 , satisfies the SDE

d Xt = b k ( Xt ) dt + σ dW t , (16)
where

b k (x) = b(x) + D[kg(x) + ∇ log h k (x)] (17)
is the optimal control drift defined for all x ∈ X . Under Assumption 1, this modified diffusion is ergodic with respect to a new invariant measure µ k , whose density is

ρ k (x) = h k (x)l k (x), (18)
where l k is the dual of

h k in L 2 (dx) satisfying L † k l k = λ(k)l k [6, Sec. 3.3]. These two functions are normalized such that X h k (x)l k (x) dx = 1, X l k (x) dx = 1. ( 19
)
Because of the change of process, the observable A T must concentrate on a new value, which can be shown to be given by λ (k), that is,

A T T →∞ ----→ λ (k) = a k , (20)
almost surely with respect to ( Xt ) t≥0 . Similarly, the control cost

C T = kA T -R T , (21)
where

R T = 1 2T T 0 [b( Xt ) -b k ( Xt )] • D -1 [b( Xt ) -b k ( Xt )]dt, (22)
reaches in the ergodic limit the value λ(k), so that

C T T →∞ ----→ λ(k) (23)
almost surely under ( Xt ) t≥0 [START_REF] Chetrite | Variational and optimal control representations of conditioned and driven processes[END_REF]Sec. 4]. Finally, it can be shown by Legendre duality that the rate function at the value a k = λ (k) is given by the ergodic limit of R T above, leading to

R T T →∞ ----→ I(a k ), (24)
almost surely with respect to ( Xt ) t≥0 [START_REF] Chetrite | Variational and optimal control representations of conditioned and driven processes[END_REF]Sec. 4]. Note that the diffusion is not modified for k = 0, so that a 0 = a * , λ (0) = a * , and I(a * ) = 0. These ergodic limits provide direct estimators of the SCGF and the rate function, based on a single trajectory of the driven process, which can be simulated for different values of the parameter k ∈ R. For the SCGF, there are in fact three possible estimators:

1. From (20): The value of A T , integrated numerically with the condition λ(0) = 0.

From (23):

The value of C T .

3. The eigenvalue returned by the algorithm proposed in Sec. III.

In practice, we find that the last estimator is more stable, although the first and second are more adapted to obtain error bars.

For the rate function, we have two possible estimators:

1. The Legendre transform [START_REF] Garrahan | Dynamical first-order phase transition in kinetically constrained models of glasses[END_REF] of the SCGF, given in parametric form by

I(a k ) = ka k -λ(k), (25)
where a k is either estimated from A T or by taking the numerical derivative of λ(k).

2. From [START_REF] Bolhuis | Transition path sampling: Throwing ropes over rough mountain passes, in the dark[END_REF]: The value of R T obtained for the value of A T , giving the couple (A T , R T ).

In practice, we find that the first estimator based on the Legendre transform is more reliable. By comparison, the computation of I(a) based on [START_REF] Juneja | Rare-event simulation techniques: An introduction and recent advances[END_REF] involves the optimal drift b k and, therefore, the logarithmic derivative of h k , which is more difficult to estimate in a stable way. In all cases, error bars can be constructed from the same trajectory by estimating, in principle, the variance of A T , C T and R T using covariance techniques for Markov processes [START_REF] Chauveau | Estimation of the asymptotic variance in the CLT for Markov chains[END_REF][START_REF] Roberts | General state space Markov chains and MCMC algorithms[END_REF][START_REF] Benveniste | Adaptive Algorithms and Stochastic Approximations[END_REF]. This is an advantage over splitting and cloning algorithms, for which the calculation of errors bars is difficult, as they involve correlated copies or "clones" of the simulated process [START_REF] Rousset | On the control of an interacting particle estimation of Schrödinger ground states[END_REF][START_REF] Nemoto | Population-dynamics method with a multicanonical feedback control[END_REF][START_REF] Nemoto | Finite-time and finite-size scalings in the evaluation of large-deviation functions: Analytical study using a birth-death process[END_REF].

In closing, it is interesting to note that the driven process can also be interpreted as the change of process in the importance sampling of the probability P(A T ∈ B) that is optimal in the sense of logarithmic or asymptotic efficiency [START_REF] Bucklew | Introduction to Rare Event Simulation[END_REF]. Therefore, it can be used not only to estimate the SCGF and rate function, but also to estimate the actual probability P(A T ∈ B) in an efficient way [START_REF] Borkar | Peformance analysis conditioned on rare events: An adaptive simulation scheme[END_REF]. The optimal change of process in this case is known to correspond to the exponential tilting of the original process [START_REF] Bucklew | Introduction to Rare Event Simulation[END_REF], which is a time-dependent process in general; see Appendix D of [START_REF] Chetrite | Nonequilibrium Markov processes conditioned on large deviations[END_REF]. In the ergodic limit, this process converges to a homogeneous process corresponding exactly to the driven process [START_REF] Cérou | Adaptive multilevel splitting for rare event analysis[END_REF]. We refer to [START_REF] Chetrite | Variational and optimal control representations of conditioned and driven processes[END_REF] for more details about these results.

III. ADAPTIVE ALGORITHM

We are now ready to present the algorithm for estimating the SCGF and the rate function of A T for continuous-time diffusions. The algorithm is based, as mentioned, on prior algorithms proposed in [START_REF] Borkar | Peformance analysis conditioned on rare events: An adaptive simulation scheme[END_REF][START_REF] Ahamed | Adaptive importance sampling technique for Markov chains using stochastic approximation[END_REF][START_REF] Basu | A learning algorithm for risk-sensitive cost[END_REF][START_REF] Borkar | Learning algorithms for risk-sensitive control[END_REF] for Markov chains and exploits the fact that e -T λ(k) (P k T ϕ)

T →∞ ----→ h k X ϕ(x) l k (x) dx, (26)
for any smooth test function ϕ. This result follows under Assumption 1; see [6, Sec. III.B] or [START_REF] Ferré | Error estimates on ergodic properties of discretized Feynman-Kac semigroups[END_REF]Sec. 6.1] for a proof. The algorithm that we propose works from this limit by approximating the action of the Feynman-Kac semi-group P k T in a stochastic way as a time average computed over a long trajectory of the diffusion. Moreover, it continuously modifies the diffusion as we estimate h k to construct the driven process ( Xt ) t≥0 , which underlies the estimators of the large deviation functions.

The algorithm is presented next. Our contribution compared to [START_REF] Borkar | Peformance analysis conditioned on rare events: An adaptive simulation scheme[END_REF][START_REF] Ahamed | Adaptive importance sampling technique for Markov chains using stochastic approximation[END_REF][START_REF] Basu | A learning algorithm for risk-sensitive cost[END_REF][START_REF] Borkar | Learning algorithms for risk-sensitive control[END_REF] is to consider general time-continuous processes, time-additive functionals of these processes that depend on both their state and increments, and, more importantly, to construct the driven process explicitly so as to estimate the SCGF and the rate function adaptively using the estimators introduced in the previous section.

A. Time discretization

The first step required in the algorithm is to discretize in time the SDE (1) and its associated Feynman-Kac semi-group ( 13) by transforming the Markov diffusion (X t ) t≥0 into a Markov chain (x n ) n∈N with a small time step. Many discretization schemes can be used for the SDE; see, for instance, [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF]. Here we use the standard Euler-Maruyama scheme with constant time step ∆t, given by

x n+1 = x n + b(x n )∆t + σ √ ∆t ξ n , (27)
where (ξ n ) n≥0 is a sequence of independent standard d-dimensional Gaussian random variables.

The corresponding discretization of the evolution semi-group P T over a time step ∆t is denoted by Q ∆t , so that

Q ∆t ϕ(x) = E [ϕ(x n+1 )|x n = x] , (28)
for any test function ϕ and x ∈ X . We refer to [START_REF] Lelièvre | Partial differential equations and stochastic methods in molecular dynamics[END_REF][START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF] for more information about the discretization of SDEs and their weak error analysis. Many discretizations also exist for the Feynman-Kac semi-group P k T . Here, we use the natural scheme where the diffusion is discretized as above and the integral of A T is discretized as a Riemann sum with the left-point rule for the integral involving f and the mid-point rule for the Stratonovich integral involving g [START_REF] Ferré | Error estimates on ergodic properties of discretized Feynman-Kac semigroups[END_REF]. The action of P k T is thus replaced by

Q k ∆t ϕ(x) = E e k f (xn)∆t+g x n+1 +xn 2 •(x n+1 -xn) ϕ(x n+1 ) x n = x . (29)
For our purposes, h k will be approximated by recursive applications of Q k ∆t , based on the following assumption.

Assumption 2. There exist a time step ∆t * > 0 and p > 0 such that, for 0 < ∆t ≤ ∆t * ,

Q k ∆t h k,∆t = e ∆tλ ∆t (k) h k,∆t , (30)
where

h k,∆t = h k + O(∆t p ), λ ∆t (k) = λ(k) + O(∆t p ). ( 31
)
This assumption means that the time-discretized operator Q k ∆t admits h k as an approximate eigenvector with approximate eigenvalue λ(k). This applies, for example, when X is compact. In this case, precise estimates for the errors in [START_REF] Borkar | Learning algorithms for risk-sensitive control[END_REF] are obtained for g = 0 in [START_REF] Ferré | Error estimates on ergodic properties of discretized Feynman-Kac semigroups[END_REF] and can be extended to g = 0.

In the following, we will drop the subscript ∆t on h k,∆t and λ ∆t (k) to simplify the notations, and will present the algorithm essentially as if h k were an exact eigenvector of Q k ∆t with exact eigenvalue λ(k). However, we should keep in mind that this is only approximately true due to the errors in ∆t. We will comment on this in Sec. IV with specific numerical examples.

B. Stochastic approximation and annealing

The main ingredient of the algorithm is the limit (26) of the Feynman-Kac semi-group, which shows that h k and λ(k) can be computed by successively applying Q k ∆t to an initial guess ϕ, so as to obtain

(Q k ∆t ) n ϕ ∼ e n∆tλ(k) h k X ϕ(x) l k (x) dx (32)
as n → ∞. To perform this iteration, which is a functional version of the well-known power method for matrices [START_REF] Demmel | Applied Numerical Linear Algebra[END_REF], we apply a stochastic approximation [START_REF] Basu | A learning algorithm for risk-sensitive cost[END_REF][START_REF] Benveniste | Adaptive Algorithms and Stochastic Approximations[END_REF][START_REF] Polyak | Acceleration of stochastic approximation by averaging[END_REF] whereby the expectation appearing in the action of Q k ∆t is replaced by the iterates of the Markov chain:

Q k ∆t ϕ(x n ) = E e k f (xn)∆t+g x n+1 +xn 2 •(x n+1 -xn) ϕ(x n+1 ) x n ≈ e k f (xn)∆t+g x n+1 +xn 2 •(x n+1 -xn) ϕ(x n+1 ), (33)
where x n+1 is a random variable distributed according to Q ∆t (x n , •). This approximation is known to reproduce the expectation as a statistical average in the ergodic limit n → ∞ [START_REF] Ahamed | Adaptive importance sampling technique for Markov chains using stochastic approximation[END_REF][START_REF] Basu | A learning algorithm for risk-sensitive cost[END_REF][START_REF] Benveniste | Adaptive Algorithms and Stochastic Approximations[END_REF].

In our case, we simulate not the Markov chain x n , but a modified chain, corresponding to the discretization of the driven process ( 16), which we express as

xn+1 = xn + [b(x n ) + F n (x n )]∆t + σ √ ∆t ξ n , (34)
where

F n = D(kg + ∇ log h n k ) ( 35
)
is the extra biasing force derived, according to [START_REF] Aristoff | Adaptive multilevel splitting in molecular dynamics simulations[END_REF], from the estimate h n k of h k at time n. In this case, the evolution ( 33) is modified by Girsanov formula [START_REF] Pavliotis | Stochastic Processes and Applications[END_REF] to

Q k ∆t ϕ(x n ) = E e k f (xn)∆t+g xn+1 +xn 2 •(x n+1 -xn) ϕ(x n+1 )R n (x n , xn+1 ) xn = x ≈ e k f (xn)∆t+g xn+1 +xn 2 •(x n+1 -xn) ϕ(x n+1 )R n (x n , xn+1 ), (36)
where

R n (x n , xn+1 ) = exp - 1 2σ 2 F 2 n (x n )∆t - ∆t σ 2 F n (x n ) • ξ n ( 37
)
is the Radon-Nikodym derivative of the transition kernel of x n with respect to that of xn .

In the end, we also apply an annealing scheme, commonly used in stochastic approximations, which consists in replacing the update rule [START_REF] Chetrite | Variational and optimal control representations of conditioned and driven processes[END_REF], defined by ϕ 0 = ϕ and ϕ n+1 = Q k ∆t ϕ n , by the scheme

ϕ n+1 = ϕ n + a n Q k ∆t ϕ n -ϕ n , (38)
where (a n ) n∈N is a decreasing sequence, often called the adaptation or learning sequence, which acts as a smoothing parameter, filtering here the noisy update of the eigenfunction. This sequence is usually chosen in such a way that n≥0

a n = ∞, n≥0 a 2 n < ∞, (39)
with the understanding that a n should not be decreased too slowly, so as to limit noise, nor too fast, so as to reach the "correct" fixed point. These conditions can be relaxed under stability assumptions or by an averaging procedure [START_REF] Polyak | Acceleration of stochastic approximation by averaging[END_REF]; see [START_REF] Benveniste | Adaptive Algorithms and Stochastic Approximations[END_REF] for more details.

C. Spatial projection

The iteration just described for approximating h k can be performed numerically by discretizing the state space X into small cells (grid discretization). For high-dimensional systems, however, it is more convenient to use a Galerkin-type approximation of the eigenvalue and eigenfunction [START_REF] Basu | A learning algorithm for risk-sensitive cost[END_REF], obtained by projecting the problem onto a set of basis functions {φ j } M j=1 with φ j : X → R [START_REF] Chatelin | Spectral Approximation of Linear Operators[END_REF]. Let us denote by H M = Span{φ j } M j=1 the space spanned by these functions. Then the M -dimensional eigenproblem that we need to solve is

Q k ∆t h = e ∆tλ h, h ∈ H M , (40)
using, for notational convenience, the same symbols h and λ for the exact and the projected spectral elements. We also drop from now on the parameter k, which will be implicit.

The eigenfunction h is expressed in that basis as

h(x) = M j=1 α j φ j (x) = α T φ(x), (41)
where α = [α 1 , . . . , α M ] T and φ(x) = [φ 1 (x), . . . , φ M (x)] T . Multiplying (40) by φ i for i ∈ {1, . . . , M } and integrating over any measure η on X yields

M j=1 α j X φ i (Q k ∆t φ j ) dη = e ∆tλ M j=1 α j X φ i φ j dη, i ∈ {1, . . . , M }. (42)
As a result, we see that the vector of coefficients α ∈ R M is the principal solution of the eigenproblem

Aα = ΛBα, λ = 1 ∆t log Λ, (43)
where

A = X φ(Q k ∆t φ T ) dη, B = X φ φ T dη. ( 44
)
Note that the matrix B is invertible as soon as the {φ j } M j=1 form a linearly independent family.

D. Algorithm

We are now ready to describe all the steps of the algorithm that estimates the principal eigenvalue of L k and its corresponding eigenfunction. For a fixed k ∈ R, we initiate the process at a position x 0 ∈ X and define a first approximation h 0 of h as

h 0 = (α 0 ) T φ, ( 45
)
where φ is the vector of basis functions and α 0 is the initial vector of coefficients, chosen such that h 0 is constant. At each iteration, we then perform the following steps:

1. Draw a new position xn+1 according to the Markov chain [START_REF] Roberts | General state space Markov chains and MCMC algorithms[END_REF].

2. Compute the extra bias F n according to [START_REF] Benveniste | Adaptive Algorithms and Stochastic Approximations[END_REF], which in the function basis takes the form

F n = D kg + M j=1 α n j ∇φ j M j=1 α n j φ j . ( 46
)
3. Compute the Girsanov weight R n according to [START_REF] Pavliotis | Stochastic Processes and Applications[END_REF].

4. Compute the matrices A n+1 and B n+1 using the formulae

A n+1 = 1 n + 1 n m=0 e k f (xm)∆t+g xm+1 +xm 2 (x m+1 -xm) φ(x m )φ(x m+1 ) T R m (x m , xm+1 ), B n+1 = 1 n + 1 n m=0 φ(x m )φ(x m ) T , (47)
which follow by projecting ( 36) and ( 44), respectively.

5. Update the coefficient vector α n , giving the decomposition of the iterate h n , as

α n+1 = α n + a n B -1 n A n h n (x 0 ) -Id α n , ( 48
)
where Id is the identity matrix of size M [START_REF] Basu | A learning algorithm for risk-sensitive cost[END_REF].

6. Estimate the eigenvalue as

λ n+1 = n n + 1 λ n + 1 (n + 1)∆t log (h n (x 0 )) , (49)
which follows from [START_REF] Gorissen | Finite size scaling of current fluctuations in the totally asymmetric exclusion process[END_REF].

Repeating these steps, it can be proved that the iterates h n and λ n converge to the solution of the spectral problem [START_REF] Bierkens | Linear PDEs and eigenvalue problems corresponding to ergodic stochastic optimization problems on compact manifolds[END_REF], following the analysis found in [START_REF] Borkar | Peformance analysis conditioned on rare events: An adaptive simulation scheme[END_REF][START_REF] Ahamed | Adaptive importance sampling technique for Markov chains using stochastic approximation[END_REF][START_REF] Basu | A learning algorithm for risk-sensitive cost[END_REF]; see also [START_REF] Benveniste | Adaptive Algorithms and Stochastic Approximations[END_REF]. Moreover, because h n → h, the Markov chain (x n ) n∈N samples in the long run the Euler-Maruyama discretization of the driven process [START_REF] Cérou | Adaptive multilevel splitting for rare event analysis[END_REF].

E. Remarks

The following are technical remarks worth noting about the algorithm:

1. The matrix A n can be updated at each step using

A n+1 = n n + 1 A n + 1 n + 1 e k f (xn)∆t+g xn+1 +xn 2 (x n+1 -xn) φ(x n )φ(x n+1 ) T R n , (50)
instead of the sum shown in [START_REF] Nemoto | Finite-time and finite-size scalings in the evaluation of large-deviation functions: Analytical study using a birth-death process[END_REF]. Similarly, the matrix B n can be updated, following [START_REF] Basu | A learning algorithm for risk-sensitive cost[END_REF], using the Shermann-Morrisson-Woodburry formula, which leads here to

B -1 n+1 = n + 1 n B -1 n - n + 1 n B -1 n φ(x n+1 )φ(x n+1 ) T B -1 n R n+1 n n+1 + φ(x n+1 ) T B -1 n φ(x n+1 )R n+1 . ( 51
)
The advantage of this formula is that the computation required for updating the coefficient α n scales with the number M of basis vectors as M 2 , as is common in the power method [START_REF] Demmel | Applied Numerical Linear Algebra[END_REF], whereas the typical cost of inverting the matrix B scales as M 3 .

2. We normalize the eigenfunction at every iteration by setting h(x 0 ) = 1 at an arbitrary location, taken here to be the initial state x0 [START_REF] Borkar | Peformance analysis conditioned on rare events: An adaptive simulation scheme[END_REF]. This prevents the norm of the eigenfunction from diverging or decaying to zero, as is common in the power method, and provides an estimate of λ(k) through [START_REF] Demmel | Applied Numerical Linear Algebra[END_REF]. This normalization step can be based on other norms, at the expense of computing integrals.

3. The algorithm can be run with the original (unbiased) process X t , but the estimation of the matrices A and B in this case typically suffers from large statistical errors due to the high variance of the underlying estimator. This is a known problem related to the estimation of Feynman-Kac functionals and exponential integrals in general [START_REF] Hartmann | Efficient rare event simulation by optimal nonequilibrium forcing[END_REF][START_REF] Hartmann | Characterization of rare events in molecular dynamics[END_REF][START_REF] Zhang | Applications of the cross-entropy method to importance sampling and optimal control of diffusions[END_REF][START_REF] Rohwer | Convergence of large deviation estimators[END_REF]. Biasing the dynamics with the driven process Xt reduces this variance in an optimal way (in the sense of asymptotic or logarithmic efficiency) by forcing the exploration of the process in important regions of the state space where the integrand of the generating function E[e kT A T ] is largest [START_REF] Chetrite | Variational and optimal control representations of conditioned and driven processes[END_REF].

4. As the dynamics is biased towards the driven process Xt , the vector α n of basis coefficients representing h k converges towards the solution of the eigenproblem [START_REF] Gorissen | Finite size scaling of current fluctuations in the totally asymmetric exclusion process[END_REF], with the matrices A and B computed as ergodic averages under Xt , so that η = µ k in [START_REF] Fleming | Controlled Markov Processes and Viscosity Solutions[END_REF]. Even for a small number of basis functions, it would be a priori impossible to compute these matrices by numerical quadrature. This shows that the algorithm can be used to obtain good approximations of h k even for high-dimensional systems, provided that enough basis functions are used to represent the support of µ k , which is typically concentrated on a subset of X .

5. The algorithm is stable despite the fact that it includes the Girsanov reweighting factor, which is exponential in time. The reason for this stability, already noted in [START_REF] Ahamed | Adaptive importance sampling technique for Markov chains using stochastic approximation[END_REF], is that the Girsanov weight is computed and accumulated incrementally over single time steps.

6. The learning sequence (a n ) n∈N is chosen here in the following way:

(a) For 0 ≤ n ≤ N 1 , we take a n = 0, so there is no adaption at the level of h and λ, although the matrices A and B are evolved. This "burn-in" period allows for a better initial guess of the various functions estimated through their ergodic averages.

(b) For N 1 < n ≤ N 1 + N iter -N 2 , we take a n = 1, that is, the full information of the process is taken into account.

(c) For

N 1 + N iter -N 2 < n ≤ N 1 + N iter , we take a n = C/(n -(N 1 + N iter -N 2 )
) with constant C > 0, so the process learns less with time, smoothing the noise in the long run. In the following, we choose C = 1.

The times N 1 and N 2 can be fixed or can be chosen dynamically according to some stopping rule.

7. In practice, we can perform independent simulations with different values of k to obtain an interpolation of the SCGF over some range, say, [k min , k max ], which can then be used to obtain the rate function by Legendre transform. Alternatively, we can do a simulation in which k is slowly increased from k = 0 in a "quasi-static" way, so as to adaptively update the biasing force over a range of values for k [START_REF] Nemoto | Computation of large deviation statistics via iterative measurement-andfeedback procedure[END_REF] or to reach some prescribed value of the SCGF for an unknown k [START_REF] Borkar | Peformance analysis conditioned on rare events: An adaptive simulation scheme[END_REF].

IV. APPLICATIONS

We apply in this section the algorithm to two simple test cases involving one-dimensional diffusions. The first is the Ornstein-Uhlenbeck process, for which the large deviation functions of the area per unit time are known exactly [START_REF] Chetrite | Nonequilibrium Markov processes conditioned on large deviations[END_REF], while the second is a driven diffusion on the circle, often used in physics as a model of nonequilibrium systems, including Josephson junctions perturbed by thermal noise and Brownian particles controlled by external forces [START_REF] Risken | The Fokker-Planck Equation: Methods of Solution and Applications[END_REF][START_REF] Reimann | Brownian motors: Noisy transport far from equilibrium[END_REF][START_REF] Ciliberto | Fluctuations in out-of-equilibrium systems: From theory to experiment[END_REF]. We discuss for both the convergence and efficiency of the algorithm.

A. Ornstein-Uhlenbeck process

The first example that we consider is the mean area or mean position

A T = 1 T T 0 X t dt (52)
of the Ornstein-Uhlenbeck process on R satisfying the SDE

dX t = -2θX t dt + √ 2 dW t , (53)
where θ > 0. In the notations of Sec. II, we thus have b

(x) = -2θx, σ = √ 2, f (x) = x and g(x) = 0, so that L = -2θx d dx + d 2 dx 2 , (54)
and

L k = L + kx.
For this process and observable, it can be checked [START_REF] Chetrite | Nonequilibrium Markov processes conditioned on large deviations[END_REF] that the SCGF, corresponding to the dominant eigenvalue of L k , is

λ(k) = k 2 4θ 2 , ( 55
)
so that I(a) = θ 2 a 2 from (25). This is expected, since the integral of a Gaussian process is also Gaussian. Moreover, the associated dominant eigenfunction is

h k (x) = e k 2θ x , (56)
leading with [START_REF] Aristoff | Adaptive multilevel splitting in molecular dynamics simulations[END_REF] to the optimal drift

b k (x) = -2θx + 2(log h k (x)) = -2θx + k θ . ( 57
)
This shows that a fluctuation of A T is created in an optimal way by adding a constant to the drift, which "moves" the Gaussian stationary density of the Ornstein-Uhlenbeck process

ρ * (x) = θ π e -θx 2 (58)
to leading to

ρ k (x) = θ π e -θ(x-m k ) 2 , m k = k 2θ 2 , (59)
A T T →∞ ----→ m k = a k ( 60
)
almost surely with respect to Xt , in agreement with [START_REF] Lecomte | A numerical approach to large deviations in continuous time[END_REF]. We compare the algorithm against these exact results using a simple mesh discretization (firstorder finite elements), defined by the basis functions

φ j (x) =    x-x j δ + 1, x ∈ [x j -δ, x j ], - x-x j δ + 1, x ∈ [x j , x j + δ], 0, otherwise, (61)
where the points x j define the centers of each "cell" of width δ > 0. In the following, we refer to this basis simply as the linear basis. To illustrate the flexibility of the algorithm, we also perform simulations using a Hermite polynomial basis, which forms a complete orthonormal basis in L 2 (ρ * ). We run the algorithm with k = 1, θ = 1/2, x 0 = 0, T = 2 × 10 4 , T 1 = N 1 ∆t = T 2 = N 2 ∆t = 2 × 10 3 with ∆t = 5 × 10 -3 , which a standard time step used in simulations relative to the basic timescale of the dynamics, corresponding here to 1/(2θ) = 1. For the linear basis, we use M = 61 equally spaced cells with δ = 0.25 around x = 0, whereas for the Hermite basis we use only M = 10 basis functions. Figure 1 illustrates the results of a typical simulation, starting on the left with the evolution of the estimate of the SCGF, given by ( 49), as time increases. We observe a very good agreement in the long run with the exact value, which for the parameters used is equal to 1, with a faster convergence observed for the Hermite basis compared to the linear basis. As mentioned before, we can also recover the SCGF by recording the stationary value of A T , which corresponds as above to m k = a k = λ (k), and numerically integrate the result in k from λ(0) = 0. The results obtained are similar to those obtained from the eigenvalue estimate (49), and are not shown for this reason.

To understand the convergence of the eigenvalue at the process level, we show in the right plot of Fig. 1 the evolution of the effective potential V k (x) associated with the modified drift b k (x) according to We show the results for both the linear basis and the Hermite basis, with the zero of the potential arbitrarily set at V k (0) = 0. We see from these that the tail of the potential takes longer to be estimated correctly, as the process starts to explore values away from m 0 = 0. After the convergence time, there is a good agreement with the parabola

b k (x) = -V k (x). ( 62
)
V k (x) = θ(x -m k ) 2 , (63)
which is the exact solution predicted by [START_REF] Tsobgni Nyawo | Large deviations of the current for driven periodic diffusions[END_REF]. The small errors are due to the finite time step ∆t used and the finite basis function set. The small "wiggles" seen in the potential obtained with the linear basis come from the fact that this basis is piecewise linear, so that V k (x) is piecewise logarithmic. The potential obtained with the Hermite basis is smoother, as expected.

Repeating the simulations for other values of k, we can recover λ(k), as shown in the left plot of Fig. 2. The result is in good agreement with the exact solution [START_REF] Nemoto | Computation of large deviation statistics via iterative measurement-andfeedback procedure[END_REF]. Estimating λ (k) with [START_REF] Lecomte | A numerical approach to large deviations in continuous time[END_REF] leads to an estimated rate function, by the Legendre transform [START_REF] Heymann | Pathways of maximum likelihood for rare events in nonequilibrium systems: Application to nucleation in the presence of shear[END_REF], which also agrees well with the exact rate function, shown in the right plot of Fig. 2. By comparison, the rate function obtained from the time average [START_REF] Bolhuis | Transition path sampling: Throwing ropes over rough mountain passes, in the dark[END_REF] of the extra biasing force (with the Hermite basis) is not as good: it lies above the exact rate function and shows a larger offset or error as k increases, which does not decrease by reducing ∆t or increasing the number of basis functions.

To understand this error, we show in the left plot of Fig. 3 the evolution of the extra biasing force F n for k = 1, estimated at the current location of the system, which should approach the constant 2, following [START_REF] Reimann | Brownian motors: Noisy transport far from equilibrium[END_REF]. The evolution is noisy, as can be seen, which is expected, since F n is estimated by the logarithmic derivative

F n (x) = 2 (h n (x)) h n (x) , (64)
computed in the function basis from [START_REF] Nemoto | Population-dynamics method with a multicanonical feedback control[END_REF]. The derivative amplifies the Monte Carlo errors inherent in the estimate h n . In addition, the denominator often takes small values away from the mean position m k , which makes the estimation of the optimal force still more difficult. The right plot of Fig. 3 shows that the noise on F n is considerably filtered out by the time average underlying the estimator R T of the rate function, although a bias remains even after the convergence time of the SCGF, which leads to the offset seen in Fig. 2. The results in both cases are more noisy for the linear basis because the derivative is not continuous across the different cells. The offset on R T remains more or less constant by running independent simulations (error bars on 30 simulations are too small to show), and so appears to be a systematic error or bias. Many factors can account for this bias. First, the rate function estimated from the limit shown in ( 24) is known to be an upper bound on the true rate function [START_REF] Chetrite | Variational and optimal control representations of conditioned and driven processes[END_REF], which is tight if and only if the modified drift estimated in the simulation is the optimal control drift [START_REF] Aristoff | Adaptive multilevel splitting in molecular dynamics simulations[END_REF]. Here, F n is not constant, as predicted from ( 57), so we expect the estimate R T to lie above its expected value. Second, we have noticed in simulations that A T underestimates a k for large k, which has the effect of further "pushing" the estimate of the rate function above I(a). This is most likely due again to F n being non-constant. Finally, the time average of F n involves the ratio of two functions, according to [START_REF] Blanchet | State-dependent importance sampling for rare-event simulation: An overview and recent advances[END_REF], which fluctuate in time via the updating of α n in [START_REF] Nemoto | Population-dynamics method with a multicanonical feedback control[END_REF]. As a result, we expect this additional randomness to artificially increase the second moment of F n , leading to a further bias.

It is difficult to isolate these factors, and all, in fact, seem to play a role. In future works, it would be interesting to study the bias observed in R T by computing, for example, the time average of F 2 n using the final estimate of h k rather than the time-evolved estimate h n . Different annealing sequences or averaging techniques could also be used to "filter" the time average of F 2 n and mitigate the bias on R T . For now, the most efficient and reliable way to compute the rate function is to use the estimator based on the SCGF and its Legendre transform.

B. Periodic diffusion

We consider for the second test a diffusion on the unit circle satisfying the SDE

dX t = (-V (X t ) + γ)dt + √ 2 dW t , (65)
with the periodic potential V (x) = cos(2πx) and γ ∈ R a constant drive. For γ = 0, the total drift cannot be expressed as the gradient of a smooth periodic function, so X t is a nonequilibrium process violating detailed balance [START_REF] Risken | The Fokker-Planck Equation: Methods of Solution and Applications[END_REF]. The observable studied for this process is the winding number calculated with the real rather than periodic state, which can be interpreted physically as the mean velocity or current of a Brownian particle moving around the circle [START_REF] Tsobgni Nyawo | Large deviations of the current for driven periodic diffusions[END_REF]. In the notation of Sec. II, we have b

A T = 1 T T 0 dX t , (66)
(x) = -V (x) + γ, σ = √ 2, f (x)
= 0, and g(x) = 1. The SCGF and the rate function of this observable are not known exactly, but Galerkin approximations can easily be found by Fourier series, using the basis functions φ j (x) = e j2πix with j ∈ Z [START_REF] Tsobgni Nyawo | Large deviations of the current for driven periodic diffusions[END_REF][START_REF] Mehl | Large deviation function for entropy production in driven onedimensional systems[END_REF][START_REF] Nemoto | Variational formula for experimental determination of high-order correlations of current fluctuations in driven systems[END_REF]. This basis is used to compute reference values for the SCGF and the rate function, devoid of ∆t errors, by projecting the spectral problem [START_REF] Aminov | Singularities in large deviation functionals of bulk-driven transport models[END_REF] in Fourier space and by ensuring that enough basis functions are used. We find in practice that M = 41 Fourier modes are sufficient. We use the same Fourier basis for the algorithm, also with M = 41 modes, in addition to T = 2 × 10 4 and T 1 = T 2 = 2 × 10 3 for the integration times, as in the first test.

We show in Fig. 4 the results of a typical simulation for γ = 1 and k = 1. The left plot in this figure shows the evolution in time of the estimated SCGF from the eigenvalue iteration (49), while the right plot shows the final estimated eigenfunction h k for different time steps ∆t, compared with the reference Fourier solution. We see that the results for the eigenvalue and the eigenvector significantly depart from the reference solutions for ∆t = 0.05, but converge to them as ∆t is decreased, in accordance with Assumption 2 and the theoretical results of [START_REF] Ferré | Error estimates on ergodic properties of discretized Feynman-Kac semigroups[END_REF]. For the same time step used for the Ornstein-Uhlenbeck process, namely, ∆t = 5 × 10 -3 , we see no notable difference between the estimated eigenfunction and the reference values, leading to a precise estimation of the SCGF. The convergence of the eigenvalue here is much faster than for the Ornstein-Uhlenbeck process because the space explored is compact, being limited to [0, 1] with periodic boundary conditions.

As for the Ornstein-Uhlenbeck process, we can repeat these simulations over a range of values for k to obtain the SCGF and the rate function. The left plot of Fig. 5 shows that the SCGF is in good agreement with the Fourier solution for ∆t = 5 × 10 -3 , and so is the rate function estimated by Legendre transform. However, as seen before, the rate function estimated with the time average R T of the biasing force shows an offset, although smaller this time, which comes from the noisy estimation of F n . As before, the estimator of the rate function that should be used is the one based on the Legendre transform of the SCGF, with a k estimated by A T .

V. CONCLUSION

We have presented a new algorithm for estimating the large deviation functions of time-integrated observables of Markov processes, which characterize the likelihood of their fluctuations in the longtime limit. The algorithm draws on earlier results on stochastic control [START_REF] Borkar | Peformance analysis conditioned on rare events: An adaptive simulation scheme[END_REF][START_REF] Ahamed | Adaptive importance sampling technique for Markov chains using stochastic approximation[END_REF][START_REF] Basu | A learning algorithm for risk-sensitive cost[END_REF][START_REF] Borkar | Learning algorithms for risk-sensitive control[END_REF], and works by adaptively estimating the principal eigenvalue and eigenfunction of a spectral problem related to the large deviation problem. The adaptive part consists in modifying the process considered, using feedback and reinforcement learning, so as to reach the so-called driven process, which is known to be optimal for the purpose of estimating large deviation functions using importance sampling [START_REF] Chetrite | Variational and optimal control representations of conditioned and driven processes[END_REF]. In this sense, the algorithm relates to many adaptive importance sampling methods that have been proposed recently for rare event simulations; see, for instance, [START_REF] Dupuis | Dynamic importance sampling for uniformly recurrent Markov chains[END_REF][START_REF] Nemoto | Finite-size scaling of a first-order dynamical phase transition: Adaptive population dynamics and an effective model[END_REF][START_REF] Blanchet | State-dependent importance sampling for rare-event simulation: An overview and recent advances[END_REF][START_REF] Kappen | Adaptive importance sampling for control and inference[END_REF]. It is also closely related to diffusion Monte Carlo methods [START_REF] Rousset | On the control of an interacting particle estimation of Schrödinger ground states[END_REF][START_REF] Foulkes | Quantum Monte Carlo simulations of solids[END_REF][START_REF] Lim | Fast randomized iteration: Diffusion Monte Carlo through the lens of numerical linear algebra[END_REF], which attempt to estimate the ground-state wavefunction of a many-body quantum system by simulating a related stochastic process.

The proposed algorithm can be applied to diffusions, as illustrated here, but also to Markov chains in discrete or continuous time, and works for both equilibrium and nonequilibrium systems, that is, reversible and non-reversible processes, respectively. Moreover, although the test cases that we have presented are simple, they clearly show that the algorithm has the potential to improve upon other simulation methods, especially cloning methods, since it runs on a single simulation of the process and provides information about how fluctuations are created in time by constructing the driven process in a non-parametric form and with no prior information. A modification of the original cloning algorithm was proposed recently [START_REF] Nemoto | Population-dynamics method with a multicanonical feedback control[END_REF] to construct the driven process, but it is based on a different feedback rule that compares two time-dependent histograms, whose estimation is noisy and requires a large number of clones. The results presented here show that a single clone, evolving over a long-enough time, is sufficient. This obviously cuts the computational complexity of estimating large deviations, but also simplifies, as mentioned, the error and convergence analyses of the algorithm.

As with any new proposal, more work is needed to understand the benefits and limitations of the algorithm, to test its applicability to realistic systems, and to benchmark it against other numerical methods. Of particular importance is to derive precise error estimates associated with the space and time discretizations of the spectral elements, large deviation functions, and the driven process.

Recent results about the discretization errors in ∆t for the SCGF can be found in [START_REF] Ferré | Error estimates on ergodic properties of discretized Feynman-Kac semigroups[END_REF]. These discretization errors are also present in cloning algorithms when applied to diffusions, so they are not specific to the algorithm presented here. The use of Galerkin discretizations is also common to both algorithms and requires further investigations, particularly in the low-noise limit and for processes involving many interacting particles.

FIG. 1 .

 1 FIG.1. Large deviations of the mean area for the Ornstein-Uhlenbeck process. Left: Evolution of the estimated SCGF in time. Right: Evolution of the biasing potential V k (x) in time for the linear basis (curves from top to bottom), compared with the final estimate for the Hermite basis and the exact result.

FIG. 2 .

 2 FIG.2. Large deviations of the mean area for the Ornstein-Uhlenbeck process. Left: Estimated SCGF compared with the exact result. Right: Estimated rate function based on the Legendre transform of the SCGF and on the averaging of the biasing force (using the Hermite basis), compared with the exact result.

FIG. 3 .

 3 FIG.3. Large deviations of the mean area of the Ornstein-Uhlenbeck process. Left: Biasing force F n at the current state xn of the process at time n. Right: Evolution of the time average R T , which gives the rate function according to[START_REF] Bolhuis | Transition path sampling: Throwing ropes over rough mountain passes, in the dark[END_REF].

FIG. 4 .

 4 FIG. 4. Current large deviations of the driven periodic diffusion. Left: Evolution of the estimated SCGF with time for different time steps ∆t. Right: Estimated eigenfunction h k for different time steps ∆t compared with the reference Fourier solution.

FIG. 5 .

 5 FIG. 5. Current large deviations for the driven periodic diffusion. Left: Estimated SCGF compared with the reference Fourier solution. Right: Estimated rate function based on the Legendre transform and on the averaging of the biasing force, compared with the reference Fourier solution.

See[START_REF] Chetrite | Nonequilibrium Markov processes conditioned on large deviations[END_REF] for a treatment of diffusions with multiplicative noise.

ACKNOWLEDGMENTS

We are grateful to Florian Angeletti for useful discussions in the initial phase of this work. G.F. is supported by the Labex Bezout. H.T. was supported by the National Research Foundation of South Africa (Grant nos 90322 and 96199) and Stellenbosch University (Project Funding for New Appointee). This research was also supported in part by the International Centre for Theoretical Sciences (ICTS) during a visit for participating in the program "Large deviation theory in statistical physics: Recent advances and future challenges" (Code: ICTS/Prog-ldt/2017/8).


Cle ´ment Feger

Laurent Mermet

Bhaskar Vira

Prue F E Addison

Richard Barker

Frank Birkin

John Burns

Stuart Cooper

Denis Couvet

Thomas Cuckston

Gretchen C Daily

Colin Dey

Louise Gallagher

Rosemary Hails

Stephen Jollands

Georgina Mace

Emily Mckenzie

Markus Milne

Paolo Quattrone

Alexandre Rambaud

Shona Russell

Marta Santamaria

William J Sutherland

Four priorities for new links between conservation science and accounting research

Introduction

Engagement with diverse social science disciplines is essential to revealing political, social, and institutional challenges that must be addressed to advance effective biodiversity conservation [START_REF] Bennett | Mainstreaming the social sciences in conservation[END_REF][START_REF] Teel | Publishing social science research in Conservation Biology to move beyond biology[END_REF]. One challenge that remains insufficiently investigated is frustration with the lack of impact of innovative information tools and systems of accounts aimed at motivating and guiding ecosystem management. The conservation community invests considerable efforts in their creation and experimentation. Species and ecosystem accounts (e.g., ABoS 2015; UNEP-WCMC 2016), general ecological indicators (e.g., [START_REF] Jørgensen | Twenty volumes of ecological indicators -an accounting short review[END_REF], and tools for ecosystem-services quantification and mapping (e.g., [START_REF] Kareiva | Natural capital: theory and practice of mapping ecosystem services[END_REF]) and ecosystem monitoring are fundamental to conservation research and practice. However, ecosystem-based tools do not always lead to the changes in decision, action, or policy conservation scientists expect (e.g., [START_REF] Ruckelshaus | Notes from the field: lessons learned from using ecosystem service approaches to inform real-world decisions[END_REF].

Often, the inability of such information systems to generate expected changes is not due to technical limitations rather than the too fragile articulation between their design and the complex realities of developing strategies and organizing management of ecosystems in a diversity of contexts. Investigating such articulation between an information system and the organizational details of its systematic use is precisely what characterizes an academic field: accounting, which belongs to management as a discipline and often intersects with social sciences or economics. Accounting has enormous but untapped potential to contribute to conservation science, practice, and goals. Accounting is often misconceived as being only the craft of producing quantitative and financially focused reports for companies. However, accounting in its broadest sense is the preparation and the framing of information (qualitative and quantitative) to assist specific organizing and decision-making processes (Jollands 2017).

We especially refer here to critical and interpretive accounting research, a field that emerged in the 1970s and developed through now well-established journals (i.e., Accounting, Organizations and Society, Accounting Auditing and Accountability Journal, Critical Perspectives on Accounting) [START_REF] Miller | Accounting, organizing, and economizing: connecting accounting research and organization theory[END_REF][START_REF] Roslender | The Routledge companion to critical accounting[END_REF]). Since 1990s, researchers have revealed and criticized the lack of consideration of sustainability issues in existing accounting systems (e.g., [START_REF] Milne | On sustainability, the environment and management accounting[END_REF] and advocated the development of new accounting approaches inspired by ecological thinking at and beyond the corporate level (e.g., [START_REF] Birkin | The ecological accountant: from the cogito to thinking like a mountain[END_REF][START_REF] Bebbington | Accounting and sustainable devel-opment: an exploration[END_REF][START_REF] Russell | Accounts of nature and the nature of accounts: critical reflections on environmental accounting and propositions for ecologically informed accounting[END_REF].

Following a recent publication that proposes a new line of inquiry focusing on developing accounting research at the ecosystem management level [START_REF] Feger | A blueprint towards accounting for the management of ecosystems[END_REF], a workshop furthered in-depth interdisciplinary dialogue between accounting scholars and conservation researchers and practitioners. Its results underline that collaboration between conservation and accounting research is essential to improve the design and the actual use of ecosystem-based information systems for accountable conservation decisions and actions. Four key areas for future joint research were identified.

What Accounting Brings to Conservation

Our call to establish new links between the accounting discipline and biodiversity conservation is not meant to be a substitute for economics, game theory, organizational psychology, or any other discipline focusing on decision making. It is an invitation to focus on questions instrumental and common to both conservation and accounting research, such as the following: How are records kept in practice and with what consequences? What languages and representations can one provide to complex organizations? Who gives and demands what kind of accounts? How are responsibilities negotiated, organized, managed, and controlled? How are explicit principles and conventions on which accounts can be developed and values defined and on which past and future actions can be assessed and compared debated and institutionalized?

The pervasive confusion in the environmental field between the disciplines of accounting and economics deserves a special comment. Although economics and accounting are somewhat related, they are essentially different disciplines [START_REF] Shiozawa | Economics and accounting: a comparison between philosophical backgrounds of the two disciplines in view of complexity theory[END_REF]. Accounting is concerned with developing and using calculative practices to support decision making as is economics. The use of economics in conservation science has brought major results, considering, for instance, the development of economic valuation of ecosystem services, analysis of environmental trade-offs, and study of incentive structures [START_REF] Helm | The economic analysis of biodiversity: an assessment[END_REF]. One of the distinctive characteristics of accounting, however, is that it focuses on the detailed analysis of the roles of information systems in the context of the concrete complexities of organizational management based on the fundamental concepts of accounts and accountability [START_REF] Burchell | The roles of accounting in organizations and society[END_REF][START_REF] Roberts | Accounting systems and systems of accountability -understanding accounting practices in their organisational contexts[END_REF][START_REF] Gray | New accounts: towards a reframing of social accounting[END_REF]. In terms of methods, accounting research combines theoretical developments that extensively draw on other social science disciplines (organizational theory, sociology, philosophy, economics, psychology, etc.) with in-depth qualitative field studies of organizations [START_REF] Ahrens | Doing qualitative field research in man-agement accounting: positioning data to contribute to theory[END_REF]. In doing so, it enriches understanding of the role of information systems and accounts in the operationalization of action and generation of intended or unintended organizational changes and wider governance transformations [START_REF] Miller | Governing by numbers: why calculative practices matter[END_REF][START_REF] Macintosh | Management accounting and control systems. an organizational and sociological approach[END_REF].

The new dialogue we advocate between conservation scientists and accounting researchers can build on a small but growing body of work in accounting research, centered on ecosystems, that aims to study the effects of varying forms of accounting on relations between human organizations and biodiversity (e.g., [START_REF] Tregidga | Biodiversity offsetting: problematisation of an emerging governance regime[END_REF][START_REF] Dey | Who speaks for the river? Exploring biodiversity accounting using an arena approach[END_REF][START_REF] Cuckston | Ecology-centred accounting for biodiversity in the production of a blanket bog[END_REF]) and develop accounting innovations adapted to the collective management of ecosystems [START_REF] Feger | A blueprint towards accounting for the management of ecosystems[END_REF].

Priorities for Development of Accounting for Ecosystem Management Studying Ecosystem-Centered Accountabilities

A first priority is to study in depth how, in diverse ecosystem management situations, stakeholders actually use or could use ecological and related social, health, economic, and financial information to assign respon-sibilities to one another and to discuss, negotiate, and manage reciprocal commitments (i.e., accountabilities) for improving environmental outcomes. This means ex-ploring questions such as what commitments have been, are being, or should be negotiated among stakeholders; who is accountable to whom and who is not regarding management of ecosystem quality; and how should in-formation be framed and exchanged to organize these accountabilities effectively? An accounting lens can illu-minate how different ways of structuring, representing, giving, and demanding environmental information can lead to creation of viable forms of ecosystem-centered management to achieve conservation goals [START_REF] Roberts | Accounting systems and systems of accountability -understanding accounting practices in their organisational contexts[END_REF][START_REF] Dey | Who speaks for the river? Exploring biodiversity accounting using an arena approach[END_REF][START_REF] Cuckston | Ecology-centred accounting for biodiversity in the production of a blanket bog[END_REF][START_REF] Feger | A blueprint towards accounting for the management of ecosystems[END_REF].

Working Collaboratively on Real-World Cases

Conservation scientists and accounting researchers need to jointly conduct in-depth studies and comparisons of realworld field cases through an accounting lens. Thus, a portfolio of case studies reflecting on past cases and observing and documenting active on-going cases (e.g., through action research interventions) need development.

Adopting a Constructive, Practical, Critical, and Reflective Approach

In working collaboratively on concrete cases, conserva-tion scientists, accounting researchers, and decision mak-ers will engage in constructive discussion to improve the design and use of ecosystem-based information tools. This call for pragmatic trial-and-error approaches that rely on action-oriented agendas and reflexive cultures that are common to conservation science (e.g., adaptive man-agement [START_REF] Gunderson | Panarchy: understanding transformations in human and natural systems[END_REF] or evidence-based conservation [START_REF] Sutherland | The need for evidence-based conservation[END_REF]) and accounting research [START_REF] Gray | The social accounting project and accounting organizations and society. Privileging engagement, imaginings, new accountings and pragmatism over critique?[END_REF].

Developing a Common Language

These 4 priority goals require intensive interdisciplinary dialogue and the development of a common language. Accounting concepts need to be adapted and enriched to analyze and discuss the organizing of ecosystem management and conservation action (e.g., ecological account, accounting entities, accounting perimeters, and accountabilities) [START_REF] Russell | Accounts of nature and the nature of accounts: critical reflections on environmental accounting and propositions for ecologically informed accounting[END_REF]. The specificities of accounting concepts, as distinct from concepts used in the field of economics or ecology, need theoretical clarification, especially when terms overlap (e.g., valu-ation and capital) [START_REF] Rambaud | The 'triple depreciation line' instead of the 'triple bottom line': towards a genuine integrated reporting[END_REF]. Finally, the formulation of new concepts and vocabularies (e.g., reciprocal commitments) has to be central to joint efforts of accountants and conservation scientists to develop an accounting approach for the management of ecosystems.

Conclusion

The 4 priorities for the development of accounting approaches centered on management of ecosystems set up an agenda that can reshape conservation practice and the way ecosystem-based information tools are designed and used in conservation and accounting research and the way accounting entities and accountabilities are understood. By collaboration and engagement across these disciplines, there is scope for contributing to constructive critical reasoning and to introduce inno-vative designs that combine insights from accounting and conservation. Ultimately, this new interdisciplinary bridge will provide a critical, theoretical, and practical addition to the already well-established collaborations of conservation research with other social science fields, such as economics, anthropology, and political ecology.

Acknowledgments

This article is the result of an interdisciplinary dialogue between conservation and accounting researchers that took place during the New Accounting for the Manage-ment of Ecosystems workshop. This event was organized at the University of Cambridge and was coconvened by C.F., B.V., and L.M. to explore opportunities and chal-lenges of developing new accounting approaches cen-tered on the collective management of ecosystems. We are grateful to the Center for Research in the Arts, Social Sciences and Humanities, and the Luc Hoffmann Insti-tute, which cofunded the event and made this dialogue possible. P.F.E.A. is supported by the Natural Environ-ment Research Council number NE/N005457/1. We also thank 2 anonymous reviewers for their insightful and constructive comments that helped improve the paper substantially.


Florian Lemarié
email: florian.lemarie@inria.fr

Hans Burchard

Laurent Debreu

Knut Klingbeil

Jacques Sainte-Marie

Florian Lemari

Advancing dynamical cores of oceanic models across all scales

come   L'archive ouverte pluridisciplinaire

Oceanic numerical models are used to understand and predict a wide range of processes from global paleoclimate scales to short-term prediction in estuaries and shallow coastal areas. One of the overarching challenges, and the main topic of the COMMODORE workshop, is the appropriate design of the dynamical cores given the wide variety of scales of interest and their interactions with atmosphere, sea-ice, biogeochemistry, and even societal processes. The construction of a dynamical core is a very long effort which takes years and decades of research and development and which requires a collaborative mixture of scientific disciplines. This work involves a significant number of fundamental choices, such as which equations to solve, which horizontal and vertical grid arrangement is adequate, which discrete algorithms allows jointly computational efficiency and sufficient accuracy, etc. Nowadays, a broad range of numerical methods are implemented in models used for realistic ocean simulations, and, owed to the advances in computational power, a meeting point has been reached between global circulation models and regional local models such that there can be mutual benefits of a cross-fertilization between communities. This report outlines an initiative to bring together the world-wide leading researchers actively contributing to the development of oceanic model dynamical cores, such that participants could network together and focus on next challenges irrespective of target applications (regional, coastal, or global). The first community for the numerical modeling of the global, regional and coastal ocean (COMMODORE) workshop (https: //commodore2018.sciencesconf.org/) has been organized in Paris in September 2018. In total, the participants represented 15 oceanic dynamical cores among the Submitted to Bull. Amer. Meteor. Soc. 1 most widely used by the research and operational community. The motivations, topics of discussion sessions, and outcomes of the workshop are summarized below.

Context

The ocean model developers community has had the tendency to be split depending on the target applications (global vs coastal) and on the type of horizontal grids (structured vs unstructured) and has been organized around relatively small modeling groups. However, the models have now reached such a high level of complexity that model development goes beyond the expertise of one given group and requires interactions between physicists, mathematicians, and computer scientists. In this context, this workshop aimed at gathering a community of model-oriented researchers to foster more regular exchanges and share expertise on outstanding issues and perspectives. During this first workshop, the emphasis was on reviewing the characteristics and diversity in the formulation of oceanic models used for realistic applications as well as on outlining upcoming challenges.

Evolution of oceanic models

Historically, global and regional ocean models have been based on the hydrostatic primitive equations (e.g. [START_REF] Griffies | Formulating the Equations of ocean models[END_REF] discretized on structured grids using a mixture of finite-difference and finite-volume techniques for the discretization in space. The time dimension is usually treated using standard predictor-corrector or two-level approaches (e.g. Lemarié et al. 2015). Those choices have been made because of their good compromise between simplicity, efficiency, and accuracy. In recent years, significant progress has been made for ocean modeling on unstructured grids, either via finite volume (e.g. [START_REF] Chen | An unstructured grid, finite-volume, three-dimensional, primitive equations ocean model: Application to coastal ocean and estuaries[END_REF][START_REF] Ringler | A unified approach to energy conservation and potential vorticity dynamics for arbitrarily-structured c-grids[END_REF]Danilov et al. 2017), or finite element (e.g. [START_REF] Zhang | Seamless cross-scale modeling with schism[END_REF][START_REF] Korn | Formulation of an unstructured grid model for global ocean dynamics[END_REF][START_REF] Kärnä | Thetis coastal ocean model: discontinuous galerkin discretization for the three-dimensional hydrostatic equations[END_REF] approach. Unstructured grid models have reached an unprecedented level of maturity at least for two reasons. First, the vertical dimension is treated in a structured way compared to earlier initiatives trying to get three-dimensional unstructured meshes working. Second, a better understanding of computational modes and dispersion properties associated to a wide range of possible choice of finite element pairs has been reached (e.g. Le [START_REF] Le Roux | Analysis of numerically induced oscillations in 2d finite-element shallow-water models part i: Inertia-gravity waves[END_REF][START_REF] Le Roux | Spurious inertial oscillations in shallow-water models[END_REF][START_REF] Eldred | Dispersion analysis of compatible galerkin schemes for the 1d shallow water model[END_REF]. For example, [START_REF] Korn | Elementary dispersion analysis of some mimetic discretizations on triangular c-grids[END_REF] have recently proposed a specific mimetic approach to control the well-known spurious mode occurring in triangular C-grids [START_REF] Wolfram | Mitigating horizontal divergence checker-board oscillations on unstructured triangular c-grids for nonlinear hydrostatic and nonhydrostatic flows[END_REF]. Unstructured grid models have been used for coastal applications for many years and they have now reached the application phase for global applications (e.g. [START_REF] Sidorenko | Towards multi-resolution global climate modeling with echam6-fesom. part i: model formulation and mean climate[END_REF][START_REF] Petersen | An evaluation of the ocean and sea ice climate of E3SM using MPAS and interannual CORE-II forcing[END_REF]. A long-standing concern is that computational cost per nominal grid point is generally much larger than for structured-grid models and this problem is further compounded by the absence of time refinement to locally adjust the time-step to the mesh resolution when explicit time-stepping is used. The test strategy presented in the next section should provide a way to quantify more rationally the difference in terms of computational costs among existing models.

Recent advances also include the development of hybrid (or generalized) vertical coordinate systems based on Arbitrary Lagrangian Eulerian (ALE) methods. The vertical distribution of Eulerian coordinate levels is predefined, whereas Lagrangian coordinate levels freely evolve with the flow. ALE methods combine the advantages of welldefined (i.e. undistorded) meshes and reduced numerical mixing, and also allow adaptation strategies (e.g. [START_REF] Bleck | An oceanic general circulation model framed in hybrid isopycnic-Cartesian coordinates[END_REF][START_REF] Burchard | Non-uniform adaptive vertical grids in one-dimensional numerical ocean models[END_REF][START_REF] White | High-order regriddingremapping schemes for continuous isopycnal and generalized coordinates in ocean models[END_REF][START_REF] Leclair | z-coordinate, an arbitrary lagrangianeulerian coordinate separating high and low frequency motions[END_REF][START_REF] Petersen | Evaluation of the arbitrary lagrangianeulerian vertical coordinate method in the mpas-ocean model[END_REF]. A difficulty in this case is the rezoning (a.k.a. regridding) phase to maintain the integrity of the grid locally and globally.

A tendency in the design of the oceanic dynamical cores is the extension to the non-hydrostatic (NH) equations. The most widely used approach nowadays in oceanic models is based on the incompressible NH system solved using a pressure correction/projection method which requires the solution of a 3D Poisson equation [START_REF] Lai | A nonhydrostatic version of fvcom: 1. validation experiments[END_REF][START_REF] Vitousek | Stability and consistency of nonhydrostatic free-surface models using the semi-implicit θ -method[END_REF][START_REF] Voltzinger | Nonhydrostatic tidal dynamics in the area of a seamount[END_REF]. Recently, [START_REF] Auclair | A non-hydrostatic non-boussinesq algorithm for free-surface ocean modelling[END_REF] have proposed the use of the compressible non-hydrostatic ocean equations with the advantage that no global algebraic system needs to be solved to compute NH pressure anomalies but with the disadvantage to permit acoustic modes. In this case a specific numerical procedure is required to maintain acceptable stability of the whole code.

Test strategy and benchmark suite

Given the wide variety of choices that need to be made during the development of dynamical cores and their overall complexity, it is crucial to define evaluation methods to compare the behavior of different models. Such effort has been made over the last decade by the global atmospheric community (Dynamical Core Model Intercomparison Project (DCMIP), [START_REF] Ullrich | Dynamical core model intercomparison project (dcmip) test case document[END_REF]. In particular, within DCMIP, a collection of test-cases that found broad acceptance in the community have been designed and applied by a large number of modeling groups. The workshop highlighted that in the context of the oceanic community, existing test-cases are scattered in the literature and not always fully documented and reproducible. The few existing examples of such effort (e.g. [START_REF] Ezer | Developments in terrain-following ocean models: intercomparisons of numerical aspects[END_REF][START_REF] Gerritsen | Validation document Delft3D-FLOW[END_REF][START_REF] Ilicak | Spurious dianeutral mixing and the role of momentum closure[END_REF][START_REF] Soufflet | On effective resolution in ocean models[END_REF] turned out to provide a valuable feedback on the consequences of model formulations. A good test-case should be easy to configure with analytical data suitable for all horizontal grids and different vertical coordinates and easy to evaluate while being relevant to test a given component of the dynamical core. The evaluation can be done either via analytical solutions (e.g. [START_REF] Bristeau | Various analytical solutions for the incompressible Euler and Navier-Stokes systems with free surface[END_REF], numerically converged solutions (provided that all models converge toward the same solution) or more subjectively based on an unambiguous physical understanding of the processes (e.g. [START_REF] Marques | Flow splitting in numerical simulations of oceanic dense-water outflows[END_REF]. Such benchmark suite is also useful to motivate communication between modeling groups and also to open room for prospective approaches from applied mathematicians to highlight their effectiveness. Existing initiatives aiming at oceanic model intercomparison based on realistic simulations (e.g. [START_REF] Chassignet | DAM ÉE-NAB: the base experiments[END_REF][START_REF] Griffies | Coordinated ocean-ice reference experiments (cores)[END_REF]) are generally too complex to clearly associate the observed differences to particular numerical choices. A way to evaluate numerical models in such complex configurations which could gain ground in the next few years is the uncertainty quantification (e.g. [START_REF] Iskandarani | An overview of uncertainty quantification techniques with application to oceanic and oil-spill simulations[END_REF]. Such approach provides tools to characterize modeling and numerical sensitivities.

Challenges

Throughout the 3 days of discussion, different current and future challenges have been identified. Addressing these challenges requires closer collaboration between the modeling groups.

Multi-resolution strategy: block-structured vs unstructured. On one hand, unstructured grid models have the ability to allow variable-resolution meshes provided an efficient mesh generation tool (e.g. Engwirda 2017), for example to adapt the resolution to follow the local Rossby radius [START_REF] Hallberg | Using a resolution function to regulate parameterizations of oceanic mesoscale eddy effects[END_REF][START_REF] Sein | Ocean modeling on a mesh with resolution following the local rossby radius[END_REF]). On the other hand structured grid models can also locally increase the resolution via nesting techniques [START_REF] Debreu | Two-way embedding algorithms: a review[END_REF][START_REF] Warner | Using a composite grid approach in a complex coastal domain to estimate estuarine residence time[END_REF][START_REF] Debreu | Two-way nesting in split-explicit ocean models: algorithms, implementation and validation[END_REF] or quadtree-octree refinement [START_REF] Popinet | A tree-based solver for adaptive ocean modelling[END_REF]. One advantage of the nesting approach is to allow the adjustment of the time-step and the physical parameters to the local resolution while unstructured models will need scale-aware parameterizations and a specific procedure for time refinement. Other points to investigate are the impact of variable resolution on propagating waves and the optimal layout to build a multiresolution mesh. A more prospective approach could be the use of an adaptive wavelet method [START_REF] Kevlahan | Adaptive wavelet simulation of global ocean dynamics using a new brinkman volume penalization[END_REF].

Energy consistency and resolved/unresolved scales coupling. Energy consistency is an important aspect for the proper interaction between resolved and parameterized scales (e.g. [START_REF] Burchard | Energy-conserving discretisation of turbulent shear and buoyancy production[END_REF][START_REF] Bachman | A scale-aware subgrid model for quasi-geostrophic turbulence[END_REF]. However as soon as a numerical core does not globally conserve energy at a discrete level (e.g. because of monotonicity enforcement, vertical remapping, or some form of upwinding), the identification of energy pathways is difficult and requires an in-depth analysis to close the energy cycle (e.g. [START_REF] Marsaleix | Energy conservation issues in sigmacoordinate free-surface ocean models[END_REF][START_REF] Eden | Closing the energy cycle in an ocean model[END_REF] which can be rather tedious (if not impossible) when advanced high-order numerics is used. An alternative could be to opt for an energy-conserving space and time discretization (e.g. [START_REF] Korn | Formulation of an unstructured grid model for global ocean dynamics[END_REF]Eldred et al. 2018) with specific care (based on explicit numerical dissipation in the dynamical core and/or on parameterization of subgrid processes) to avoid instability issues of the existing approaches [START_REF] Bell | Numerical instabilities of vector-invariant momentum equations on rectangular cgrids[END_REF]. This consistent coupling between dynamical cores and subgrid processes (known as physics-dynamics coupling) is an increasingly important topic for the building of geophysical models in general [START_REF] Gross | Physics-dynamics coupling in weather, climate and earth system models: Challenges and recent progress[END_REF].

Vertical coordinates & spurious numerical mixing. Spurious mixing (especially spurious dianeutral mixing) is a long-standing issue for oceanic dynamical cores (e.g. [START_REF] Griffies | Spurious diapycnal mixing associated with advection in a z-coordinate ocean model[END_REF]. The use of an ALE vertical coordinate is a way to mitigate this issue. Despite the fact that [START_REF] Gibson | Attribution of horizontal and vertical contributions to spurious mixing in an arbitrary lagrangianeulerian ocean model[END_REF] suggest that the vertical component of spurious mixing dominates as horizontal resolution increases it should not overshadow that many components of dynamical cores can be a source of numerical mixing (e.g. horizontal advection, time-stepping, the stabilization of the mode-splitting procedure). There is still a need to better understand the implications of different choice of momentum / tracers advection schemes, rezoning and remapping procedures on numerical mixing. Idealized test-cases and efficient diagnostic tools are important to tackle this issue (e.g. [START_REF] Ilicak | Spurious dianeutral mixing and the role of momentum closure[END_REF][START_REF] Klingbeil | Quantification of spurious dissipation and mixing discrete variance decay in a finite-volume framework[END_REF]. It would also be instructive to keep investigating the improvement of quasi-Eulerian vertical coordinates (e.g. [START_REF] Berntsen | A perfectly balanced method for estimating the internal pressure gradients in σ -coordinate ocean models[END_REF]. In this context schemes for internal pressure gradient and for isoneutral tracer diffusion should also be considered (e.g. [START_REF] Shao | Improvements to an Extrema-diminishing, Density-preserving Lateral Diffusion Algorithm[END_REF].

Non-hydrostatic pressure contribution. As discussed earlier, many oceanic dynamical cores have the possibility to account for non-hydrostatic effects. A difficulty is now to clearly identify under which conditions relaxing the hydrostatic assumption is necessary and which resolution is required for proper NH modeling. Another challenge is the possibility to locally account for NH effects within a primitive equations model either in the form of a superparameterization (e.g. [START_REF] Campin | Superparameterization in ocean modeling: Application to deep convection[END_REF] or in the form of 2-way nesting between coarse hydrostatic and fine nonhydrostatic meshes (e.g. [START_REF] Blayo | About interface conditions for coupling hydrostatic and nonhydrostatic navier-stokes flows[END_REF]. Ultimately, further investigations of the merits and flaws of the incompressible vs compressible NH approaches would be worthwhile, also for global applications [START_REF] Losch | How sensitive are coarse general circulation models to fundamental approximations in the equations of motion?[END_REF].

Coupling with other Earth-System compartments. Oceanic dynamical cores are often used as a component of larger coupled model systems. Coupling to several other Earth system compartments is common, such as to surface wave, sea-ice, atmospheric, biogeochemical, benthic, and hydrological models. The numerical implementation of such coupling can become an issue (e.g. Lemarié et al. 2015;[START_REF] Beljaars | On the numerical stability of surface-atmosphere coupling in weather and climate models[END_REF] especially at high coupling frequency and/or spatial resolution. More systematic analysis of the coupling stability and consistency using simplified equation sets and the design of simplified coupled test-cases must be encouraged.

Vanishing layers, wetting and drying and shockresolving numerics. An accurate treatment of wetting and drying is essential for coastal simulations as well as for climate simulations of under-ice-shelf cavities. At a numerical level this requires the non-negativity of the water height and an adequate volume conserving treatment of dry states (a.k.a. vacuum states). Standard numerical methods used in oceanic dynamical cores do not have the ability to handle vacuum, or equivalently shocks. Instead, approaches based on some pre-defined minimum water depth and specific ad-hoc manipulation of discrete fluxes are often used (e.g. Sec. 5.2 in [START_REF] Klingbeil | The numerics of hydrostatic structured-grid coastal ocean models: State of the art and future perspectives[END_REF]). However there is a vast literature dedicated to the design of numerical schemes preserving positivity and able to correctly treat vacuum states which furthermore satisfy an entropy-preserving property; i.e. the nonlinear solution is physically relevant even in the presence of discontinuities (e.g. [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for Shallow Water flows[END_REF][START_REF] Audusse | Kinetic entropy inequality and hydrostatic reconstruction scheme for the Saint-Venant system[END_REF]. Considering advectiondiffusion equations e.g. for tracers (temperature, salinity) there is an obvious benefits in using numerical schemes preserving the maximum principle. There could be an interest in comparing these more advanced approaches with the usual treatment adopted in dynamical cores.

Conclusions

The workshop gave a broad and fresh overview of existing numerical methods used in realistic ocean models as well as some examples of alternatives from the applied maths community. The participants have been enthusiastic and very positive about the possibility to sustain this type of workshop into a biennial workshop series. A collective article is currently in preparation to summarize the challenges and prospects for oceanic numerical cores across all scales. Moreover, a particular effort will be directed toward the formation of an active community of model developers with international collaborations, starting e.g. with the standardization of existing idealized test cases as the basis for model/methods intercomparison studies. The next meeting will be organized either in Fall 2019 or Winter 2020 in Hamburg, Germany.

Acknowledgments. The authors thank Gurvan Madec for reading and checking the content of this summary as well as all the workshop's participants 1 . The workshop was hosted by the French National Institute for computer science and applied mathematics (Inria), Paris, France. We gratefully acknowledge the support from Inria, from 1 https://commodore2018.sciencesconf.org/resource/ listeparticipants the Coastal and Regional Ocean Community Model research group (Gdr Croco), and from the french program Earths Fluid Envelopes (LEFE).