Mélanie Otténio

Michel Destrade

Raymond W Ogden

Acoustic waves at the interface of a pre-stressed incompressible elastic solid and a viscous fluid

Acoustic waves at the interface of a pre-stressed incompressible elastic solid and a viscous fluid

Introduction

Seismic records show that underground rocks and ocean beds are subject to stress and strain and that surrounding fluids are viscous and under high pressures. Clinical ultrasonic measurements indicate that arteries can undergo large strains in service and are most sensitive to changes in blood pressure. Many moving and vibrating parts of automotive devices are made of loaded elastomers in contact with highly viscous fluids. These are a few examples of situations where it is crucial to model and understand the motions and the stability of the interface between a deformed elastic solid and a viscous fluid. Yet only a handful of studies can be found on the subject, especially when compared with the abundant literature on waves at the interface between an elastic solid and an inviscid fluid, which goes from the pioneering works of Galbrun [START_REF] Galbrun | Propagation d'une Onde Sonore dans l'Atmosphère et Théorie des Zones de Silence[END_REF], Cagniard [START_REF] Cagniard | Réflexion et Eéfraction des Ondes Sismiques Progressives[END_REF], Scholte [START_REF] Scholte | On the Stoneley wave equation[END_REF], and Biot [START_REF] Biot | The interaction of Rayleigh and Stoneley waves in the ocean bottom[END_REF] to the definitive treatment of the acoustoelastic effect by Sinha et al. [START_REF] Sinha | Stoneley and flexural modes in a pressurized borehole[END_REF]; see also Poirée [6] and Degtyar and Rohklin [START_REF] Degtyar | Stress effect on boundary conditions and elastic wave propagation through an interface between anisotropic media[END_REF]. Waves at the interface between a viscous fluid and an undeformed isotropic elastic solid were examined by Vol'kenshtein and Levin [START_REF] Vol'kenshtein | Structure of a Stoneley wave at an interface between a viscosus fluid and a solid[END_REF] and the corresponding problem for an anisotropic elastic solid by Wu and Wu [START_REF] Wu | Surface waves in coated anisotropic medium loaded with viscous liquid[END_REF]. To the best of our knowledge, only Bagno, Guz, and their co-workers have studied the title problem (see, for example, [START_REF] Bagno | Stonely waves on the contact boundary between a prestressed incompressible rigid half-space and a viscous compressible fluid[END_REF], [START_REF] Bagno | Elastic waves in prestressed bodies interacting with a fluid[END_REF]). Their analytical treatment is, however, quite succinct and we therefore aim to shed new light on the problem by re-examining it on the basis of recent developments in the theory of small-amplitude waves, linearized in the neighbourhood of a finite, static, homogeneous deformation.

It turns out that one of the trickiest aspects of the study is the derivation of proper incremental boundary conditions at the interface because these are usually written in terms of the nominal stress in a deformed solid (Lagrangian formulation, Section 2) and in terms of the Cauchy stress in a fluid (Eulerian formulation, Section 3). These equations are combined in an appropriate way for a general interface in Section 4. We then specialize the analysis to principal wave propagation for an arbitrary (incompressible, isotropic) strainenergy function in Section 5. In the course of the analysis in Section 5, by way of application of the theory, we show that in respect of an abdominal aortic aneurysm it is appropriate to neglect the curvature and finite thickness for ultrasonic waves (10 MHz), i.e. to treat the aneurysm locally as a half-space. It also shown that it is not appropriate to use waves with a real frequency to study the stability of compressed solids in contact with a viscous fluid. Finally, in Section 6, by specializing to the neo-Hookean solid, we consider the propagation of non-principal waves for both tension and compression of the half-space in order to illustrate the influence of the fluid.

Basic equations for the solid

For the solid material we denote by F the deformation gradient relating the stress-free reference configuration, denoted B 0 , to the finitely deforming configuration, denoted B. This has the form F = Gradx, where x = χ(X, t) is the position vector in B at time t of a material point located at X in B 0 , χ is the deformation mapping, and Grad is the gradient operator relative to B 0 .

We consider the material to be elastic with a strain-energy function, defined per unit volume, denoted by W = W (F ). Furthermore, we restrict attention to incompressible materials so that the constraint

det F = 1 (2.1)
is in force. The nominal stress tensor, here denoted by S, and the Cauchy stress tensor σ are then given by

S = ∂W ∂F -pF -1 , σ = F ∂W ∂F -pI, (2.2)
where p is a Lagrange multiplier associated with the constraint (2.1) and I is the identity tensor.

The equation of motion is

DivS = ρẍ, (2.3)
where Div is the divergence operator relative to B 0 , ρ is the mass density of the material, and a superposed dot signifies the material time derivative. Next, we consider a small motion superimposed on the finite deformation. Let u(x, t) be the displacement vector relative to B and v(x, t) = u the associated particle velocity (the material time derivative of u). Then, on taking the increment of equation (2.3) and thereafter changing the reference configuration from B 0 to B, we obtain

divs = ρü ≡ ρ v, (2.4)
where s is the increment in S (referred to B) and div the divergence operator relative to B. The (linearized) incremental version of the constitutive relation (2.2) is written s = A(gradu) + p(gradu) -pI, (

where A is a fourth-order tensor of elastic moduli, grad is the gradient operator relative to B, and p is the increment in p. In component form, this is written

s ij = A ijkl u l,k + pu i,j -pδ ij , (2.6)
where ,k denotes ∂/∂x k and δ ij is the Kronecker delta. In terms of W the components of A are given by

A ijkl = F ip F kq ∂ 2 W ∂F jp ∂F lq . (2.7)
For details of these derivations (in a slightly different notation) we refer to Dowaikh and Ogden [START_REF] Dowaikh | On surface waves and deformations in a pre-stressed incompressible elastic solid[END_REF]. We now consider the material to be isotropic, so that W = W (λ 1 , λ 2 , λ 3 ) is a symmetric function of the principal stretches, λ 1 , λ 2 , λ 3 (the positive square roots of the principal values of F F T , where T signifies the transpose), subject to the constraint

λ 1 λ 2 λ 3 = 1, (2.8)
which follows from (2.1). Then, on noting that for an isotropic material σ is coaxial with F F T and specializing equation (2.2) 2 , we obtain the principal Cauchy stresses (see, for example, Ogden [START_REF] Ogden | Non-Linear Elastic Deformations[END_REF])

σ i = -p + λ i W i , i = 1, 2, 3 (no sum over i), (2.9)
where

W i = ∂W/∂λ i , i = 1, 2, 3.
When referred to the same principal axes, the only non-zero components of A are

A iijj = λ i λ j W ij , i, j ∈ {1, 2, 3}, (2.10)
A ijij = λ i W i -λ j W j λ 2 i -λ 2 j λ 2 i , i, j ∈ {1, 2, 3}, i = j, (2.11)
A ijji = A jiij = A ijij -λ i W i , i, j ∈ {1, 2, 3}, i = j, (2.12)
where W ij = ∂ 2 W/∂λ i ∂λ j . For subsequent convenience, we adopt the notations defined by

γ ij = (λ i W i -λ j W j )λ 2 i /(λ 2 i -λ 2 j ), β ij = (λ 2 i W ii + λ 2 j W jj )/2 -λ i λ j W ij + (λ i W j -λ j W i )λ i λ j /(λ 2 i -λ 2 j ), (2.13) noting that γ ji λ 2 i = γ ij λ 2 j and β ji = β ij .

The pre-stressed elastic half-space

We now consider B to be independent of time and to correspond to a pure homogeneous strain of a half-space defined by x 2 0. The half-space is maintained in this configuration so that its boundary x 2 = 0 is a principal plane of strain. We denote by x 1 and x 3 the other two principal directions of strain and by λ 1 , λ 2 , λ 3 the principal stretches in the x 1 , x 2 , x 3 directions, respectively. The corresponding principal Cauchy stresses are then as given by (2.9). In particular, the boundary x 2 = 0 is subject to a normal stress σ 2 and, after elimination of p, the other two principal Cauchy stresses are then given by

σ 1 = σ 2 + λ 1 W 1 -λ 2 W 2 , σ 3 = σ 2 + λ 3 W 3 -λ 2 W 2 .
(

We are interested in the propagation of incremental (small amplitude) acoustic waves along the boundary plane x 2 = 0, in a direction making an angle θ with the principal direction x 1 . The incremental velocity and nominal stress fields v and s are then considered as superimposed on this finite static configuration. We examine inhomogeneous time-harmonic plane waves of the form

{v, s}(x 1 , x 2 , x 3 , t) = {v(x 2 ), -(k/ω)ŝ(x 2 )}e ik(c θ x 1 +s θ x 3 ) e -iωt , (2.15)
where we have introduced the notations c θ = cos θ, s θ = sin θ, k is the wave number, ω is the wave frequency, and v, ŝ are functions of x 2 only, such that

v(∞) = 0, ŝ(∞) = 0. (2.16)
Using the results of Destrade et al. [START_REF] Destrade | Nonprincipal surface waves in deformed incompressible materials[END_REF] (see also Chadwick [START_REF] Chadwick | The application of the Stroh formalism to prestressed elastic media[END_REF]), we find that the incremental equations of motion can be written as a first-order differential system of six equations, namely

ξ ′ (x 2 ) = ikN ξ(x 2 ), (2.17)
where the notation ξ is defined by

ξ = [v 1 , v2 , v3 , ŝ21 , ŝ22 , ŝ23 ] T , (2.18)
and the 6 × 6 matrix N has the block structure

N = N 1 N 2 N 3 + ρI N T 1 , (2.19)
in which the 3×3 matrices N 1 , N 2 , N 3 are real and their components depend on the material parameters γ ij and β ij given in (2.13), and the notation ρ = ρω 2 /k 2 has been introduced. Here I represents the 3 × 3 identity matrix.

Explicitly, -N 1 , N 2 , -N 3 are   0 c θ (γ 21 -σ 2 )/γ 21 0 c θ 0 s θ 0 s θ (γ 23 -σ 2 )/γ 23 0   ,   1/γ 21 0 0 0 0 0 0 0 1/γ 23   ,   η 0 -κ 0 ν 0 -κ 0 µ   , (2.20)
respectively, where

η = 2c 2 θ (β 12 + γ 21 -σ 2 ) + s 2 θ γ 31 , ν = c 2 θ [γ 12 -(γ 21 -σ 2 ) 2 /γ 21 ] + s 2 θ [γ 32 -(γ 23 -σ 2 ) 2 /γ 23 ], µ = c 2 θ γ 13 + 2s 2 θ (β 23 + γ 23 -σ 2 ), κ = c θ s θ (β 13 -β 12 -β 23 -γ 21 -γ 23 + 2σ 2 ).
(2.21) Equation (2.17) provides the general expression for the equations of motion, for arbitrary θ and W . Now, in seeking a decaying partial-mode solution of the form

ξ(x 2 ) = e -ksx 2 ζ, ℜ(ks) > 0, (2.22)
where ζ is a constant vector and s an unknown scalar, we arrive at the eigenvalue problem N ζ = isζ. In general, the associated propagation condition, det(N -isI) = 0, is a cubic in s 2 [START_REF] Rogerson | Harmonic wave propagation along a non-principal direction in a pre-stressed elastic plate[END_REF], where now I is the 6 × 6 identity matrix. Its analytical resolution is too cumbersome to be of practical interest, and so we specialize the general equations to the following, simpler, situations: (i) principal wave propagation (θ = 0) for arbitrary W ; (ii) nonprincipal wave propagation (θ = 0) for the neo-Hookean material, for which

W = C(λ 2 1 + λ 2 2 + λ 2 3 -3)/2, ( 2.23)
where C > 0 is a constant (the shear modulus of the material in the reference configuration).

In Case (i), the equations of motion decouple the system v′

3 ŝ′ 23 = ik 1/γ 21 0 0 1/γ 23 v3 ŝ23 , (2.24)
(for which the trivial solution may be chosen) from a system of four differential equations for v1 , v2 , ŝ21 , ŝ22 . Hence, in this case, the wave is elliptically polarized, in the (x 1 , x 2 ) plane. The corresponding propagation condition is a quadratic in s 2 , which can be solved explicitly.

In Case (ii), we also find that the wave is two-partial, polarized in the plane containing the directions of propagation and attenuation (the saggital plane); there, the corresponding propagation condition involves the product of the factor (s 2 -1) and a term linear in s 2 , which simplifies the analysis.

Before embarking on the details of these cases, we complete the description of the boundary-value problem by considering the behaviour of the wave in the fluid in the half-space x 2 ≤ 0.

The fluid half-space

Adjoining the deformed solid half-space is a half-space x 2 0 filled with an incompressible viscous Newtonian fluid, for which all mechanical fields are denoted by a superscripted asterisk. In the static state the fluid is subject only to a hydrostatic stress σ * = -P * I, and by continuity of traction across the boundary x 2 = 0 we must have

-P * = σ 2 . (3.1)
The constitutive law for the fluid associated with the motion is then written in terms of a superimposed Cauchy stress tensor, denoted here by s * and given by

s * = -p * I + 2µ * D * , trD * = 0, (3.2)
where µ * is the viscosity of the fluid,

D * = 1 2 [gradv * + (gradv * ) T ], (3.3)
v * is the fluid velocity, and p * = p * (x, t).

We seek inhomogeneous waves with the same structure as in the solid, that is

{v * , s * }(x 1 , x 2 , x 3 , t) = {v * (x 2 ), -(k/ω)ŝ * (x 2 )}e ik(c θ x 1 +s θ x 3 ) e -iωt , (3.4)
where v * , ŝ * are functions of

x 2 only, such that v * (-∞) = 0, ŝ * (-∞) = 0. (3.5)
We find that the equations of motion, divs * = ρ * v * (where ρ * is the mass density of the fluid), linearized in v * , can be cast as

ξ * ′ (x 2 ) = ikN * ξ * (x 2 ), (3.6)
where

ξ * = [v * 1 , v * 2 , v * 3 , ŝ * 21 , ŝ * 22 , ŝ * 23 ] T , (3.7)
and the constant complex matrix N * has the block structure

N * = N * 1 N * 2 N * 3 + ρ * I N * 1 , (3.8)
N * 1 , N * 2 , N * 3 being real symmetric matrices, and the notation ρ * = ρ * ω 2 /k 2 has been adopted. If we write μ * = µ * ω, then, respectively,

-N * 1 , -iμ * N * 2 , -iN * 3 /μ * are   0 c θ 0 c θ 0 s θ 0 s θ 0   ,   1 0 0 0 0 0 0 0 1   ,   4c 2 θ + s 2 θ 0 3c θ s θ 0 0 0 3c θ s θ 0 c 2 θ + 4s 2 θ   . (3.9)
Again, when we seek a decaying partial-mode solution, this time in the form

ξ * (x 2 ) = e ks * x 2 ζ * , ℜ(ks * ) > 0, (3.10)
where ζ * is a constant vector and s * an unknown scalar, we end up with an eigenvalue problem, here

N * ζ * = -is * ζ * . The associated propagation condition is det(N * + is * I) = 0, which here simplifies to (s ⋆2 -1)(s ⋆2 -1 + iǫ) 2 = 0, (3.11) with roots ±1, ± √ 1 -iǫ (repeated), (3.12)
where ǫ = ρ * /μ * = ρ * ω/(µ * k 2 ). The roots are independent of θ, as expected, because the fluid is isotropic. Corresponding to each of the four roots, there are four eigenvectors and therefore potentially four partial-modes. However, two of these must be discarded since their amplitudes do not decay with distance from the interface x 2 = 0. The two remaining modes form the basis for the general solution of the equations of motion that is needed for matching with the two-partial wave in the solid. We now give the general boundary conditions at the deformed solid/viscous fluid interface.

The interface

In order to match the incremental tractions across the boundary it is necessary to work in terms of the Cauchy stress since the nominal stress is not defined inside the fluid. Towards this end we first calculate the incremental traction in the solid in terms of the Cauchy stress. Continuity of traction requires

S T N dA = σnda = σ * nda, (4.1)
where dA and da are the area elements in B 0 and B, respectively. Taking the increment of this and updating the reference configuration to B yields

s T nda ≡ σnda + σ nda = s * nda + σ * nda, (4.2)
where a superposed tilde indicates an increment. Note that, after updating, F = I and S = σ in the configuration B. Now, according to Nanson's formula (applied to the boundary of the solid), we have nda = F -T N dA, from which it follows, again after updating, that nda = -(gradu) T nda.

Hence, the incremental traction continuity condition can be written

s T n ≡ [σ -σ(gradu) T ]n = [σ * -σ * (gradu) T ]n, (4.4)
and we recall that σ * = -P * I.

Since n is in the x 2 direction for the considered half-space we may write the continuity condition in component form as

s 2i = s * 2i + P * u 2,i , i = 1, 2, 3, on x 2 = 0. (4.5)
Additionally, the velocity must be continuous, i.e.

v * i = v i , i = 1, 2, 3, on x 2 = 0. (4.6)
In terms of the functions v(x 2 ), ŝ(x 2 ) and their counterparts in the fluid, the boundary conditions become

v * i (0) = vi (0), i = 1, 2, 3, (4.7)
and, noting that

v i = -iωu i , ŝ * 12 (0) + c θ P * v2 (0) = ŝ21 (0), ŝ * 32 (0) + s θ P * v2 (0) = ŝ23 (0), (4.8)
and ŝ * 22 (0) -

i k P * v′ 2 (0) = ŝ22 (0). (4.9)
5 Principal waves: no restriction on W

General solution in the solid

Here we take θ = 0 and place no restriction on the form of W . When θ = 0, the in-plane mechanical fields in the solid satisfy the equations of motion ξ ′ = ikN ξ, where now ξ(x 2 ) = [v 1 , v2 , ŝ21 , ŝ22 ] T and

N =     0 -1 + σ 2 1/γ 21 0 -1 0 0 0 ρ -η 0 0 -1 0 ρ -ν -1 + σ 2 0     , (5.1)
with σ 2 = σ 2 /γ 21 and η and ν now reduced to

η = 2[β 12 + γ 21 (1 -σ 2 )], ν = γ 12 + γ 21 (1 -σ 2 ) 2 .
Here the propagation condition is a quadratic in s 2 [START_REF] Dowaikh | On surface waves and deformations in a pre-stressed incompressible elastic solid[END_REF] and given by

s 4 -2βs 2 + α 2 = 0, (5.2)
where 2β = 2β 12 - ρ γ 21 , α 2 = γ 12 - ρ γ 21 . (5.3)
We recall that ρ = ρω 2 /k 2 . Formally, the roots of the quartic (5.2) are

± β + β 2 -α 2 , ± β -β 2 -α 2 .
(5.4)

We pause the analysis to highlight a particular feature of the present interface waves. Because the fluid is viscous, the wave number k is complex and so, therefore, are β and α. It follows that, in contrast to the purely elastic case [START_REF] Dowaikh | Interfacial waves and deformations in pre-stressed elastic media[END_REF], it is not clear a priori which two of these four roots are such that the decay condition Eq. (2.22) 2 is satisfied. Let s 1 and s 2 be two such roots. We note first that

s 2 1 + s 2 2 = 2β, s 2 1 s 2 2 = α 2 , (5.5)
and then, depending on which value in Eq. ( 5.4) they correspond to, one of the following four possibilities may arise:

s 1 s 2 = α, s 1 + s 2 = 2(β + α), s 1 s 2 = α, s 1 + s 2 = -2(β + α), s 1 s 2 = -α, s 1 + s 2 = 2(β -α), s 1 s 2 = -α, s 1 + s 2 = -2(β -α). (5.6)
Now, the eigenvectors ζ 1 and ζ 2 associated with s 1 and s 2 are columns of the matrices adjoint to N -is 1 I and N -is 2 I, respectively. Taking, for instance, the fourth column gives

ζ i = [is i , -1, -γ 21 (1 + s 2 i -σ 2 ), -iγ 21 s i (1 + 2β -s 2 i -σ 2 )] T , i = 1, 2.
(5.7) Finally, the general solution in the solid is

ξ(x 2 ) = A 1 e -ks 1 x 2 ζ 1 + A 2 e -ks 2 x 2 ζ 2 , (5.8)
where A 1 and A 2 are constants.

General solution in the fluid

In the fluid, a two-partial solution is required for adequate matching at the interface. Taking θ = 0 in Section 3, we find that the in-plane equations of motion are ξ * ′ = ikN * ξ * , where now

ξ * (x 2 ) = [v * 1 , v * 2 , ŝ * 21 , ŝ * 22 ] T and N * =     0 -1 i/μ * 0 -1 0 0 0 4iμ * + ρ * 0 0 -1 0 ρ * -1 0     .
(5.9)

We assume that the wave propagates and is attenuated in the direction x 1 0. From Eq. (3.4) with θ = 0, we see that these assumptions lead to

ℜ(k) > 0, ℑ(k) > 0.
(5.10)

On the other hand, the wave is also attenuated with distance from the interface; it follows that we can discard the root -1 from Eq. (3.12) and retain the root +1. The choice to be made for the remaining two roots in Eq. (3.12) is not so clear cut and for the time being we call s * the suitable root; hence s * is such that s * 2 = 1 -iǫ, ℜ(ks * ) > 0.

(5.11)

Finally, the general solution in the fluid is

ξ * (x 2 ) = A * 1 e kx 2     i 1 -2iμ * -μ * (1 + s * 2 )     + A * 2 e ks * x 2     is * 1 -iμ * (1 + s * 2 ) -2μ * s *     , (5.12)
where A * 1 and A * 2 are constants.

Dispersion equation for the interface wave

When we specialize the general boundary conditions (4.7), (4.8) and (4.9) to the present context, we find a linear homogeneous system of four equations for the four unknowns A 1 , A 2 , A * 1 , A * 2 , for which the associated determinant must be zero. After some manipulations, using Eq. (5.5), we find that

-s 1 -s 2 1 s * 1 1 1 1 γ 21 (1 + s 2 1 ) γ 21 (1 + s 2 2 ) -2iμ * -iμ * (1 + s ⋆2 ) γ 21 s 1 (1 + s 2 2 ) γ 21 s 2 (1 + s 2 1 ) iμ * (1 + s ⋆2 ) 2iμ * s * = 0. ( 5.13)
We see at once that the normal load σ 2 does not appear explicitly in this equation. This feature highlights a major difference between a wave at the interface of a loaded solid half-space and a vacuum (Rayleigh wave) and a wave at the interface of a loaded solid half-space and a viscous fluid, as considered here; Chadwick and Jarvis [START_REF] Chadwick | Interfacial waves in a pre-strained neo-Hookean body I. Biaxial states of strain[END_REF] also noted this peculiarity for waves at the interface of two loaded solid half-spaces (Stoneley wave). Of course, σ 2 still plays an important role, in particular in the determination of the pre-stretch ratios and of the amplitudes of the tractions in the solid.

Once the determinant is expanded and the factors (s 1 -s 2 )(1 -s * ) are removed, we end up with the exact dispersion equation,

γ 2 21 γ 12 - ρ γ 21 + 2β 12 + 2γ 21 - ρ γ 21 s 1 s 2 -1 -iγ 21 μ * [2(1 -s 1 s 2 )(1 -s * ) + (s 1 + s 2 )(1 + s * )(s * + s 1 s 2 )] + μ * 2 [iǫ + (iǫ -4) s * ] = 0. (5.14)
This equation is fully explicit because the terms s 1 + s 2 , s 1 s 2 , and s * are given by Eqs. (5.6) and Eq. (5.11). Of course, as noted earlier, there are several possibilities for these terms, which generate in total eight different dispersion equations. In each case the resulting root(s) in k must be checked for validity against the propagation and decay conditions (2.22) 2 , (5.10), and (5.11) 2 , which we summarize here as

ℜ(k) > 0, ℑ(k) > 0, ℜ(ks * ) > 0, ℜ(ks 1 ) > 0, ℜ(ks 2 ) > 0. (5.15)
Notice that in the special case of solids whose strain-energy function W is such that 2β 12 = γ 12 + γ 21 , which includes the neo-Hookean solid (2.23), the biquadratic (5.2) factorizes as (s 2 -1)(s 2 -α 2 ) = 0. Hence s 1 = 1, s 2 = ±α, and the exact secular equation (5.14) simplifies accordingly, leading this time to four different exact dispersion equations.

Application: modelling of intravascular ultrasound

In recent years, intravascular ultrasound (IVUS) has proved to be a most promising tool of investigation for measuring and assessing abdominal aortic aneurysms. Its accuracy is as good as that of computed tomography (CT) scans and it has obvious non-radiative advantages [START_REF] Lederle | Variability in measurement of abdominal aortic aneurysms[END_REF], [START_REF] Wolf | Duplex ultrasound scanning versus computed tomographic angiography for postoperative evaluation of endovascular abdominal aortic aneurysm repair[END_REF], [START_REF] Garret | Intravascular ultrasound aids in the performance of endovascular repair of abdominal aortic aneurysm[END_REF]. We now apply the results of this section to an IVUS context. First we recall that medical ultrasound imaging devices operate in the 1-10 MHz range; accordingly we take ω = 10 7 Hz. We argue, and will check a posteriori, that at such high frequency the wavelength is small compared with the radius and thickness of an artery so that, as far as the propagation of localized waves is concerned, an aortic aneurysm can be modelled as a half-space. Here we take the axis of the artery along the x 1 direction and consider that the half-space x 2 0 is filled with blood.

For the solid, we use a strain-energy function devised by Raghavan and Vorp [START_REF] Raghavan | Toward a biomechanical tool to evaluate rupture potential of abdominal aortic aneurysm: identification of a finite strain constitutive model and evaluation of its applicability[END_REF] to fit experimental data collected on uniaxial tension tests of aortic aneurysms, namely

W = C 1 (λ 2 1 + λ 2 2 + λ 2 3 -3) + C 2 (λ 2 1 + λ 2 2 + λ 2 3 -3) 2 , (5.16)
where, typically, C 1 = 0.175 MPa, C 2 = 1.9 MPa. We assume that the aneurysm corresponds to a region where the tissue undergoes a strain of 20%, and that the end-systolic blood pressure is high, 150 mmHg (= 20 kPa) say, so we take σ 2 = 20 kPa. We consider that the tissue is free to expand or contract in the x 3 direction with σ 3 = 0. From (2.14) we find the remaining stresses and strains. Summarizing, we have λ 1 = 1.2, λ 3 ≃ 0.908, σ 1 ≃ 71.45 kPa, σ 2 = 20 kPa, σ 3 = 0, (5.17)

with λ 2 calculated from the incompressibility condition. For the mass density, we take ρ = 1000 kg/m 3 . For the blood, we use typical values of viscosity [START_REF] Spring | Decreased wall shear stress in the common carotid artery of patients with peripheral arterial disease or abdominal aortic aneurysm: Relation to blood rheology, vascular risk factors, and intima-media thickness[END_REF] and mass density [START_REF] Kenner | The continuous highprecision measurement of the density of flowing blood[END_REF]: µ * = 3.5 × 10 -3 Ns/m 2 , ρ * = 1050 kg/m 3 .

Here, the frequency ω is fixed as a real quantity a priori. We then replace k everywhere by k = ωS, where S = S + + iS -is the (complex) slowness, and the only unknown in the dispersion equation (5.14). The propagation and decay conditions (5.15) are satisfied and may be written as

ℜ(S) > 0, ℑ(S) > 0, ℜ(Ss * ) > 0, ℜ(Ss 1 ) > 0, ℜ(Ss 2 ) > 0. (5.18)
We find that the only qualifying root is S = 2.721 × 10 -2 + 4.424 × 10 -4 i, (5.19) from which we deduce the phase speed v = 1/S + , the damping factor γ = S -, and the wavelength λ 0 = 2π/(ωS + ): v = 36.75 m/s, γ = 4.424 × 10 -4 m -1 , λ 0 = 23.1 µm.

(5.20)

Then we plot the depth profiles of the wave. Its amplitude is the real part of v in the solid and of v * in the fluid. Explicitly,

ℜ (v(x 1 , x 2 , t)) = e -γx 1 [ℜ (v(x 2 )) cos ω(x 1 /v -t)] .
-ℑ (v(x 2 )) sin ω(x 1 /v -t)] , (5.21) in the solid, and similarly for the fluid, for which v and v are replaced by their asterisk counterparts. Clearly, the particle velocity is elliptically polarized, and the lengths of the ellipse semi-axes decay with distance away from the interface and also with increasing x 1 . Figure 1 shows the variations of the normal (continuous curve) and tangential (dotted curve) velocity components in the fluid (x 2 0) and in the solid (x 2 0), normalized with respect to v2 (0), as functions of x 2 /λ 0 . Note that the components are in phase quadrature. Note also, as is clear from the zooms shown in Figure 1, that the components are continuous across the interface, as expected, but their first derivatives are discontinuous. The wave is elliptically polarized near the interface; the major axis is normal to the interface and more than 12 times the length of the minor axis; the ellipse is described in a retrograde manner. This shape is carried through the depth of the solid whereas in the fluid, the wave becomes rapidly circularly polarized at a depth of about 0.06 wavelengths, and remains nearly so through the rest of the half-space. The localization is greater in the fluid than in the solid: the amplitude has almost vanished after one wavelength into the former and after five wavelengths in the latter. An aneurysm is typically 1 mm thick, which, with the numerical values used here, is more than 50 wavelengths; thus, the assumption of a semi-infinite solid is justified, as is the assumption of a flat interface (aneurysms are typically of diameter 5 cm).

Finally, we note that if the fluid is absent then the corresponding surface wave would travel with speed 40.734 m/s; hence the viscous fluid not only dampens the wave but also slows it down noticeably. Of course, blood flows in an artery, and creates shear deformation and stress in the solid. The influence of wall shear stress on the waves will be treated elsewhere.

Example: compressive stresses

Here we consider the behaviour of the elastic half-space (in contact with the viscous fluid) when it is compressed in the x 1 direction (λ 1 < 1). In their pioneering works, Green and Zerna [START_REF] Green | Theoretical Elasticity[END_REF] and Biot [START_REF] Biot | Surface instability of rubber in compression[END_REF], [START_REF] Biot | Mechanics of Incremental Deformations[END_REF] showed that when a highly elastic half-space with a free surface (no fluid contact) is compressed, a surface instability may develop. Focusing on neo-Hookean solids (2.23) they showed that the critical stretches are λ cr = 0.666 for tangential equibiaxial compression (λ 1 = λ, λ 2 = λ -2 , λ 3 = λ), λ cr = 0.544 for plane strain compression (λ 1 = λ, λ 2 = λ -1 , λ 3 = 1), and λ cr = 0.444 for normal equibiaxial compression (λ 1 = λ, λ 2 = λ -1/2 , λ 3 = λ -1/2 ). Theirs was a static stability analysis, later also included in a wider dynamical context by Flavin [START_REF] Flavin | Surface waves in pre-stressed Mooney material[END_REF], Willson [START_REF] Willson | Plate waves in Hadamard materials[END_REF], Chadwick and Jarvis [START_REF] Chadwick | Surface waves in a pre-stressed elastic body[END_REF], Dowaikh and Ogden [START_REF] Dowaikh | On surface waves and deformations in a pre-stressed incompressible elastic solid[END_REF], and others. We now consider this problem.

A localized small-amplitude wave propagates over the free surface of a deformed Mooney-Rivlin or neo-Hookean half-space with normalized squared speed ρv 2 /γ 12 = 1 -σ 2 0 λ 2 2 /λ 2 1 , where σ 0 is the real root of σ 3 + σ 2 + 3σ -1 = 0 (σ 0 = 0.2956). Clearly, in the examples of plane strain and equi-biaxial strain above, the squared wave speed increases when λ increases and decreases when λ decreases. At λ = 1 there is no pre-strain and ρv 2 /C = 0.9126, the value found by Lord Rayleigh [START_REF] Rayleigh | On waves propagated along the plane surface of an elastic solid[END_REF] for isotropic incompressible linearly elastic solids. As λ → ∞, the squared wave speed tends to the squared wave speed γ 12 /ρ of a transverse bulk wave. As λ decreases, there is a critical stretch λ cr at which the squared speed is zero and below which v 2 < 0. Since the wave time dependence is of the form e ik(x 1 -vt) , with k > 0, it follows that the amplitude grows without bound in time when λ < λ cr , and that the surface becomes unstable (at least, in the linearized theory). It is therefore appropriate to ask what happens when the compressed half-space is in contact with a viscous fluid.

Bagno and co-workers addressed this question in a series of articles (see the review by Bagno and Guz [START_REF] Bagno | Elastic waves in prestressed bodies interacting with a fluid[END_REF] and references therein). For a neo-Hookean solid they found that the wave speed drops to zero when λ = 0.544 in plane strain compression and when λ = 0.444 in normal equi-biaxial compression, that is at the same critical stretches as for surface (solid/vacuum) instability, irrespective of the viscous fluid characteristics. On the other hand, Bagno and Guz [START_REF] Bagno | Stonely waves on the contact boundary between a prestressed incompressible rigid half-space and a viscous compressible fluid[END_REF] find that the wave speed falls to zero at a critical stretch that does depend upon the viscosity of the fluid. To address this disparity, we now compute the speed of the interfacial wave when the half-space is in compression.

In order to minimize the number of parameters, we use the neo-Hookean solid, with W given by (2.23). We take a normal equi-biaxial pre-strain (λ 1 = λ, λ 2 = λ -1/2 , λ 3 = λ -1/2 ) and, following Bagno and Guz, we take the frequency ω to be real. A dimensional analysis of the resulting dispersion equation shows that (5.14) now depends on just three non-dimensional parameters: a measure of the pre-strain, λ; a measure of the dynamic viscosity of the fluid compared with the shear modulus of the solid, µ * ω/C; and the ratio of the densities, ρ * /ρ. Once these quantities are specified, the dispersion equation may be solved for the non-dimensional complex unknown x defined by

x := ρ γ 12 ω k = ρ Cλ 2 S -1 , (5.22)
where S is the (complex) scalar slowness. The dispersion equation can now be solved numerically for x, and the interfacial wave speed, normalized with respect to the transverse bulk shear wave speed in the deformed solid, is c = ℜ(1/x).

For Figure 2(a), we fix ρ * /ρ at 1.0 and we take in turn µ * ω/C = 0.2, 0.04, 0.02. The first choice (ρ * /ρ = 1.0, µ * ω/C = 0.2) is roughly that obtained for the blood/artery interface of the previous section with C 1 = C and C 2 = 0. We see clearly that as µ * ω/C decreases the wave speed decreases towards zero as λ tends to 0.444, the critical compression stretch for the solid/vacuum interface (the thick curve gives the solid/vacuum interface wave speed). We find that in the extension to moderate compression range, the solid/fluid interface wave speed is significantly lower than the solid/vacuum wave speed. In the strongly compressive range (as λ → 0.444), the speed plot dips towards zero, and dips further as µ * ω/C decreases, without ever reaching that value; the plot gets close to the plot for the solid/vacuum interface wave speed as λ decreases, but then crosses it, and the fluid/solid interface wave speed increases again. Note that we checked the validity of the solution at all compressive stretches using (5.18). For Figure 2(b), we fixed µ * ω/C at 1.0 and took ρ * /ρ = 0.1, 0.05, 0.01. In the extension to moderate compression range, the plot for the solid/fluid interface wave speed gets closer to that for the solid/vacuum interface wave speed as ρ * /ρ decreases. In the strongly compressive range, similar comments to those made for Figure 2(a) apply.

The conclusion is that as both µ * ω/C and ρ * /ρ tend to zero, the speed tends to zero when the stretch tends to the critical compression stretch of the solid/vacuum surface instability. This is to be expected because this double limit corresponds to the vanishing of the fluid. However, as emphasized earlier, the speed never reaches the zero limit and the bifurcation criterion is therefore never met. For instance, the typical values used by Bagno and co-workers [START_REF] Bagno | Elastic waves in prestressed bodies interacting with a fluid[END_REF] put µ * ω/C at about 0.0002 and ρ * /ρ as low as 0.1, giving c smaller than 10 -4 , but not zero. In effect a localized damped wave exists for the whole compressive range, with speed values starting below the solid/vacuum interface wave speed in the moderate compressive range, reaching a minimum (above the solid/vacuum interface wave speed) as λ approaches 0.444, and then rising rapidly to infinity as λ decreases below 0.444. This result suggests that when a neo-Hookean half-space is in contact with a viscous fluid it becomes completely stable. This, however, is an incorrect deduction, because it is based on special motions, for which the frequency ω is assumed real. Other motions might be unstable, as is illustrated below.

We now take the wave number k to be real (k > 0) and let the speed v be complex,

k > 0, v = v + + iv -, (5.23)
so that the motion is now proportional to e kv -t e ik(x 1 -v + t) . Clearly in this case, the conditions for a stable, localized wave, propagating in the x 1 > 0 direction, are

ℜ(v) 0, ℑ(v) 0, ℜ(s * ) > 0, ℜ(s 1 ) > 0, ℜ(s 2 ) > 0. ( 5.24)
(note that ℜ(s 1 ) > 0 is automatically satisfied in a neo-Hookean solid, because s 1 = 1).

Then, an analysis of the dispersion equation (5.14), in the case of a neo-Hookean solid with normal equi-biaxial pre-strain, reveals three nondimensional quantities: λ, µ * k/ √ ρC, and ρ * /ρ. Once they are specified, we solve the dispersion equation for the non-dimensional quantity x defined as

x := ρ γ 12 ω k = ρ Cλ 2 v.
(5.25)

The interfacial wave speed, normalized with respect to the transverse bulk shear wave speed in the deformed solid, is c = ℜ(x).

For Figure 3, we take ρ * /ρ = 1.0, and µ * k/ √ ρC = 0.2, 0.002 in turn, and we plot v against λ in the compressive range. We find that at a compressive stretch close to the critical compressive stretch of surface stability for the solid/vacuum interface, the normalized speed c drops to zero. From the figure it is not clear that the values of λ at this point are different for 0.2 and 0.002, but the zoom in Figure 3 shows, however, that the value of the compressive stretch at which c = 0 depends on the material parameters. For comparison, the curve corresponding to µ * k/ √ ρC = 1 is also shown. We emphasize that the situation at c = 0 does not correspond to a static solution of the equations of motion (which would be impossible in the fluid), but to a non-propagating damped motion, proportional to e kv -t e ikx 1 . Moreover, this situation does not correspond to an instability because at that point, and a bit beyond, the requirements (5.24) still hold. For instance, in the case ρ * /ρ = 1.0, µ * k/ √ ρC = 0.2, the speed drops to zero at λ ≃ 0.44539 but there still exist non-propagating, localized motions for λ > 0.42212. Beyond that stretch value, however, all solutions of the considered type grow unboundedly with time and/or with space, indicating instability. From this example we see that the viscous load stabilizes slightly the solid half-space, because it can now be compressed by an extra 2%, from 0.444 to 0.422.

Non-principal wave: neo-Hookean solid

In this section we consider that the solid half-space is a deformed neo-Hookean material, with strain-energy function given by Eq. (2.23). For this material, Flavin [START_REF] Flavin | Surface waves in pre-stressed Mooney material[END_REF] noticed that Rayleigh waves are plane-polarized for any direction of propagation in a principal plane; see [START_REF] Braun | Rayleigh waves in a prestressed neo-Hookean material[END_REF], [START_REF] Dowaikh | Surface waves propagating in a half-space of neo-Hookean elastic material subject to pure homogeneous strain[END_REF]. Chadwick and Jarvis [START_REF] Chadwick | Interfacial waves in a pre-strained neo-Hookean body I. Biaxial states of strain[END_REF] showed the same result for Stoneley waves. Indeed, we now show that the saggital motion is always decoupled from the 'anti-saggital' motion, for any θ and any type of (bulk or interfacial) inhomogeneous wave.

First, we note that for neo-Hookean bodies, the parameters γ ij and β ij have the simple expressions

γ ij = Cλ 2 i , β ij = C(λ 2 i + λ 2 j )/2. (6.1)
It follows that the matrices N 1 , N 2 , N 3 , given by Eq. (2.20), are greatly simplified. Next, we recall that the direction of propagation and the normal to the interface define the saggital plane, with unit normal [-s θ , 0, c θ ] T . Now consider the new unknown functions w i , t 2i (i = 1, 2, 3), defined by

w i = Ω ij vj , t 2i = Ω ij ŝ2j , (6.2)
where

Ω ij =   c θ 0 s θ 0 1 0 -s θ 0 c θ   . (6.3)
Some algebraic manipulations reveal that the equations of motion (2.17), written for w i , t 2i , decouple the 'anti-saggital' motion [w with roots s = -1, which is discarded because it does not lead to a decaying wave,

s 1 = 1, s 2 = ± [C(c 2 θ λ 2 1 + s 2 θ λ 2
3 ) -ρ]/(Cλ 2 2 ). (6.7)

The analysis leading to the derivation of the dispersion equation is by and large the same as that conducted for the principal wave in Section 5. The end result is that the dispersion equation is again (5.13), where now s 1 and s 2 are given by the expressions above, and γ 21 is replaced by Cλ 2 2 , but the other quantities remain unchanged.

As an illustration, we take the same numerical values as in Section 5.4, but with C 2 = 0 in (5.16). Hence the solid is neo-Hookean with C = 0.175 MPa and ρ = 1000 kg/m 3 . It is under the pre-stress σ 2 = 20 kPa, σ 3 = 0. The fluid has Newtonian viscosity µ * = 3.5 × 10 -3 Ns/m 2 , and mass density ρ * = 1050 kg/m 3 . The frequency is ω = 10 7 Hz. In turn, we take the solid to be compressed by 20% in the x 1 direction (so that λ 1 = 0.8 and then λ 2 ≃ 1.131, λ 3 ≃ 1.105), to be unstretched in the x 1 direction (so that λ 1 = 1, and λ 2 ≃ 1.014, λ 3 ≃ 0.986), and to be under an extension of 20% in the x 1 direction (so that λ 1 = 1.2, and λ 2 ≃ 0.929, λ 3 ≃ 0.897).

Figure 4 shows the dependence of the interfacial wave speed v + = ℜ(ω/k) on the angle of propagation θ. When the solid is almost unstrained, the speed hardly varies with the angle; when it is strained to ±20%, the induced anisotropy causes speed changes of more than ±25% in some directions. The figure also shows clearly that the wave travels at its fastest in the direction of greatest stretch, and at its slowest in the direction of greatest compression, indicating that it could be used for the acoustic determination of these directions.

15 -Figure 1 :

 151 Figure 1: Depth profiles as functions of x 2 /λ 0 of the velocity components of the acoustic wave: x 2 (normal) component -continuous curve; x 1 (tangential) component -dotted curve. The zooms show continuity of these components and discontinuity of their derivatives.

Figure 2 :

 2 Figure 2: Normalized wave speed for a neo-Hookean solid under normal equibiaxial pre-strain and a viscous fluid: (a) ρ * /ρ = 1.0 and µ * ω/C = 0.2, 0.04, 0.02; (b) ρ * /ρ = 1.0, 0.4, 0.2 and µ * ω/C = 0.1. The thick curves are for an unloaded solid half-space.

Figure 3 :

 3 Figure 3: Normalized wave speed for a neo-Hookean solid under normal equi-biaxial pre-strain and a viscous fluid with ρ * /ρ = 1.0 and µ * k/ √ ρC = 0.2, 0.002, 1 (indicated by the arrows) for the compressive range 0.4 < λ < 1.The zoom shows the details for µ * k/ √ ρC = 0.2, 0.002 for 0.442 < λ < 0.445 range. The thick curve corresponds to the solid half-space with no fluid loading.

Figure 4 :

 4 Figure 4: Influence of propagation angle θ on the speed v + of a wave at the interface between a viscous fluid and a deformed neo-Hookean solid (thin curves). The thick curves represent the wave speed in the absence of fluid loading.

  3 , t 23 ] from its saggital counterpart. For the latter we find [w ′ 1 , w ′ 2 , t ′ 21 , t ′ 22 ] T = ikN [w 1 , w 2 , t 21 , t 22 ] T , A search for partial-mode solutions in the form [w 1 , w 2 , t 21 , t 22 ] T = ζe -ksx 2 , where ζ is a 4-component constant vector, leads to an eigenvalue problem. The associated characteristic equation is the propagation condition (s 2 -1)[Cλ 2 2 s 2 -C(c 2

			(6.4)
	with N in the form (5.1), where now
	η = C(c 2 θ λ 2 1 + s 2 θ λ 2 3 + 3λ 2 2 ) -2σ 2 , ν = C(c 2 θ λ 2 1 + s 2 θ λ 2 3 ) + Cλ 2 2 (1 -σ 2 ) 2 ,	γ 21 = Cλ 2 2 , σ 2 = σ 2 /(Cλ 2 2 ).	(6.5)
	θ λ 2 1 + s 2 θ λ 2 3 ) + ρ] = 0,	(6.6)


Charles-Edouard Bréhier
email: brehier@math.univ-lyon1.fr

Tony Lelièvre
email: lelievre@cermics.enpc.fr

On a new class of score functions to estimate tail probabilities of some stochastic processes with Adaptive Multilevel Splitting

We investigate the application of the Adaptive Multilevel Splitting algorithm for the estimation of tail probabilities of solutions of Stochastic Differential Equations evaluated at a given time, and of associated temporal averages.

We introduce a new, very general and effective family of score functions which is designed for these problems. We illustrate its behavior on a series of numerical experiments. In particular, we demonstrate how it can be used to estimate large deviation rate functionals for the longtime limit of temporal averages.

I. INTRODUCTION

Fast and accurate estimation of rare event probabilities, and the effective simulation of these events, is a challenging computational issue, which appears in many fields of science and engineering. Since rare events are often the ones that matter in complex systems, designing efficient and easily implementable algorithms is a crucial question which has been the subject of many studies in the recent years.

Since the pionnering works on Monte-Carlo methods, several classes of algorithms have been developed, see for instance the monographs [START_REF] Asmussen | Stochastic simulation: algorithms and analysis[END_REF][START_REF] Bucklew | Introduction to rare event simulation[END_REF][START_REF] Rubino | Introduction to rare event simulation[END_REF] . The most popular strategies are importance sampling and splitting. On the one hand, importance sampling consists in changing the probability distribution, such that under the new probability distribution the events of interest are not rare anymore. Appropriate reweighting then yields consistent estimators. This strategy has for instance been applied recently to simulate rare events in climate models [START_REF] Ragone | Computation of extreme heat waves in climate models using a large deviation algorithm[END_REF] . On the other hand, splitting techniques, consist in writing the rare event probability as a product of conditional probabilities which are simpler to estimate, and in using interacting particle systems in order to estimate these conditional probabilities.

In this manuscript, a class of splitting algorithms is considered. Splitting techniques have been introduced in the 1950s [START_REF] Kahn | Estimation of particle transmission by random sampling[END_REF] , and have been studied extensively in the last two decades [START_REF] Cérou | Sequential Monte Carlo for rare event estimation[END_REF][START_REF] Dean | Splitting for rare event simulation: a large deviation approach to design and analysis[END_REF][START_REF] Garvels | On the importance function in splitting simulation[END_REF][START_REF] Glasserman | Multilevel splitting for estimating rare event probabilities[END_REF] . Many variants have appeared in the literature: Generalized multilevel splitting [START_REF] Botev | Efficient Monte Carlo simulation via the generalized splitting method[END_REF][START_REF] Botev | Static network reliability estimation via generalized splitting[END_REF] , RESTART [START_REF] Villén-Altamirano | RESTART: A method for accelerating rare events simulations[END_REF][START_REF] Villén-Altamirano | RESTART: a straightforward method for fast simulation of rare events[END_REF] , Subset simulation [START_REF] Au | Estimation of small failure probabilities in high dimensions by subset simulation[END_REF] , Nested sampling [START_REF] Skilling | Nested sampling for general Bayesian computation[END_REF][START_REF] Skilling | Nested sampling for Bayesian computations[END_REF] , Reversible shaking transformations with interacting particle systems [START_REF] Agarwal | Study of new rare event simulation schemes and their application to extreme scenario generation[END_REF][START_REF] Gobet | Rare event simulation using reversible shaking transformations[END_REF] , genealogical particle analysis [START_REF] Del Moral | Genealogical particle analysis of rare events[END_REF][START_REF] Wouters | Rare event computation in deterministic chaotic systems using genealogical particle analysis[END_REF] , etc... The Adaptive Multilevel Splitting (AMS) algorithm [START_REF] Cérou | Adaptive multilevel splitting for rare event analysis[END_REF] is designed to estimate rare event probabilities of the type P(τ B < τ A ), where τ A and τ B are stopping times associated with a Markov process X, typically the entrance times of X in regions A and B of the state space. In many applications, A and B are metastable states for the process. The algorithm is based on selection and mutation mechanisms, which leads to the evolution of a system of interacting replicas. The selection is performed using a score function, which is often referred to as a reaction coordinate when dealing with metastable systems.

The objective of this article is to design and test new score functions, using the AMS strategy, to estimate probabilities of the type P(Φ(X T ) > a) or P 1 T T 0 φ(X t )dt > a , where a is a threshold, Φ and φ are real-valued functions, and X t 0≤t≤T is a Markov process. In fact, as will be explained below, the probabilities of interest can be rewritten as P(τ B < τ A ), associated with an auxiliary Markov process. Our main contribution is the identification of appropriate score functions related to this interpretation, and which return a non-zero value for the estimator of the probability of interest. By using then an AMS algorithm which fits in the Generalized Adaptive Multilevel Splitting framework [START_REF] Bréhier | Unbiasedness of some generalized adaptive multilevel splitting algorithms[END_REF] , one can construct unbiased estimators of the probability (and possibly of other quantities of interest).

The efficiency of the approach is investigated with numerical experiments, using several test cases taken from the literature on rare events. First, validation is performed on onedimensional Gaussian models (Brownian Motion [START_REF] Vanden-Eijnden | Rare event simulation of small noise diffusions[END_REF] , Ornstein-Uhlenbeck process [START_REF] Wouters | Rare event computation in deterministic chaotic systems using genealogical particle analysis[END_REF] ). More complex test cases then illustrate the efficiency of the approach and of the new score functions introduced in this article: drifted Brownian Motion, three-dimensional Lorenz model [START_REF] Beck | Rare-event simulation[END_REF][START_REF] Lorenz | Deterministic nonperiodic flow[END_REF] . Estimations of probabilities depending on temporal averages 1 T T 0 φ(X s )ds are also considered for two models [START_REF] Ferré | Adaptive sampling of large deviations[END_REF] : the one-dimensional Ornstein-Uhlenbeck process and a driven periodic diffusion. In these examples, values of large deviations rate functionals for the longtime limit T → ∞ are estimated.

In the last decade, many works have been devoted to the analysis and applications of AMS algorithms. A series of work has been devoted to the analysis of the so-called ideal case [START_REF] Bréhier | Large deviations principle for the adaptive multilevel splitting algorithm in an idealized setting[END_REF][START_REF] Bréhier | Central limit theorem for adaptive multilevel splitting estimators in an idealized setting[END_REF][START_REF] Bréhier | Analysis of adaptive multilevel splitting algorithms in an idealized case[END_REF][START_REF] Guyader | Simulation and estimation of extreme quantiles and extreme probabilities[END_REF] , namely when the AMS algorithm is applied with the optimal score function (namely the so-called committor function). In practice, this optimal score function is unkown. Beyond the ideal case, consistency [START_REF] Bréhier | Unbiasedness of some generalized adaptive multilevel splitting algorithms[END_REF] (unbiasedness of the small probability estimator) and efficiency [START_REF] Cérou | On the asymptotic normality of adaptive multilevel splitting[END_REF] (variance of the small probability estimator) have been studied. Moreover, the adaption of the original algorithm to the discrete-in-time setting has been studied in details in [START_REF] Bréhier | Unbiasedness of some generalized adaptive multilevel splitting algorithms[END_REF] . It can be used to compute transition times between metastable states [START_REF] Cérou | A multiple replica approach to simulate reactive trajectories[END_REF] , return times [START_REF] Lestang | Computing return times or return periods with rare event algorithms[END_REF] , or other observables associated with the rare event of interest [START_REF] Louvin | Development of an adaptive variance reduction technique for Monte Carlo particle transport[END_REF] . The AMS algorithm has been successfully applied in many contexts: the Allen-Cahn stochastic partial differential equation [START_REF] Bréhier | CEMRACS 2013-modelling and simulation of complex systems: stochastic and deterministic approaches[END_REF][START_REF] Rolland | Computing transition rates for the 1-D stochastic Ginzburg-Landau-Allen-Cahn equation for finite-amplitude noise with a rare event algorithm[END_REF] , the simulation of Bose-Einstein condensates [START_REF] Poncet | Méthodes numériques pour la simulation d'équations aux dérivées partielles stochastiques non-linéaires en condensation de Bose-Einstein[END_REF] , molecular dynamics and computational chemistry [START_REF] Aristoff | CEMRACS 2013-modelling and simulation of complex systems: stochastic and deterministic approaches[END_REF][START_REF] Cérou | A multiple replica approach to simulate reactive trajectories[END_REF][START_REF] Lopes | Analysis of the adaptive multilevel splitting method with the alanine di-peptide's isomerization[END_REF][START_REF] Teo | Adaptive multilevel splitting method for molecular dynamics calculation of benzamidine-trypsin dissociation time[END_REF] , nuclear physics [START_REF] Louvin | Adaptive multilevel splitting for Monte Carlo particle transport[END_REF][START_REF] Louvin | Development of an adaptive variance reduction technique for Monte Carlo particle transport[END_REF][START_REF] Louvin | Three-dimensional neutron streaming calculations using adaptive multilevel splitting[END_REF] and turbulence [START_REF] Bouchet | A rare event algorithm links transitions in turbulent flows with activated nucleations[END_REF][START_REF] Rolland | Extremely rare collapse and build-up of turbulence in stochastic models of transitional wall flows[END_REF] , for example. This article is organized as follows. Section II presents the precise mathematical setting, in particular the rare event probability of interest is defined by (7). A general formulation of the AMS algorithm designed to estimate this quantity is provided in Section III, in particular see Section III B for the full algorithmic description. Examples of appropriate score functions are discussed in Section IV. To overcome the limitations of a vanilla strategy, Section IV A, our main contribution is the construction of the score functions presented in Section IV B.

Finally, numerical experiments are reported in Section V.

II. SETTING

We consider stochastic processes, with values in R d , in dimension d ∈ N, which are solutions of Stochastic Differential Equations (SDEs) of the type: for 0 ≤ t 0 ≤ t ≤ T and

x 0 ∈ R d , dX t 0 ,x 0 t = f (t, X t 0 ,x 0 t )dt + σ(t, X t 0 ,x 0 t )dW (t) (1)
where

X t 0 ,x 0 t ∈ R d , with initial condition, X t 0 ,x 0 t 0 = x 0 . (2)
The noise W (t) t≥0 is given by a standard Wiener process with values in R D , for some

D ∈ N. The coefficients f : [0, T ] × R d → R d and σ : [0, T ] × R d → R
d×D are assumed to be sufficiently smooth to ensure global well-posedness of the SDE.

In this work, two types of rare events associated with (X t 0 ,x 0 t ) 0≤t≤T are considered. Let a ∈ R denote a threshold, and let Φ, φ : R d → R be two measurable functions. First, we are interested in tail probabilities for the random variable Φ(X t 0 ,x 0 T ), namely in

P Φ(X t 0 ,x 0 T ) > a . (3)
Second, we are interested in tail probabilities for temporal averages, defined as

P 1 T -t 0 T t 0 φ(X t 0 ,x 0 t )dt > a . (4)
We will investigate numerically the performance of AMS estimators for both (3) and ( 4)

on various examples. In particular, we will consider the regime T → ∞ for (4) in order to estimate large deviation rate functionals.

Notice that the case of temporal averages (4) can be rewritten in the form of (3). Indeed, the probability (4) may be written as (3) for the auxiliary process defined by Xt 0 ,x

0 t = X t 0 ,x 0 t , Y t 0 ,x 0 t
, where Y t 0 ,x 0 t = φ(x 0 ), and for t > t 0

Y t 0 ,x 0 t = 1 t -t 0 t t 0 φ(X t 0 ,x 0 s )ds,
and with Φ(x, y) = y. The process Y t 0 ,x 0 t t 0 ≤t≤T is solution of the following ODE,

dY t 0 ,x 0 t = 1 t -t 0 φ(X t 0 ,x 0 t ) -Y t 0 ,x 0 t , t > t 0 , Y t 0 ,x 0 t 0 = φ(x 0 ),
coupled with the SDE for the diffusion process X t 0 ,x 0 t t 0 ≤t≤T . This trick will be used for our numerical experiments below. Therefore, in the following, we present the AMS algorithm and discuss its theoretical properties only for (3).

For future purposes, observe that the target probability (3) may be written as u a (t 0 , x 0 ), where

u a (t, x) = P Φ(X t,x T ) > a (5)
is the solution (under appropriate regularity assumptions) of the backward Kolmogorov

equation      ∂ua(t,x) ∂t + L t u a (t, x) = 0 for t ∈ [0, T ] and x ∈ R d , u a (T, x) = 1 Φ(x)>a for x ∈ R d , (6)
where the infinitesimal generator L t is defined by: for all test functions ϕ,

L t ϕ(x) = f (t, x) • ∇ϕ(x) + 1 2 σ(t, x)σ(t, x) : ∇ 2 ϕ(x).
Approximating the solutions of PDEs of this type using deterministic methods is in general possible only when the dimension d is small. Instead, Monte Carlo methods may be used. However, naive Monte Carlo algorithms are not efficient in the rare event regime, e.g. when a → ∞ or when the diffusion coefficient is of the type σ = √ σ and → 0.

In practice, discrete-time approximations are implemented. Let ∆t > 0 denote the timestep size of the integrator (for instance the standard Euler-Maruyama method), with T = N ∆t and t 0 = n 0 ∆t, where n 0 ∈ N 0 , N ∈ N, n 0 ≤ N . With a slight abuse of notation, let us denote the discrete-time process obtained after discretization of (1) by X n 0 ,x 0 n n 0 ≤n≤N . The time-discrete counterpart of (3) is then

P Φ(X n 0 ,x 0 N ) > a . (7)
The algorithms presented below are used to estimate probabilities of the type (7).

Remark 1. It is assumed that the initial condition is deterministic: X n 0 ,x 0 n 0 = x 0 . The adaptation of the algorithms presented below to the case of a random initial condition is straightforward, by simply using the Markov property:

P Φ(X N ) > a = P Φ(X n 0 ,x 0 N ) > a dµ 0 (x 0 )
where µ 0 denotes the law of X n 0 .

III. GENERAL FORMULATION OF THE ADAPTIVE MULTILEVEL SPLITTING ALGORITHM A. Context

The goal is to estimate the probability p given by (7), in the regime where p is small, which is for example the case when a is large.

It is convenient to introduce an auxiliary process Z n n 0 ≤n≤N , such that Z n = n∆t, X n 0 ,x 0 n . Indeed, let

A = {(T, x); Φ(x) ≤ a} , B = {(T, x); Φ(x) > a} ,
and define the associated stopping times

τ A = inf {n ≥ n 0 , n ∈ N 0 ; Z n ∈ A} , τ B = inf {n ≥ n 0 , n ∈ N 0 ; Z n ∈ B} .
Then the probability p given by ( 7) can be rewritten as

p = P Φ(X n 0 ,x 0 N ) > a = P(τ B < τ A ). (8)
We are then in position to build algorithms which fit in the Generalized Adaptive Multilevel Splitting framework developed in [START_REF] Bréhier | Unbiasedness of some generalized adaptive multilevel splitting algorithms[END_REF] , which ensure that the obtained estimators of the probability (8) are unbiased.

For that, a score function, or reaction coordinate, ξ, needs to be given. Following the interpretation above, ξ may depend on z = (n∆t, x).

To run the algorithm and define simple unbiased estimators of p, only one requirement is imposed on the function ξ: there exists ξ max such that

B ⊂ {z; ξ(z) > ξ max } ,
which in the context of this article is rephrased as

Φ(x) > a =⇒ ξ(T, x) > ξ max . (9)
The principle of splitting algorithm is then to write

P(τ B < τ A ) = P(τ ζ 1 < τ A )P(τ ζ 2 < τ A |τ ζ 1 < τ A )P(τ ζ 3 < τ A |τ ζ 2 < τ A ) . . . P(τ B < τ A |τ ξmax < τ A )
for an increasing sequence of levels (ζ q ) q≥1 , where

τ ζ = inf{n ≥ n 0 ; ξ(Z n ) > ζ}.
If the levels are well chosen, then the successive conditional probabilities P(τ ζ q+1 < τ A |τ ζq < τ A ) are easy to compute. The principle of the adaptive multilevel splitting algorithm 18 is to choose the levels adaptively, so that the successive conditional probabilities P(τ ζ q+1 < τ A |τ ζq < τ A ) are constant and fixed. The levels constructed in the algorithm are then random.

B. The Adaptive Multilevel Splitting algorithm

Before giving the detailed algorithm, let us roughly explain the main steps (we also refer to [START_REF] Bréhier | Unbiasedness of some generalized adaptive multilevel splitting algorithms[END_REF] for more details and intuition on the algorithm). In the initialization, one samples n rep trajectories following (1)-( 2) and compute the score of each trajectory, namely the maximum of ξ attained along the path. Then the algorithm proceeds as follows: one discards the trajectory which has the smallest score and in order to keep the number of trajectories fixed, a new one is created by choosing one of the remaining trajectories at random, copying it up to the score of the killed trajectory, and sampling the end of trajectory independently from the past. This is called the partial resampling. One thus obtains a new ensemble of n rep trajectories on which one can iterate by again discarding the the trajectory which has the smallest score. As the iteration goes, one thus obtains trajectories with largest and largest scores, and an estimate of the probability of interest is obtained as (1

-1/n rep ) Q iter P (τ B < τ A |τ ξmax < τ A ) (notice that (1 -1/n rep
) is an estimate of the conditional probability to reach level ζ q+1 conditionally to the fact that level ζ q has been reached), where Q iter is the number of iterations required to reach the maximum level ξ max . In practice,

P (τ B < τ A |τ ξmax < τ A )
is estimated by the proportion of trajectories which reach B before A at the last iteration of the algorithm, namely when all the trajectories satisfy τ ξmax < τ A .

Actually, the algorithm has to be adapted in order to take into account situations when more than one particle has the smallest score, which happens with non zero probability for Markov chains. Let us now give the details of the AMS algorithm.

To simplify notation, in the sequel, the initial condition x 0 and the time n 0 are omitted in the notation of the replicas.

a. Input

• n rep ∈ N, the number of replicas,

• a score function z = (n∆t, x) → ξ(n∆t, x) ∈ R and a stopping level ξ max ∈ R such that ( 9) is satisfied.

b. Initialization

• Sample n rep independent realizations of the Markov process

X j = X j m n 0 ≤m≤N , 1 ≤ j ≤ n rep
following the dynamics (1)-( 2).

• Compute the score of each replica, M j = max

n 0 ≤m≤N ξ(m∆t, X j m ). • Compute the level Z = min 1≤j≤nrep M j . • Define K = {j ∈ {1, . . . , n rep } ; M j = Z}. • Set q = 0, p = 1, B = 1. c. Stopping criterion If Z ≥ ξ max or card(K) = n rep , then set B = 0. d. While B = 1 • Update -q ← q + 1 and p ← p • 1 -card(K) n .
• Splitting -Reindex the replicas, such that

     M j = Z if j ∈ {1, . . . , card(K)} M j > Z if j ∈ {card(K) + 1, . . . , n rep } .
-For replicas with index j ∈ {1, . . . , card(K)}, sample labels 1 , . . . , card(K) , independently and uniformly in {card(K) + 1, . . . , n rep }.

• Partial resampling -Remove the replicas with label j ∈ {1, . . . , card(K)}.

-

For j ∈ {1, . . . , card(K)}, define m j = inf m ∈ {n 0 , . . . , N } ; ξ(X j m ) > Z .
-For m ∈ {n 0 , . . . , m j }, set X j m = X j m .

-Sample a new trajectory X j m m j ≤m≤N with the Markov dynamics (1) driven by independent realizations of the Brownian motion.

• Level computation -Compute the scores M j = max n 0 ≤m≤N ξ(m∆t, X j m ).
-Compute the level Z = min 1≤j≤nrep M j .

-Define the set K = {j ∈ {1, . . . , n rep } ; M j = Z}.

e. Stopping criterion

If Z ≥ ξ max or card(K) = n rep , then set B = 0.
f. End while g. Update:

p ← p 1 nrep nrep j=1 1 Φ(X j N )≥a . h. Output: p and Q iter = q.
Remark 2. We presented the algorithm in its simplest form. There are many variants, see [START_REF] Bréhier | Unbiasedness of some generalized adaptive multilevel splitting algorithms[END_REF] . For example, the killing level Z can be defined as Z = M (k) where M [START_REF] Agarwal | Study of new rare event simulation schemes and their application to extreme scenario generation[END_REF] ≤ M (2) ≤ . . . ≤ M (nrep) denotes an increasing relabelling of the scores (M j ) 1≤j≤nrep (order statistics).

C. Consistency result

Let p and Q iter be the outputs of the realization of the algorithm above. We quote the following result [START_REF] Bréhier | Unbiasedness of some generalized adaptive multilevel splitting algorithms[END_REF] , which states that the output p of the algorithm described in Section III B above is an unbiased estimator of the probability given by (7).

Theorem 1. Let ξ be a score function and ξ max ∈ R be such that (9) is satisfied. Let n rep ∈ N be a given number of replicas. Assume that almost surely the algorithm stops after a finite number of iterations:

Q iter < ∞ almost surely.
Then p is an unbiased estimator of the probability p given by (8):

E[p] = P Φ(X n 0 ,x 0 N ) > a . Note that if ψ : R d → R is a function with support included in {x; Φ(x) > a}, i.e. ψ(x) = 0 if Φ(x) ≤ a, then an unbiased estimator of E[ψ(X n 0 ,x 0 N
)] is given by replacing the final update in the algorithm above (see step g.), by

p ← p 1 n rep nrep j=1 ψ(X j N ).
The unbiasedness property is crucial in practice for the following two reasons. First, it is very easy to parallelize the estimation of rare events using this property. Indeed, since the estimator is unbiased whatever the value of n rep , to get a convergent estimation, one has simply to fix n rep to a value which enables the computation of p on a single CPU, and then to sample M independent realizations of this p, run in parallel. In the large M limit, one obtains a convergent estimator of the quantity of interest by simply considering the average of the realizations of p. Second, the practical interest of Theorem 1 is that since E(p) is the same whatever the choice of the numerical parameters (namely n rep and ξ), one can compare the results obtained with different choices to get confidence in the result. For example, one can consider the confidence intervals obtained with M independent realizations of p for two different choices of ξ, and check whether these confidence intervals overlap or not.

Remark 3. In the algorithm described in Section III B, the set K defined in the initialization and in the level computation steps may have a cardinal strictly larger than 1, even if the level Z is defined as the minimum of the scores over the replicas. This simply means than more than one replica has a score which is the smallest among the replicas. In the discrete-time setting (namely for Markov chains), this happens with non zero probability, and it requires an appropriate modification of the original AMS algorithm, as described above, see [START_REF] Bréhier | Unbiasedness of some generalized adaptive multilevel splitting algorithms[END_REF] for more details.

Notice that in particular, there is a possibility that the algorithm stops if card(K) = n rep , in which case there is an extinction of the system of replicas.

IV. CHOICES OF THE SCORE FUNCTION

Let us now describe the various score functions we will consider in order to estimate (7).

A. Vanilla score function and limitations

The simplest choice consists in choosing the score function as given by

ξ std (n∆t, x) = Φ(x), (10)
with ξ max = a. In this case, the score function does not depend on the time variable.

If the conditional probability

q = P Φ(X n 0 ,x 0 N ) > a max n 0 ≤n≤N Φ(X n 0 ,x 0 n ) > a
is small, the performance of this vanilla strategy may be poor. Indeed, only a small proportion of the replicas have a non-zero contribution to the value of the estimator of the probability (7). It may even happen that all the replicas satisfy max

n 0 ≤n≤N Φ(X n 0 ,x 0 n ) > a but that none of them satisfies Φ(X n 0 ,x 0 N ) > a.
In that situation, the algorithm returns p = 0 to estimate p > 0. One of the goals of this work is to construct score functions which circumvent that issue: with the score functions introduced below, almost surely p = 0.

B. Time-dependent score functions

As discussed above, it is natural to design score functions ξ, which satisfy the following condition:

{Φ(X n 0 ,x 0 N ) > a} = max n 0 ≤n≤T ξ(n∆t, X n 0 ,x 0 n ) > 1 . ( 11
)
The choice of the value ξ max = 1 on the right-hand side above is arbitrary, but no generality is lost. Indeed, if a score function ξ satisfying ( 11) is used in the AMS algorithm above, with ξ max = 1, at the last update, the ratio

1 nrep nrep j=1 1 Φ(X j N
)≥a is identically equal to 1, since all replicas satisfy max

n 0 ≤n≤T ξ(n∆t, X n 0 ,x 0 n ) > 1.
In particular, by construction p = 0 (provided

Q iter < ∞).
As will be seen below, in practice it is more natural to identify functions ξ taking values in (-∞, 1], which satisfy the condition

{Φ(X n 0 ,x 0 N ) > a} = max n 0 ≤n≤T ξ(n∆t, X n 0 ,x 0 n ) = 1 , (12)
instead of (11). To justify the use of the algorithm in this case, observe that ξ(t, x) = ξ(t, x)+1 ξ(t,x)=1 then satisfies (11). In addition, when running the algorithm, choosing either ξ or ξ exactly yields the same result. We are thus in the setting where the unbiasedness result Theorem 1 applies. The score functions presented below will satisfy (12) instead of (11).

One of the novelties of this article is the introduction of the following score function:

ξ new (n∆t, x) = Φ(x) -a 1 Φ(x)≤a + n∆t N ∆t 1 Φ(x)>a . (13)
Observe that ξ new takes values in (-∞, 1], and that ξ new (n∆t, x) = 1 if and only if n = N and Φ(x) ≥ a. Thus the condition ( 12) is satisfied. We refer to Figure 1 for a schematic representation of this score function. Note that the score function defined by ( 13) only depends on the function Φ, on the threshold a, and on the final time T = N ∆t. It may thus be applied in any situation, but in some cases better score functions may be built upon using more information on the dynamics. The practical implementation is very simple.

Let us explain how the AMS algorithm proceeds when used with the score function (13).

Observe that max

n 0 ≤n≤N ξ new (n∆t, X n 0 ,x 0 n ) ≥ 0 = max n 0 ≤n≤N Φ(X n 0 ,x 0 n ) > a .
The first iterations of the algorithm, up to reaching level 0, are thus devoted to construct n rep replicas which satisfy the weaker condition max

n 0 ≤n≤N Φ(X n 0 ,x 0 n ) > a .
In other words, if the stopping level ξ max in the algorithm is set equal to 0 instead of 1, one thus recovers the vanilla AMS algorithm described above, applied with the score function ξ(t, x) = Φ(x) (independent of time t).

Compared with the vanilla score function, the AMS algorithm the new score function does not stop when max

n 0 ≤n≤N Φ(X n 0 ,x 0 n ) > a .
In terms of splitting, observe that this consists in writing

P Φ(X n 0 ,x 0 N ) > a = P Φ(X n 0 ,x 0 N ) > a | max n 0 ≤n≤N Φ(X n 0 ,x 0 n ) > a P max n 0 ≤n≤N Φ(X n 0 ,x 0 n ) > a ,
and the remaining effort consists in estimating the conditional probability above.

More generally, observe that for every n 1 ∈ {n 0 , . . . , N },

max n 0 ≤n≤N ξ new (n∆t, X n 0 ,x 0 n ) ≥ n 1 N = max n 1 ≤n≤N Φ(X n 0 ,x 0 n ) > a .
The construction of the score function ( 13) is associated to the following family of nested events:

max

n 0 ≤n≤N ξ new (n∆t, X n 0 ,x 0 n ) ≥ 1 ⊂ . . . ⊂ max n 0 ≤n≤N ξ new (n∆t, X n 0 ,x 0 n ) ≥ n 1 N ⊂ . . . ⊂ max n 0 ≤n≤N ξ new (n∆t, X n 0 ,x 0 n ) ≥ 0 ,
which equivalently may be rewritten as

{Φ(X n 0 ,x 0 N ) > a} ⊂ . . . ⊂ max n 1 ≤n≤N Φ(X n 0 ,x 0 n ) > a ⊂ . . . ⊂ max n 0 ≤n≤N Φ(X n 0 ,x 0 n ) > a .
In the selection procedure, the intervals [n 1 ∆t, N ∆t] are iteratively reduced (by increasing the left end point of the interval), until they ultimately contain only the point N ∆t (the right end point of the interval which remains fixed).

To conclude, we mention that the construction given by ( 13) can be generalized as follows.

Let a : [0, T ] → R be a non-decreasing function, such that a(T ) = a. Define

ξ a new (n∆t, x) = Φ(x) -a(n∆t) 1 Φ(x)≤a(n∆t) + n∆t N ∆t 1 Φ(x)>a(n∆t) . (14)
The score function ξ new defined by ( 13) is a particular case of ( 14), with a(t) = a. Optimizing the choice of the function a may help improve the efficiency of the algorithm. Notice that condition ( 12) is satisfied with ξ = ξ a new . We refer to Figure 2 for a schematic representation of this score function.

C. The optimal score function: the committor function

For the general setting presented in Section III A where ones want to estimate P(τ B < τ A ), the committor function is defined as z → P z (τ B < τ A ), where the upperscript z refers to the initial condition of the process Z. In 17 , it is shown that, in a continous-time setting, the asymptotic variance (as the number of replicas n rep goes to infinity) of AMS algorithm is minimized when using the committor function as the score function. It is thus interesting to look at what the committor function looks like in our context.

In our context, the committor function is given by For the discussion, it is more convenient to consider the continuous-time version ξ com (t, x) = P Φ(X t,x T ) > a , for t ∈ [0, T ], that we still denote ξ com with a slight abuse of notation. Recall that u a (t, x) = P Φ(X t,x T ) > a = ξ com (t, x) satisfies the Kolmogorov backward equation ( 6), as explained in Section II.

ξ com (n∆t, x) = P Φ(X n,x N ) > a . (15)
As mentioned above, the asymptotic variance (as the number of replicas n rep goes to infinity) of AMS algorithm is minimized when using the committor function as the score function [START_REF] Cérou | On the asymptotic normality of adaptive multilevel splitting[END_REF] . The asymptotic variance is then -p 2 log(p) nrep , where p is the probability which is estimated. The analysis of the AMS algorithm in the ideal case, i.e. when using the committor function as a the score function, has been performed in many works [START_REF] Bréhier | Large deviations principle for the adaptive multilevel splitting algorithm in an idealized setting[END_REF][START_REF] Bréhier | Central limit theorem for adaptive multilevel splitting estimators in an idealized setting[END_REF][START_REF] Bréhier | Analysis of adaptive multilevel splitting algorithms in an idealized case[END_REF][START_REF] Guyader | Simulation and estimation of extreme quantiles and extreme probabilities[END_REF] .

Of course, in practice, the committor function is unknown and the asymptotic variance depends on the chosen score function. It has been proved [START_REF] Cérou | On the asymptotic normality of adaptive multilevel splitting[END_REF] that the asymptotic variance is always bounded from above by 2p (1-p) nrep , for any choice of the score function, where we recall that the asymptotic variance of the vanilla Monte-Carlo method is p (1-p) nrep . This can be seen as a sign of the robustness of the AMS approach to estimate rare event probability (contrary to importance sampling method which may result in a dramatic increase of the asymptotic variance compared with the vanilla Monte-Carlo method).

For simple Gaussian models, namely when X is a Brownian Motion, an Ornstein-Uhlenkeck process, or a drifted Brownian Motion, it is possible to compute analytically the committor function. This is useful to validate algorithms on test cases, as will be illustrated in Section V. Figure 3 represents the level lines of the committor function for a one-dimensional Brownian Motion (with T = 1 and a = 1). In that case,

ξ com (t, x) = 1 -F a -x √ T -t ,
where F is the cumulative distribution function of the standard Gaussian distribution N (0, 1), to be compared to the level sets of ξ new and ξ a new on Figures 1 and2. This form leads to define other families of appropriate score functions:

ξ(t, x) = 1 -F φ(t, x)
where φ(t, x) → t→∞ (-∞)1 Φ(x)>a +(+∞)1 Φ(x)<a . But the efficiency depends a lot on the choice of φ. In practice, we did not observe much gain in our numerical experiments, compared to the score function ξ new introduced in the previous section.

Let us mention that various techniques have been proposed in order to approximate the committor function, in particular in the context of importance sampling techniques for rare events, since the committor function also gives the optimal change of measure. If diffusions with vanishing noise are considered [START_REF] Dupuis | Importance sampling for multiscale diffusions[END_REF][START_REF] Dupuis | Escaping from an attractor: Importance sampling and rest points i[END_REF][START_REF] Vanden-Eijnden | Rare event simulation of small noise diffusions[END_REF] , solutions of associated Hamilton-Jacobi equations are good candidates to estimate the committor function. See also [START_REF] Schütte | Optimal control of molecular dynamics using markov state models[END_REF] for approximations based on coarse-grained models. Whether such constructions are possible when considering temporal averages, instead of the terminal value of the process is unclear.

V. NUMERICAL SIMULATIONS

Let p denote the rare event probability of interest. An estimator of p is calculated as the empirical average over independent realizations of the AMS algorithm, given a choice of score function ξ. The main objective of this section is to investigate the behavior of the algorithm when choosing ξ = ξ new given by (13). A comparison with the vanilla score function ξ = ξ std given by ( 10 In some of the numerical experiments below, the conditional probability

q = P Φ(X N ) > a max 0≤n≤N Φ(X n ) > a , (16)
is also estimated by q = p pmax , where pmax = 1 M M m=1 pmax,m is the estimator of the probability

p max = P max 0≤n≤N Φ(X n ) > a , (17)
which is estimated using the vanilla score function ξ = ξ std .

A. Validation using two Gaussian models

In this section, we validate the AMS alogorithm with various score functions on simple models for which the probability of the rare event is known with arbitrary precision.

Brownian Motion

We follow here numerical experiments from [START_REF] Vanden-Eijnden | Rare event simulation of small noise diffusions[END_REF] . Let d = 1, and consider the diffusion process given by

dX(t) = 2β -1 dW (t), X(0) = 0.1,
where W (t) t≥0 is a standard real-valued Wiener process.

The dynamics is discretized using the explicit Euler-Maruyama method, with time-step size ∆t = 10 -3 (notice that the numerical scheme gives here the exact solution):

X n+1 = X n + 2β -1 ∆tζ n ,
with X 0 = X(0) = 0.1, where ζ n 0≤n≤N are independent standard Gaussian random variables.

The goal is to estimate the probability

p = P |X N | > 1 .
This corresponds with the choice Φ(x) = |x|, a = 1, T = 1 so that N = T /∆t = 10 3 .

Since X(t) = X(0) + W (t), the law of X(t) is a Gaussian distribution, and the value of

P(|X(1)| ≥ 1
) can be computed exactly in terms of the cumulative distribution function of the standard Gaussian distribution.

Two numerical experiments are reported below, using ξ = ξ new . First, in Table I, the number of replicas is set equal to n rep = 10 2 , and the empirical average is computed over These numerical experiments validate the algorithm using ξ = ξ new in the case of a onedimensional Brownian Motion. The empirical variance σ2 is much smaller than p(1-p) nrep which would be obtained using a naive Monte-Carlo strategy (using n rep independent replicas). It is observed that the ratio between the empirical and the optimal variances increase when p decreases, but in practice this increase only has a limited impact and the new algorithm remains effective.

M = 10

Ornstein-Uhlenbeck

We consider here an example taken from [START_REF] Wouters | Rare event computation in deterministic chaotic systems using genealogical particle analysis[END_REF] . Let d = 1, and consider the diffusion process given by

dX(t) = -X(t)dt + dW (t), X(0) = 0.
The dynamics is discretized using the explicit Euler-Maruyama method, with time-step size ∆t = 10 -3 :

X n+1 = X n -∆tX n + √ ∆tζ n ,
with X 0 = X(0) = 0, where ζ n 0≤n≤N are independent standard Gaussian random variables.

The goal is to estimate the probability

p = P X N > a ,
for different values of a. This corresponds with the choice Φ(x) = x. The value of T is set to T = 2 so that N = T /∆t = 2000.

In this numerical experiment, the number of replicas is set equal to n rep = 10 2 , and the empirical average is computed over M = 10 4 independent realizations of the algorithm. The estimator pnew of p and the empirical variance σ2 new are obtained using the score function ξ = ξ new . The estimator pstd and the empirical variance σ2 std are obtained using the vanilla splitting strategy, with reaction coordinate ξ = ξ std . The value of the probability p for the continuous time process, and the optimal variance -p 2 log(p) nrep are also reported for comparison. is also reported, when the vanilla score function is used. This is the proportion of the independent realizations of the algorithm which contribute in the empirical average. This proportion depends on the conditional probability q (see ( 16)): it may happen that the n rep replicas obtained at the final iteration all satisfy X N ≤ 1, even if by construction they all satisfy max 0≤n≤N X n > 1. However, by construction (except if extinction of the system of replicas happens, which has not been observed in this experiment), if the new score function is used, r is identically equal to 1. Observe that, when T goes to infinity, by ergodicity of the process, p max → 1 (see (17) for the definition of p max ), whereas q and p converge to a non trivial probability. Thus, when T goes to infinity, it is expected that r will be equal to 0 if M is too small, when using the vanilla strategy with ξ = ξ std . The conditional probability q is also estimated: q = 0.07 when T = 2, q = 0.02 when T = 4, q = 0.01 when T = 8. To conclude this section, note that on this example, the AMS algorithms applied with the vanilla and the new score functions have a similar quantitative behavior in terms of asymptotic variance. However, their qualitative properties are different. When the conditional probability q gets small, the advantage of the new score function is the fact that the output p is always positive, so that even with a few realizations, one gets a rough but informative approximation of the target probability.

B. Drifted Brownian Motion

We here considers numerical examples taken from [START_REF] Bréhier | Unbiasedness of some generalized adaptive multilevel splitting algorithms[END_REF][START_REF] Rolland | Statistical behaviour of adaptive multilevel splitting algorithms in simple models[END_REF] . Let d = 1, and consider the diffusion process given by

dX t = -αdt + 2β -1 dW (t), X(0) = 0.
The dynamics is discretized using the explicit Euler-Maruyama method, with time-step size ∆t = 10 -2 (which gives again the exact solution in this simple case):

X n+1 = X n -α∆t + 2β -1 ∆tζ n ,
with X 0 = X(0) = 0, where ζ n 0≤n≤N are independent standard Gaussian random variables.

The goal is to estimate the probability

p = P X N > 1 ,
thus Φ(x) = x, a = 1. One considers the final time T = 1, so that N = T /∆t = 10 2 . The value of α is set equal to α = 4. As in the previous example, the value of p is easy to get using the fact that X N is Gaussian.

In this numerical experiment, see Table VI, three choices of score functions are considered.

The number of replicas is n rep = 10 3 . First, the estimator pnew and the empirical variance σ2 new are obtained using ξ = ξ new , with a sample size M = 4.10 4 . Second, the estimator pnew,a and the empirical variance σ2 new,a are obtained using ξ = ξ a new with a(t) = at T , with a sample size M = 4.10 5 . Finally, the estimator pstd and the empirical variance σ2 std are obtained using ξ = ξ std , with a sample size M = 4.10 4 . The sample sizes are chosen such that the total computational cost is of the same order for the three methods. Since the sample size is not the same for the three examples of score functions in Table VI, the values of the empirical variances σ2 should be taken with care when comparing the methods. One would rather compare the values of σ2 M . Then one concludes that the best performance is obtained when using ξ = ξ a new . The vanilla strategy, with ξ = ξ std , seems to behave quantitatively the same as when ξ = ξ new . However, the values of the proportion r of realizations such that pm = 0 is not zero is also reported, when ξ = ξ std (by construction, r = 1 for the first two cases). This means that if M was decreased (for instance, M of the order 10 2 for β = 4), then the output of the experiment would be pstd = 0.

As a consequence, the new algorithm clearly overcomes the limitation of the vanilla strategy when ξ = ξ std . However, the score functions are far from being optimal, as revealed by the comparison with the optimal variance.

C. Temporal averages for an Ornstein-Uhlenbeck process

In this section, we consider an example taken from [START_REF] Ferré | Adaptive sampling of large deviations[END_REF] . Consider the one-dimensional Ornstein-Uhlenbeck process X, which is the solution of the SDE

dX(t) = -X(t)dt + 2β -1 dW (t), X(0) = 0,
and define the temporal average

Y (T ) = 1 T T 0 X(s)ds.
More generally, set Y (t) = 1 t t 0 X(s)ds, for t ∈ (0, T ] and Y (0) = 0. The discretization is performed using the explicit Euler-Maruyama method, with timestep size ∆t = 5. 10 -3 : for n ∈ {0, . . . , N } with N ∆t = T ,

X n+1 = (1 -∆t)X n + 2β -1 ∆tζ n , Y n = 1 n n m=1 X m , with X 0 = Y 0 = 0. Note that Y satisfies a recursion formula Y n+1 = (1 -1 n+1 )Y n + 1 n+1 X n+1 .
The number of replicas is set equal to n rep = 10 3 and the sample size to compute empirical averages is M = 10 2 .

In this section, the probability which is estimated is

p(T, a) = P(Y N > a).
The associated estimator is denoted by p(T, a) and the empirical variance by σ2 (T, a).

In the large time limit T → ∞, since the law of Y (T ) converges to a centered Gaussian with variance 1, Y (T ) satisfies a large deviation principle, with rate function I defined by: for all a > 0, lim

T →∞ - 1 T log P Y (T ) > a = I(a) = a 2 4 .
In the numerical experiment, we illustrate the potential of the AMS algorithm to estimate the large deviations rate function. Notice that in the large T limit, the probability is extremely small and in practice cannot be estimated by the vanilla splitting strategy. The estimate of the rate function Î(a) is obtained by a regression procedure, see Figure 4. In addition to statistical error, two sources of numerical error are identified: values of T may not be sufficiently large, and the discretization of the dynamics and of the computation of temporal averages introduces a bias. The results, reported in Table VII, show the interest of this approach to estimate large deviations rate functionals.

D. Lorenz model

We consider the following stochastic version of the 3-dimensional Lorenz system, see [START_REF] Beck | Rare-event simulation[END_REF][START_REF] Lorenz | Deterministic nonperiodic flow[END_REF] for similar numerical experiments: Consider first the deterministic case, i.e. β = ∞. Then the system admits three unstable equilibria, and one of them is

           dX β 1 (t) = σ X 2 (t) -X 1 (t) dt + 2β -1 dW (t), dX β 2 (t) = rX 1 (t) -X 2 (t) -X 1 (t)X 3 (t) dt, dX β 3 (t) = X 1 (t)X 2 (t) -bX 3 (t),
x = b(r -1), b(r -1), r -1 = (5, 5 , 25).
Let the initial condition be given by X ∞ (0) = x + 1 2 (1, 1, 1). Then, one has the following stability result 5 : for all t ≥ 0, Φ X ∞ (t) ≤ 1, where

Φ(x) = x 2 1 (r + σ) 2 b σ + x 2 2 (r + σ) 2 b + x 3 -(r + σ) 2 (r + σ) 2 .
When noise is introduced in the system, i.e. β < ∞, we are interested in the estimation of the probability

P Φ(X β (T )) > 1 ,
with threshold a = 1.

In the numerical experiments, 2β -1 = 3, and the discretization is performed using the explicit Euler-Maruyama method, with time-step size ∆t = 10 -2 . The sample size is M = 10 4 , and the number of replicas is n rep = 10 3 . This numerical experiment thus illustrates the potential of the adaptive multilevel splitting algorithms introduced in this article, for applications to complex, nonlinear, stochastic models.

E. Driven periodic diffusion

We finally consider an example taken from [START_REF] Ferré | Adaptive sampling of large deviations[END_REF] . In this section, we consider the SDE on the unit circle, i.e. on the torus T,

dX(t) = -V (X(t)) + γ dt + √ 2dW (t),
where the potential energy function V (x) = cos(2πx) is periodic, and γ ∈ R. If γ = 0, this is called a non-equilibrium process since the drift term -V (x) + γ is not the derivative of a function defined on the torus T. In the remainder of this section, let γ = 1, and let the initial condition in the simulation be X 0 = 0. The discretization is performed using the Euler-Maruyama method, with time-step size ∆t = 10 -2 .

We are interested in the behavior of X T T , when T → ∞, more precisely we apply the AMS algorithm to estimate p(T, a) = P X(T ) > aT = P X(T ) T > a .

Following the same approach as for the temporal averages of the Ornstein-Uhlenbeck process, a large deviations rate function where comp.time(X) comp.time(Y ) is the ratio of the total computational times for the experiments using X and Y respectively. The values of Eff(Y |X) in this numerical experiment are reported in Table IX. We observe that Eff(Y |X) > 1 which means that the algorithm is more efficient using the process Y than the process X. To have a comparison with the committor score function, since the value of p(T, a) is not known, an approximation of the optimal variance is computed using the estimator pY (T, a).

I(a) = lim

FIG. 1 .

 1 FIG. 1. Level lines of the score function (t, x) → ξ new (t, x), with Φ(x) = x, a = 1, T = 1.

FIG. 2 .

 2 FIG. 2. Level lines of the score function (t, x) → ξ a new (t, x), with a(t) = atT -1 , and Φ(x) = x, a = 1, T = 1.

8 FIG. 3 .

 83 FIG. 3. Level lines of the committor function (t, x) → ξ com (t, x), in the Brownian Motion case X(t) = B(t), with a = 1 and T = 1.

1 M - 1 M

 11 ) is provided.Let n rep ∈ N denote the number of replicas, and let M ∈ N, and pm 1≤m≤M be the output probabilities of M independent realizations of the AMS algorithm. We report the values of the empirical average p of the empirical variance σ2 = m=1 (p m -p)[START_REF] Aristoff | CEMRACS 2013-modelling and simulation of complex systems: stochastic and deterministic approaches[END_REF] . Confidence intervals are computed as pnumber of realizations M is sufficiently large to use the Gaussian, Central Limit Theorem, regime.Recall that E[p] = E[p 1 ] = p, whatever the choice of the score function ξ and of the number of replicas n rep , thanks to Theorem 1. The variance of the estimator and thus the efficiency strongly depends on ξ. In the experiments below, the empirical variance σ2 is compared with the optimal asymptotic variance -p 2 log(p) nrep for (adaptive) multilevel splitting algorithms, which is obtained in the regime n rep → ∞, when choosing the (unknown in general) committor function ξ = ξ com as the score function. The difference between the empirical variance and the optimal one can be seen as a measure of how far the chosen score function is from the committor.

  10 -5 3.252 10 -5 3.213 10 -5 2.377 10 -10 3.327 10 -10 1.068 10 -10 2.9 1.756 10 -5 1.728 10 -5 1.742 10 -5 7.917 10 -10 1.009 10 -10 3.325 10 -11 3.0 9.341 10 -6 9.300 10 -6 9.260 10 -6 2.411 10 -11 5.832 10 -11 9.93710 -12

  10 -5 3.789 10 -5 3.740 10 -5 7.106 10 -10 1.089 10 -9 1.426 10 -10 0.85 2.9 2.052 10 -5 2.069 10 -5 2.049 10 -5 2.438 10 -10 3.484 10 -10 4.532 10 -11 0.83 3.0 1.104 10 -5 1.124 10 -5 1.101 10 -5 7.540 10 -11 1.147 10 -10 1.384 10 -11 0.81 3.1 5.822 10 -6 5.856 10 -6 5.805 10 -6 2.412 10 -11 3.584 10 -11 4.062 10 -12 0.78 3.2 3.022 10 -6 3.016 10 -6 3.002 10 -6 7.452 10 -12 9.964 10 -12 1.146 10 -12 0.75 TABLE IV. Ornstein-Uhlenbeck, T = 4. Comparison of the new and of the vanilla splitting algorithms. n rep = 10 2 , M = 10 4 .

  10 -5 3.714 10 -5 3.751 10 -5 1.696 10 -9 2.435 10 -9 1.434 10 -10 0.54 2.9 2.036 10 -5 2.071 10 -5 2.055 10 -5 5.439 10 -10 8.289 10 -10 4.557 10 -11 0.51 3.0 1.103 10 -5 1.128 10 -5 1.104 10 -5 1.843 10 -10 2.745 10 -10 1.392 10 -11 0.49 3.1 5.915 10 -6 5.968 10 -6 5.824 10 -6 5.599 10 -11 8.457 10 -11 4.089 10 -12 0.47 3.2 2.978 10 -6 3.022 10 -6 3.013 10 -6 1.639 10 -11 2.368 10 -11 1.154 10 -12 0.43 TABLE V. Ornstein-Uhlenbeck, T = 8. Comparison of the new and of the vanilla splitting algorithms. n rep = 10 2 , M = 10 4 .

a 2 4

 2 p(T = 25, a) p(T = 50, a) p(T = 100, a) p(T = 200, a) Î(a) a σ2 (T = 25, a) σ2 (T = 50, a) σ2 (T = 100, a) σ2 (T = 200, a) 0.4 7.28 10

3 and M = 10 2

 2 which depends on parameters σ, r, b and β. The parameters are given the following values in this section: σ = 3, r = 26 and b = 1.

2 FIG. 4 . 12 10 2 . 12 15 2 . 12 20 2 .

 24122122122 FIG. 4. Evolution of log p(T, a) as a function of T , for different values of a ∈ {0.4, 0.6, 0.8, 1.0, 1.2}, see Table VII.

T →∞ - 1 T

 1 log p(T, a) , is estimated, based on estimators of the probability p(T, a) for several values of T . In this numerical experiment, we compare two ways of applying the AMS algorithm, with the new score function ξ new but with different processes: considering either the process X(t) 0≤t≤T with the threshold aT , or the process Y (t) = X(t) t 0<t≤T , with the threshold a. Numerical values for different choices of a, T and n rep , of the associated estimators pX (T, a) and pY (T, a), and of the empirical variances σ2,X (T, a) and σ2,Y (T, a) are reported in Table IX below. It is observed that σ2,Y (T, a) < σ2,X (T, a), but a fair comparison requires to take into account the (average) computational cost. Thus the relative efficiency Eff(Y |X) of using the process Y instead of X, is computed as the ratio Eff(Y |X) = σ2,X comp.time(X) σ2,Y comp.time(Y ) ,

  Estimators Î(a) of the large deviations rate function I(a) are estimated by a regression procedure (with respect to T ) using the estimators pY (T, a), for several values of a. The numerical values are in excellent agreement with the numerical experiments in[START_REF] Ferré | Adaptive sampling of large deviations[END_REF] . The AMS algorithm introduced in this article can thus be an efficient tool to estimate large deviations rate functions. a T n rep pX (T, a) pY (T, a) σ2,X (T, a) σ2,Y (T, a) -(p Y (T,a)) 2 log(p Y (T,a)) nrep Eff(Y |X) Î(a) 0.8 100 10 2 8.483 10 -2

  -1 3.197 10 -1 [3.192 10 -1 , 3.205 10 -1 ] 1.256 10 -3 1.166 10 -3 4 1.613 10 -1 1.614 10 -1 [1.608 10 -1 , 1.617 10 -1 ] 5.275 10 -4 4.751 10 -4 8 4.978 10 -2 4.983 10 -2 [4.958 10 -2 , 4.998 10 -2 ] 1.030 10 -4 7.447 10 -5 16 6.395 10 -3 6.386 10 -3 [6.353 10 -3 , 6.436 10 -3 ] 4.449 10 -6 2.061 10 -6 32 1.631 10 -4 1.645 10 -4 [1.611 10 -4 , 1.651 10 -4 ] 1.052 10 -8 2.358 10 -9 64 1.787 10 -7 1.782 10 -7 [1.721 10 -7 , 1.854 10 -7 ] 1.139 10 -13 4.935 10 -15 128 2.501 10 -13 3.011 10 -13 [1.766 10 -13 , 3.235 10 -13 ] 1.403 10 -23 2.614 10 -26

	β	p	p (∆t = 0)	confidence interval	σ2	-p 2 log(p) nrep
	2 3.199 10

[START_REF] Au | Estimation of small failure probabilities in high dimensions by subset simulation[END_REF]

independent realizations of the algorithm. Second, in Table

II

, the number of replicas is set equal to n rep = 10 3 , and the empirical average is computed over M = 10 3 independent realizations of the algorithm.

TABLE I .

 I Brownian Motion, n rep = 10 2 , M = 10 4 .

	β	p	p (∆t = 0)	confidence interval	σ2

-p 2 log(p) nrep 2 3.196 10 -1 3.197 10 -1 [3.188 10 -1 , 3.203 10 -1 ] 1.362 10 -4 1.166 10 -4 4 1.617 10 -1 1.614 10 -1 [1.612 10 -1 , 1.621 10 -1 ] 5.456 10 -5 4.751 10 -5 8 4.983 10 -2 4.983 10 -2 [4.963 10 -2 , 5.002 10 -2 ] 9.815 10 -6 7.447 10 -6 16 6.411 10 -3 6.386 10 -3 [6.371 10 -3 , 6.451 10 -3 ] 4.242 10 -7 2.061 10 -7 32 1.634 10 -4 1.645 10 -4 [1.614 10 -4 , 1.655 10 -4 ] 1.063 10 -9 2.358 10 -10 64 1.800 10 -7 1.782 10 -7 [1.740 10 -7 , 1.860 10 -7 ] 9.360 10 -15 4.935 10 -16 128 3.045 10 -13 3.011 10 -13 [2.426 10 -13 , 3.664 10 -13 ] 9.986 10 -25 2.614 10 -27

TABLE II .

 II Brownian Motion, n rep = 10 3 , M = 10 3 .

  3.1 4.857 10 -6 4.826 10 -6 4.827 10 -6 7.011 10 -12 8.482 10 -12 2.852 10 -12 3.2 2.486 10 -6 2.449 10 -6 2.468 10 -6 1.984 10 -12 2.475 10 -12 7.864 10 -13

TABLE III

 III The results of numerical experiments with T = 4 and T = 8 are reported below. The

	quantity
	r =	1 M	M m=1	1 pm>0

. Ornstein-Uhlenbeck, T = 2. Comparison of the new and of the vanilla splitting algorithms. n rep = 10 2 , M = 10 4 .

  -4 2.036 10 -4 2.033 10 -4 2.035 10 -4 1.572 10 -9 5.734 10 -9 3.525 10 -9 3.519 10 -10 0.99 2 2.843 10 -7 2.870 10 -7 2.878 10 -7 2.867 10 -7 4.714 10 -14 2.348 10 -12 9.527 10 -14 1.238 10 -15 0.69 3 4.613 10 -10 4.325 10 -10 4.705 10 -10 4.571 10 -10 1.817 10 -18 2.197 10 -17 3.084 10 -18 4.493 10 -21 0.11 4 7.620 10 -13 6.975 10 -13 7.582 10 -13 7.687 10 -13 5.034 10 -23 4.388 10 -22 8.343 10 -23 1.648 10 -26 0.01 TABLE VI. Drifted Brownian Motion, β ∈ {1, 2, 3, 4}, Comparison of two versions of the new splitting algorithm and of the vanilla splitting algorithm.

	β	pnew	pnew,a	pstd	p	σ2 new	σ2 new,a	σ2 std	-p 2 log(p) nrep	r
	1 2.037 10

TABLE VII .

 VII Temporal averages for an Ornstein-Uhlenbeck process. n rep = 10

		-2	2.12 10 -2	2.22 10 -3	2.75 10 -5	0.045 0.040
	-	2.02 10 -5	3.25 10 -6	1.27 10 -7	3.85 10 -11	-	-
	0.6 1.45 10 -2	1.16 10 -3	1.16 10 -5	7.28 10 -10 0.096 0.090
	-	1.36 10 -6	5.85 10 -8	2.21 10 -10	6.44 10 -19	-	-
	0.8 1.76 10 -3	2.57 10 -5	6.12 10 -9	2.73 10 -16 0.169 0.160
	-	8.01 10 -8	3.32 10 -10	1.44 10 -16	1.98 10 -31	-	-
	1.0 1.37 10 -4	1.67 10 -7	3.06 10 -13	1.71 10 -24 0.261 0.250
	-	1.47 10 -9	2.07 10 -14	1.16 10 -25	9.32 10 -48	-	-
	1.2 6.21 10 -6	4.83 10 -10	2.71 10 -18	2.88 10 -34 0.373 0.360
	-1.40 10 -11	6.69 10 -19	3.49 10 -35	5.21 10 -67	-	-

TABLE VIII .

 VIII Lorenz model. n rep = 10 3 and M = 10 4

  8.489 10 -2 7.136 10 -4 2.487 10 -4 1.777 10 -4 1.0 -200 -2.647 10 -2 2.776 10 -2 1.832 10 -4 3.519 10 -5 2.762 10 -5 1.7 0.0112 1 50 10 3 1.529 10 -2 1.505 10 -2 7.046 10 -6 1.513 10 -6 9.504 10 -7 3.0 -100 -1.026 10 -3 1.085 10 -3 2.586 10 -7 3.879 10 -8 8.036 10 -9 2.8 0.0526 1.25 50 10 3 1.374 10 -4 1.311 10 -4 1.227 10 -8 1.355 10 -9 1.537 10 -10 4.9 -100 -8.941 10 -7 1.017 10 -7 1.048 10 -13 1.585 10 -15 1.665 10 -16 35 0.189 TABLE IX. Estimates of P(X T > aT ) and of I(a) for the Periodic driven diffusion. The sample size is M = 100.

ACKNOWLEDGMENTS

The authors would like to thank G. Ferré and H. Touchette for discussions concerning the numerical experiments, and F. Cérou, A. Guyader and M. Rousset for stimulating discussions at early stages of this work. The work of T. Lelièvre is supported by the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement number 614492.