Pierre Étoré
email: pierre.etore@imag.fr

Miguel Martinez
email: miguel.martinez@univ-mlv.fr

Exact simulation of one-dimensional stochastic differential equations involving the local time at zero of the unknown process

Keywords: Exact simulation methods, Skew Brownian motion, One-dimensional diffusion, Local Time

In this article we extend the exact simulation methods of Beskos et al. in [4] to the solutions of one-dimensional stochastic differential equations involving the local time of the unknown process at point zero. In order to perform the method we compute the law of the skew Brownian motion with drift. The method presented in this article covers the case where the solution of the SDE with local time corresponds to a divergence form operator with a discontinuous coefficient at zero. Numerical examples are shown to illustrate the method and the performances are compared with more traditional discretization schemes.

Introduction

Presentation

Exact simulation methods for trajectories of one-dimensional SDEs has been a subject of much interest in the last years : see for example [START_REF] Beskos | Exact simulation of diffusions[END_REF], [START_REF] Beskos | Retrospective exact simulation of diffusion sample paths with applications[END_REF], [START_REF] Beskos | A factorisation of diffusion measure and finite sample path constructions[END_REF], [START_REF] Reutenauer | Exact simulation of prices and greeks: application to cir[END_REF], [START_REF] Sbai | Modélisation de la dépendance et simulation de processus en finance[END_REF]. Unlike the classical simulation methods which all involve some kind of discretization error (we mention [START_REF] Bally | The law of the Euler scheme for stochastic differential equations: I. Convergence rate of the distribution function[END_REF] for the Euler Scheme), the exact simulation methods are constructed in such a way that they do not present any discretization error under the strong hypothesis that the diffusion coefficient is constant and equal to one. In the last years, the original method presented in the fundamental article [START_REF] Beskos | Retrospective exact simulation of diffusion sample paths with applications[END_REF] has been extended to overcome various limitations of the initial algorithm ; it has been generalized to include the cases of unbounded drifts ( [START_REF] Beskos | A factorisation of diffusion measure and finite sample path constructions[END_REF], [START_REF] Beskos | MCMC methods for diffusion bridges[END_REF]).

On another hand, the numerical simulation of SDEs corresponding to divergence form operators

L = ∇. 1 2 a(x)∇
involving a discontinuous coefficient a has been also the subject of various studies in the last years. Indeed these operators are of great importance since they appear in a wide range of modelling problems involving diffusion phenomena in discontinuous media. Among applications, we can mention ecology ( [START_REF] Cantrell | Diffusion models for population dynamics incorporating individual behavior at boundaries: Applications to refuge design[END_REF]), geophysics ( [START_REF] Lejay | Simulating a diffusion on a graph. Application to reservoir engineering[END_REF], [START_REF] Lejay | Monte Carlo methods for fissured porous media: a gridless approach[END_REF]), astrophysics ( [START_REF] Zhang | Calculation of diffusive shock acceleration of charged particles by skew brownian motion[END_REF]), or magneto/electroencephalography ( [START_REF] Faugeras | The Inverse EEG and MEG Problems : The Adjoint State Approach I: The Continuous Case[END_REF]).

In the one-dimensional context, various Random Walks and an Euler Scheme have been studied for the simulation of the solution of such SDEs : for Random Walks we mention [START_REF] Etoré | On random walk simulation of one-dimensional diffusion processes with discontinuous coefficients[END_REF], [START_REF] Etoré | Approximation de processus de diffusion à coefficients discontinus en dimension un et applications à la simulation[END_REF], [START_REF] Etoré | A Donsker theorem to simulate one-dimensional processes with measurable coefficients[END_REF], [START_REF] Lejay | A scheme for simulating one-dimensional diffusion processes with discontinuous coefficients[END_REF] ; for the Euler Scheme see [START_REF] Martinez | Interprétations probabilistes d'opérateurs sous forme divergence et analyse des méthodes numériques probabilistes associées[END_REF], [START_REF] Martinez | Discrétisation d'équations différentielles stochastiques unidimensionnelles à générateur sous forme divergence avec coefficient discontinu[END_REF], [START_REF] Martinez | One-dimensional parabolic diffraction equations: pointwise estimates and discretization of related stochastic differential equations with weighted local times[END_REF] in the case where the discontinuity of the coefficient in the divergence operator appears at point 0. Of course, for such SDEs, the order of discretization error of these discretization schemes is usually greater than those obtained in a more classical context.

An important problem comes from the fact that SDEs corresponding to divergence form operators do not enter the classical scope of SDEs covered by the exact simulation methods.

The main difficulty is that these SDEs include an additional term, which involves in dimension one the local time of the unknown process (in dimension greater than one, it involves the local time of a one-dimensional auxiliary process; see [START_REF] Bossy | Probabilistic interpretation and random walk on spheres algorithms for the Poisson-Boltzmann equation in molecular dynamics[END_REF]). In fact, the laws of the solution of such one-dimensional SDEs are no longer absolutely continuous with respect to the Wiener measure.

In this paper we present a first attempt for the adaptation of the exact simulation methods of [START_REF] Beskos | Retrospective exact simulation of diffusion sample paths with applications[END_REF] to one-dimensional SDEs with an additional term that involves the local time of the unknown process at point 0. Namely, our object of study is (X t ) t≥0 solution of dX t = σ(X t )dW t + b(X t )dt + βdL 0 t (X),

X 0 = x, (1)
where |β| < 1 (with β = 0) and L 0 t (X) is the symmetric local time of X in zero at time t ; the reason why we only deal with |β| < 1 is that there is no solution to (1) when β > 1 (the case where |β| = 1 corresponds to a reflected diffusion and we do not include it in our discussion). The reason why we restrict ourselves to the case β = 0 is made clear in Section 3 (see Remark 2 and also the conclusion of the paper).

Note that when σ is identically equal to 1 and b is identically equal to 0, the solution (X t ) of ( 1) is a standard Skew Brownian Motion (SBM in short). Under mild assumptions concerning b and standard ellipticity conditions on σ, it is known that there exists a unique strong solution (X t ) t≥0 to (1) (see [START_REF] Gall | One-dimensional stochastic differential equations involving the local times of the unknown process[END_REF] for details).

Let us emphasize that this work includes the situation where b may be discontinuous at 0. So that the results of this paper are also suited for the situation stated in [START_REF] Martinez | Discrétisation d'équations différentielles stochastiques unidimensionnelles à générateur sous forme divergence avec coefficient discontinu[END_REF], [START_REF] Martinez | One-dimensional parabolic diffraction equations: pointwise estimates and discretization of related stochastic differential equations with weighted local times[END_REF] where the solution of (1) corresponds to a divergence form operator whose coefficient is discontinuous at 0 (and is sufficiently smooth elsewhere). We show a numerical example to illustrate this interesting case.

Let us now briefly explain our main idea. When σ ≡ 1, we show that the law of (X t ) t≥0 (solution of [START_REF] Appuhamillage | Occupation and local times for skew brownian motion with applications to dispersion across an interface[END_REF]) is absolutely continuous with respect to the law of some Skew Brownian Motion (SBM) with a drift component. The reason why the SBM with drift appears naturally in our computations is explained in Section 3 (see Remark 1).

So, contrary to the already mentioned discretization schemes where the standard SBM is used in force, we do not longer deal with a simple SBM but with a SBM that possesses a drift component. As a consequence, in order to adapt the method of [START_REF] Beskos | Retrospective exact simulation of diffusion sample paths with applications[END_REF] in this setting, we have to be able to simulate bridges of the SBM with drift.

In the last section of the paper, we discuss the limitations of the initial algorithm. The main issue is to relax the boundedness assumptions made on the drift function b, as is done in [START_REF] Beskos | A factorisation of diffusion measure and finite sample path constructions[END_REF] for "classical" SDEs. In [START_REF] Beskos | A factorisation of diffusion measure and finite sample path constructions[END_REF], the authors use some kind of factorisations for the sample state space of the standard Brownian Bridge, which are consequences of William's decomposition theorem for Brownian Motion. Proving similar factorisations for the Skew Brownian Bridge with drift seems difficult to us. Nevertheless, we have been able to apply a result stated in Pitman-Yor [START_REF] Pitman | Decomposition at the maximum for excursions and bridges of onedimensional diffusions[END_REF] in the case of the standard Skew Brownian Bridge, which gives a first partial result. Unfortunately, we have not been able to relax the boundedness assumption on the drift function b and we think that much remains to do in this direction.

In our concluding remarks we discuss the particular problem of being able to produce an exact simulation algorithm in the case β = 0, which in our opinion should be regarded as a special separate problem. We also draw bold lines for further investigation.

Organisation of the paper

In Section 2 we precise the hypotheses and define the problem we will deal with. We also introduce notations used in the sequel. In Section 3 we present the exact simulation algorithm, adapted from Beskos et al. to our situation. It turns out that, in order to use the algorithm, we need to sample bridges of a Skew Brownian Motion with drift. Section 4 is devoted to the computation of the transition probability density of the Skew Brownian Motion with drift. Then Section 5 explains how to sample bridges of a Skew Brownian Motion with drift, using rejection sampling with brownian bridges as proposals. Section 6 presents numerical experiments, including a divergence form case. Finally, we discuss possible extensions in Section 7.

2. Exposition of the exact simulation problem.

Exact simulation problem and first assumptions

Denote C = C([0, T ], R) the set of continuous mappings from [0, T ] to R and C the Borel σ-field on C induced by the supreme norm.

Let P be a probability measure on (C, C) and W a Brownian motion under P together with its completed natural filtration (F t ) t≥0 . We will denote

P x = P (• | W 0 = x).
Throughout the whole paper, we will make the following assumptions

-|β| < 1.
-The function b : R → R is bounded and differentiable on R * ,+ and R * ,-with a possible discontinuity at point {0}. We suppose that both limits lim We seek for an exact simulation algorithm of the paths of the solution of the one-dimensional Stochastic Differential Equation

dX t = dW t + b(X t )dt + βdL 0 t (X), X 0 = x, (2)
where L 0 t (X) is the symmetric local time of X in zero at time t.

Some recalls on

dY t = 1/(g -1 ± ) ′ (Y s )dW s + b • g -1 /(g -1 ±
) ′ (Y s )ds without local time. Here (g -1 ± ) ′ denotes the half sum of the right and left derivatives of g -1 . Since (g -1 ± ) ′ is bounded from below by a strictly positive constant, this transformed equation makes sense and well-known results for classical one dimensional SDEs (see for example [START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF] Chap 5. Section 5.5) ensure that under the assumptions of Subsection 2.1, it possesses a weak solution. Then, g -1 (Y ) gives a weak solution of (2).

The difficulty concerns strong uniqueness. Strong uniqueness for the solutions of equation ( 2) is proved with a direct application of Theorem 1.3 p.55 in the fundamental article [START_REF] Gall | One-dimensional stochastic differential equations involving the local times of the unknown process[END_REF], which deals with a broader class of stochastic differential equations involving the local time of the unknown process.

Note that when β = -1 or β = +1, equation (2) possesses a unique strong solution, which is a reflected diffusion at 0, either reflected below 0 (β = -1) or above 0 (β = +1).

Let us now briefly explain why there is no solution to equation (2) when |β| > 1. Remember that we are working with the symmetric local time L 0 t (X). Let us denote by L 0,r t (X) (resp. L 0,l t (X)) the right-hand sided local time of the process X (resp. the left hand sided local time of X). It is an exercise to prove that if X is a solution of (2), then L 0 t (X) = 1+β 2 L 0,r t (X) and L 0 t (X) = 1-β 2 L 0,l t (X). In particular, we see that when |β| > 1 there is no solution to (2) (otherwise the symmetric local time of X would be negative !).

Strong Markov property

The proof of the strong Markov property for solutions of equation ( 2) is a separate problem from the one of existence and uniqueness. We refer to [START_REF] Kulik | On the solution of a one-dimensional stochastic differential equation with a singular drift coefficient[END_REF] for a rigorous proof. In the multidimensional context of diffusion processes with generalized drift, the proof of the Markov property may be found in [START_REF] Ludmila | On the Markov property of strong solutions to SDE with generalized coefficients[END_REF]. Throughout the whole paper, we use the following notations :

Exact simulation

-We note

N c (x) := 1 √ 2π ∞ x e -y 2 2 dy.
-We set

µ := 1 + β 2β b(0+) - 1 -β 2β b(0-). (3)
and define b(z) := b(z)µ.

We set

z → φ(z) := b 2 (z) + b ′ (z) + 2µb(z) 2 ½ R * ,+ ∪R * ,- (z).
-We set φ(z) := φ(z)m with m = inf z∈R φ(z); the constant K denotes an upper bound of the function φ.

-We set B(u

) := u 0 b(y)dy, u ∈ R.
-B β,µ will denote the SBM of parameter β and drift µ. That is to say B β,µ is the strong solution of (2) in the case b ≡ µ, namely :

dB β,µ t = dW t + µdt + βdL 0 t (B β,µ ). (4)
We will denote p β,µ (t, x, y) the transition probability density of B β,µ .

Note that, with this notation, p 0,µ (t, x, y) is the transition probability density of the Brownian motion with constant drift µ ∈ R, namely

p 0,µ (t, x, y) = 1 √ 2πt exp - (y -x -µt) 2 2t . (5)
Note also that p β,0 (t, x, y) is the transition probability density of the SBM of parameter β (without drift), see [START_REF] Lejay | On the constructions of the skew Brownian motion[END_REF].

-We set τ 0 := inf t ≥ 0 : B β,µ t = 0 with the convention inf(∅) = +∞.

We will denote by h(x, .) the density of τ 0 under P x .

Additional assumptions

In this section, we will make the following additional assumptions :

β = 0 (for the reason of this last hypothesis, see Remark 2).

-We assume that the function z → φ(z) is bounded.

-We assume that the function u → exp[B(u) -(ux) 2 /2T ] is integrable.

3.2. Presentation of the exact simulation algorithm 3.2.1. Change of probability Recall that in this case b(z) := b(z)µ where µ is the constant defined by

µ := 1 + β 2β b(0+) - 1 -β 2β b(0-).
Note that since β = 0 by assumption, this constant is well-defined. In the case where b is continuous at point {0}, observe that µ reduces to b(0).

We have

dX t = dW t + b(X t )dt + µdt + βdL 0 t (X).
In particular, we may perform Girsanov's theorem (see Theorem 3.5.1 in [START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF]) and we write

dX t = dW SD t + µdt + βdL 0 t (X), (6)
where

W SD t := W t + t 0 b(X s )
ds is a Brownian motion under the new probability W SD defined by

dP dW SD = exp T 0 b(X t )dW SD t - 1 2 T 0 b 2 (X t )dt . (7)
From our assumptions on b, we are in position to apply the symmetric Itô-Tanaka formula to the function B(u) := u 0 b(y)dy and (X t ) t≥0 .

Applying the occupation's time formula, we obtain

B(X T ) -B(x) = T 0 b(X t +) + b(X t -) 2 dX t + 1 2 T 0 b ′ (X t )½ Xt =0 dt + b(0+) -b(0-) 2 L 0 T (X) = T 0 b(X t )½ Xt =0 dW SD t + µ T 0 b(X t )½ Xt =0 dt + 1 2 T 0 b ′ (X t )½ Xt =0 dt + b(0+) + b(0-) 2 β + b(0+) -b(0-) 2 = 0 L 0 T (X) (8)
where the last line comes from the definitions of b and µ (and the property dL 0 t (X) = 1 Xt=0 dL 0 t (X)). From the fact that ℓ{t ∈ [0, T ] : X t = 0} = 0 (where ℓ stands for the Lebesgue measure), we see that

B(X T ) -B(x) = T 0 b(X t )dW SD t + µ T 0 b(X t )dt + 1 2 T 0 b ′ (X t )dt. (9)
Thus, [START_REF] Beskos | Exact simulation of diffusions[END_REF] implies that for any functional F (X) of the path up to time T , one has :

E P [F (X)] = E W SD F (X) exp B(X T ) -B(x) - T 0 φ(X t )dt , where φ(z) = b 2 (z) + b ′ (z) + 2µb(z) 2 .
Remark 1 Note that, because of the definition of b, there is no local time appearing in equality [START_REF] Cantrell | Diffusion models for population dynamics incorporating individual behavior at boundaries: Applications to refuge design[END_REF] after the application of the Itô-Tanaka formula. This ensures that there is no local time involved in the exponential martingale of Girsanov's theorem, which makes it tractable for a numerical perspective.

Retrospectively, this explains why in the sequel we have to deal with a Skew Brownian Motion with drift instead of a simple standard SBM.

Remark 2 We now explain our assumption β = 0.

Note that in the case β = 0, the constant µ is no more defined. In fact, in the case β = 0, because of the discontinuity of b, it is no longer possible to get rid of the local time as in [START_REF] Bossy | Probabilistic interpretation and random walk on spheres algorithms for the Poisson-Boltzmann equation in molecular dynamics[END_REF]. More precisely, there is no constant θ such that proceeding as the computations in ( 6) and ( 7) with b(x) := b(x)θ we can cancel the local time term appearing in the exponential weight.

For a more detailed discussion on the case β = 0 the interested reader is invited to read the conclusion at the end of this paper.

Exact simulation algorithm (after Beskos and al)

Considering ( 6), we see that the law of X under W SD is given by p β,µ (t, x, y)dy.

Following the lines of Beskos et al. in [START_REF] Beskos | Retrospective exact simulation of diffusion sample paths with applications[END_REF], and considering the computations performed in the above section, we give an algorithm that returns an exact simulation of a skeleton of (X t ) t∈[0,T ] solution of (2) starting from x 0 .

EXACT SIMULATION ALGORITHM FOR A SOLUTION OF (2) starting from x 0 .

1. Simulate a random variable Z according to the density

h(y) = C exp (B(y) -B(x 0 )) p β,µ (T, x 0 , y),
where C is the normalizing constant such that h(y)dy = 1. Keep in memory the value z of Z.

Simulate a Poisson Point Process with unit density on

[0, T ] × [0, K]. The result is a random number n of points of coordinates (t 1 , z 1 ), . . . , (t n , z n ). 3. Simulate a skeleton (B β,µ t 1 , . . . , B β,µ tn ) conditioned on B β,µ 0 = x 0 and B β,µ T = z. 4. If ∀i ∈ {1, . . . , n} φ(B β,µ t i ) ≤ z i accept the skeleton. Else return to step 1.
This algorithm returns an exact sample of (X t 1 , . . . , X tn , X T ) (in particular we get an exact simulation of X T , it is the value z of Z used for an accepted trajectory).

Note that in order to apply the methodology of [START_REF] Beskos | Retrospective exact simulation of diffusion sample paths with applications[END_REF] we have to be able to generate bridges of a drifted Skew Brownian Motion B β,µ . Indeed, this is the key one has to reach for in order to perform the Step 3.

Computation of the law of the SBM with drift

Recalls on the construction of the SBM using a "random flipping" of excursions

In this paragraph, we present a construction of the Skew Brownian Motion -solution of (4) with µ = 0 and starting from x = 0 -that gives an understanding of its relation with the standard Brownian motion. The construction is made out from a reflecting Brownian Motion with a change of sign of each excursion with probability (1β) /2. It is explained in [START_REF] Revuz | Continuous martingales and Brownian motion[END_REF] page 487 exercise 2.16 (we use the same notations as [START_REF] Revuz | Continuous martingales and Brownian motion[END_REF] in the explanations below).

Suppose that x = 0. The construction is as follows : let (Y n ) n≥0 be a sequence of independent r.v.'s taking the values 1 and -1 with probabilities (1 + β) /2 and (1β) /2 and independent of some Brownian Motion B. Let us set F B t t≥0 the natural filtration generated by B and satisfying the usual right continuous and completeness conditions and F B := t≥0 F B t . We also denote H := σ(Y n : n ≥ 0) the corresponding σ-algebra generated by the whole sequence (Y n ) n≥0 and E := σ(e s : s ≥ 0) the σ-algebra generated by all the excursions of B so that E ⊂ F B . The "good" time clock for the point-wise excursion process (e s ) s≥0 is the process of time-change (τ t ) t≥0 defined as the r.c.l inverse of the local time L 0 t (B) t≥0 : so that the excursion process (e s ) s≥0 may be viewed as a Point Poisson Process on C 0→0 running in the local time scale.

For each ω in the set on which B is defined, the set of excursions e s (ω) is countable and may be ordered. Define a process B β by putting

B β t (ω) = Y ns(e) (ω)|e s (t -τ s-(ω), ω)| if τ s-≤ t ≤ τ s
and where e s is the n s (e)-th excursion in the above ordering. The random number n s (e) is a random variable measurable w.r.t. E, which depends on the whole excursion process e = (e u ) u>0 and the time variable s in the local time scale. It may be proved that the process thus obtained is a Markov process and that it is a Skew Brownian motion of parameter β starting at x = 0.

By construction the sigma algebras H = σ(Y n : n ≥ 0) and F B are independent ; in particular, if we denote by R(e) the end point of excursion e, then we have τ t (ω) = s≤t R(e s (ω)) and thus τ t is measurable w.r.t. F B . Note that this construction implies that almost surely,

L 0 t (B β ) = L 0 t (|B|) = L 0 t (B)
for any t ≥ 0. So L t (B β ) t≥0 may be recovered as the r.c.l inverse of (τ t ) t≥0 . Consequently, it is adapted to F B t t≥0 and it is independent of H.

Remark 3 For a possible extension of this "random flipping of excursions" method for the construction of the solution to the more general equation ( 2) (and possibly solutions of ( 1)), we refer to the article of Lejay [START_REF] Lejay | On the decomposition of excursions measures of processes whose generators have diffusion coefficients discontinuous at one point. Markov Process[END_REF], which gives a decomposition of the Itô measure associated to X in the general context of solutions of (1). At least in the context of equation ( 2), the result stated in [START_REF] Lejay | On the decomposition of excursions measures of processes whose generators have diffusion coefficients discontinuous at one point. Markov Process[END_REF] should allow to perform a construction along the same lines as above, flipping the excursions of some reflected process whose law should be the same as the law of |X|. However and up to our knowledge, such construction has never been explicitly written down in the literature, even in the context of equation (2).

Computation of the joint law of SBM and its local time

Let us begin with a direct consequence of the construction explained above in Subsection 4.1.

Lemma 1 Let W be a Brownian motion defined on (C, C, P) and B β the strong solution of (4) with µ = 0.

We have for all t > 0,

P 0 |B β t | ∈ dy; L 0 t (B β ) ∈ dl = P 0 |W t | ∈ dy; L 0 t (W ) ∈ dl . ( 10
)
Remark 4 We even have that the process

(|B β t |, L 0 t (B β )) t≥0 is distributed as (|W t |, L 0 t (W )
) t≥0 under P 0 . This is a clear consequence of the construction explained in paragraph 4.1. Another way to prove this fact is to check that their common distribution is the one of (M W t -W t , M W t ) t≥0 under P 0 , where M W t = max 0≤s≤t W s . This is related to the Lévy theorem, as stated for instance in Theorem 3.6.17 in [START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF], where it is proved by using the Skorokhod method.

Let us now state an intuitive result, which is somewhat not so easy to prove without using the construction explained in Subsection 4.1. The difficulty comes from the presence of the local time in the equalities below.

Lemma 2

We have for all t > 0,

P 0 B β t ∈ dy; L 0 t (B β ) ∈ dl = 1 + β 2 P 0 |B β t | ∈ dy; L 0 t (B β ) ∈ dl + 1 -β 2 P 0 -|B β t | ∈ dy; L 0 t (B β ) ∈ dl .
Proof. We start from the construction of SBM using a random flipping of excursions coming from a reflected Brownian Motion as explained at paragraph 4.1.

Let S be the space of real sequences (a k ) k∈N and denote Φ : S × N → R the coordinate function defined by Φ((a k ) k∈N , n) = a n . From the independence of H and F B t and since L 0 t (B β ) t≥0 is adapted w.r.t. F B t t≥0 , from the properties of the conditional expectation,

1 y≥0 P 0 B β t ∈ dy; L 0 t (B β ) ∈ dl = 1 y≥0 E 0 P 0 B β t ∈ dy; L 0 t (B β ) ∈ dl | F B = E 0 P 0 Y ns(e) > 0; |e s (t -τ s-(ω), ω)| ∈ dy; L 0 t (B β ) ∈ dl | F B = E 0 P 0 Φ((Y k ) k∈N (ω), n s (e(ω))) > 0 | F B ; |e s (t -τ s-(ω), ω)| ∈ dy; L 0 t (B β ) ∈ dl = E 0 P 0 Y n > 0 | n=ns(e(ω)) ; |e s (t -τ s-(ω), ω)| ∈ dy; L 0 t (B β ) ∈ dl = 1 + β 2 P 0 |B β t | ∈ dy; L 0 t (B β ) ∈ dl .
Proceeding similarly on R * -we get,

1 y<0 P 0 B β t ∈ dy; L 0 t (B β ) ∈ dl = 1 -β 2 P 0 -|B β t | ∈ dy; L 0 t (B β ) ∈ dl ,
therefore the result.

Using the two last lemmas we can prove the following result.

Proposition 1 Let W be a Brownian motion defined on (C, C, P) and B β the strong solution of (4) with µ = 0.

We have for all t > 0, x ≥ 0,

P x B β t ∈ dy; L 0 t (B β ) ∈ dl = 1 y≥0 1 l>0 (1+β)(l+y+x) √ 2πt 3 exp -(l+y+x) 2 2t dydl +1 y≥0 1 √ 2πt exp{-(y-x) 2 2t } -exp{-(y+x) 2 2t } dyδ 0 (dl) +1 y<0 1 l≥0 (1-β)(l-y+x) √ 2πt 3 exp -(l-y+x) 2 2t
dydl.

Remark 5 The result of Proposition 1 appears as a corollary of a more general result stated in [START_REF] Ludmila | On a multidimentional brownian motion with partly reflecting membrane on a hyperplane[END_REF]. Note also that the result of Proposition 1 (and consequently the result stated in Proposition 2 in the next Section) differ slightly from results published by T. Appuhamillage et al. in the recent article [START_REF] Appuhamillage | Occupation and local times for skew brownian motion with applications to dispersion across an interface[END_REF] where there is a computational error (see also [START_REF] Appuhamillage | Occupation and local times for skew Brownian motion with applications to dispersion across an interface[END_REF] for a discussion). We will detail the computations for the sake of completeness and clarification.

Proof of Proposition 1.

Step 1. Combining the results of the Lemmas 1 and 2 we have

1 y≥0 P 0 B β t ∈ dy; L 0 t (B β ) ∈ dl = 1+β 2 P 0 |B β t | ∈ dy; L 0 t (B β ) ∈ dl = 1+β 2 P 0 |W t | ∈ dy; L 0 t (W ) ∈ dl .
Step 2. Let x > 0. As 1 y≥0 1 l>0 P x [B β t ∈ dy; L 0 t (B β ) ∈ dl; t < τ 0 ] = 0, we have, using the strong Markov property,

1 y≥0 1 l>0 P x [ B β t ∈ dy; L 0 t (B β ) ∈ dl ] = 1 y≥0 1 l>0 P x [ B β t ∈ dy; L 0 t (B β ) ∈ dl; t ≥ τ 0 ] = 1 y≥0 1 l>0 E x [1 {t≥τ 0 } P x [ B β t ∈ dy; L 0 t (B β ) ∈ dl |F τ 0 ]] = 1 y≥0 1 l>0 t 0 P 0 [B β t-s ∈ dy; L 0 t-s (B β ) ∈ dl ]h(x, s)ds,
where h(x, .) is the density of τ 0 under P x . But h(x, .) is also the density of T 0 = inf(t ≥ 0 : W t = 0). And using the first step of the proof we have

1 y≥0 1 l>0 P 0 [B β t-s ∈ dy; L 0 t-s (B β ) ∈ dl ] = 1 l>0 1 + β 2 P 0 [ |W t-s | ∈ dy; L 0 t-s (W ) ∈ dl ].
Using again the strong Markov property we get

1 y≥0 1 l>0 P x [ B β t ∈ dy; L 0 t (B β ) ∈ dl ] = 1 l>0 E x [1 {t≥T 0 } 1+β 2 P x [ |W t | ∈ dy; L 0 t (W ) ∈ dl |F T 0 ]] = 1 l>0 1 + β 2 P x [ |W t | ∈ dy; L 0 t (W ) ∈ dl; t ≥ T 0 ] = 1 l>0 1 + β 2 P x [ |W t | ∈ dy; L 0 t (W ) ∈ dl ] = 1 y≥0 1 l>0 (1 + β)(l + y + x) √ 2πt 3 exp - (l + y + x) 2 2t dydl.
Step 3. It is a consequence of the reflection principle that

1 y≥0 P x [ B β t ∈ dy; L 0 t (B β ) = 0 ] = 1 y≥0 P x [ B β t ∈ dy; t < τ 0 ] = 1 y≥0 1 √ 2πt exp{- (y -x) 2 2t } -exp{- (y + x) 2 2t } .
Using Step 1 to 3 we have the result for x ≥ 0 on R + × R + . In order to retrieve the result on R * -× R + we use Step 1 and 2 with 1+β 2 replaced by 1-β 2 and 1 y≥0 replaced by 1 y<0 , and the fact that for x ≥ 0,

1 y<0 P x [ B β t ∈ dy; L 0 t (B β ) = 0 ] = 1 y<0 P x [ B β t ∈ dy; t < τ 0 ] = 0.

The law of the Skew Brownian motion with drift

We have the following proposition.

Proposition 2 We have that,

p β,µ (t, x, y) =                                                                            1 √ 2πt exp{µ(y -x) -1 2 µ 2 t} exp{-(y-x) 2 2t } -exp{-(y+x) 2 2t } + 1+β √ 2πt exp -(x+y) 2 2t + µ(y -x) -1 2 µ 2 t × 1 -βµ √ 2πt exp (x+y+tβµ) 2 2t N c ( βµt+x+y √ t ) , if x ≥ 0, y ≥ 0, 1-β √ 2πt exp -(x-y) 2 2t + µ(y -x) -1 2 µ 2 t × 1 -βµ √ 2πt exp (x-y+tβµ) 2 2t N c ( βµt+x-y √ t ) , if x ≥ 0, y < 0. 1 √ 2πt exp{µ(y -x) -1 2 µ 2 t} exp{-(y-x) 2 2t } -exp{-(y+x) 2 2t } + 1-β √ 2πt exp -(x+y) 2 2t + µ(y -x) -1 2 µ 2 t × 1 -βµ √ 2πt exp (-x-y+tβµ) 2 2t N c ( βµt-y-x √ t ) , if x < 0, y < 0, 1+β √ 2πt exp -(x-y) 2 2t + µ(y -x) -1 2 µ 2 t × 1 -βµ √ 2πt exp (y-x+tβµ) 2 2t N c ( βµt+y-x √ t ) , if x < 0, y ≥ 0. Remark 6 It can be shown that the quantity 1 -βµ √ 2πt exp (|x|+|y|+tβµ) 2 2t N c ( βµt+|x|+|y| √ t
) involved in p β,µ (t, x, y) remains strictly positive, whatever the sign of βµ (see Remark 7).

Proof of Proposition 2. We have

dB β,µ t = dW µ t + βdL 0 t (B β,µ ),
with W µ t = W t + µt a Brownian motion under Q µ defined by dQ µ dP = exp{-µW t -1 2 µ 2 t}. Note that under Q µ the process B β,µ starting from 0 is distributed as B β starting from 0 under P.

For any bounded continuous function f and any t ≥ 0, we have

E x P [f (B β,µ t )] = E 0 P [f (B β,µ t + x)] = E 0 Q µ [f (B β,µ t + x) exp{µW µ t -1 2 µ 2 t}] = R 2 f (y + x) exp{µw -1 2 µ 2 t}P 0 [B β t ∈ dy; W t ∈ dw] = R 2 f (y) exp{µw -1 2 µ 2 t}P x [B β t ∈ dy; W t -x ∈ dw] (11)
Suppose β > 0.

We set

Φ x (z, l) = (z, z -x -βl) which defines a bijection Φ x : R × R + → D x where D x = {(y, w) ∈ R 2 : y -x ≥ w}. Note that (B β t , W t -x) = Φ x (B β t , L 0 t (B β )). Besides, almost surely, (B β , L 0 (B β )) ∈ R × R + and (B β , W -x) ∈ D x .
For x > 0, Proposition 1 ensures that the measure P x [ B β t ∈ dy; L 0 t (B β ) ∈ dl ] has a density with respect to dy dl on R × R * ,+ , and gives mass to the segments of R + × {0} with the density

P x [ B β t ∈ dy; L 0 t (B β ) = 0 ] = 1 √ 2πt exp{- (y -x) 2 2t } -exp{- (y + x) 2 2t } dy. (12)
Let us denote

∆ x := {(y, w) ∈ R + × R : y = w + x} = Φ x (R + × {0}
). The measure P x [B β t ∈ dy; W tx ∈ dw] has a density g x B β ,W (y, w) with respect to dy dw on D x \ ∆ x = Φ x (R × R * ,+ ). But it gives mass to the segments of the line ∆ x . Let us denote Φ -1 x (y, w) := (Φ -1 1 (y, w), Φ -1 2 (y, w)) and notice that Φ -1 1 (y, w) = y.

Let A 1 ⊂ R + and A = {(y, w) ∈ R 2 : y ∈ A 1 , y = w + x} ⊂ ∆ x . As Φ -1 x (A) ⊂ R + × {0} we have P x [(B β t , W t -x) ∈ A] = P x [(B β t , L 0 t (B β )) ∈ Φ -1 x (A)] = P x [B β t ∈ Φ -1 1 (A); L 0 t (B β ) = 0 ] = P x [B β t ∈ A 1 ; L 0 t (B β ) = 0 ].
Using this and ( 12) in [START_REF] Etoré | On random walk simulation of one-dimensional diffusion processes with discontinuous coefficients[END_REF] we get

E x P [f (B β,µ t )] = Dx\∆x f (y) exp{µw -1 2 µ 2 t}P x [B β t ∈ dy; W t -x ∈ dw] + ∆x f (y) exp{µw -1 2 µ 2 t}P x [B β t ∈ dy; W t -x ∈ dw] = R f (y) y-x -∞ exp{µw -1 2 µ 2 t}g x B β ,W (y, w)dw dy + R + f (y) 1 √ 2πt exp{µ(y -x) -1 2 µ 2 t} exp{-(y-x) 2 2t } -exp{-(y+x) 2 2t } dy.
We now compute y-x -∞ exp{µw -1 2 µ 2 t}g x B β ,W (y, w)dw with a change of variable and an integration by parts. We have for y ≥ 0,

y-x -∞ exp{µw- 1 2 µ 2 t}g x B β ,W (y, w)dw = e -1 2 µ 2 t β y-x -∞ e µw (1 + β)( y-w-x β + x + y) √ 2πt 3 e - ( y-w-x β +x+y) 2 2t dw.
And,

y-x -∞ e µw ( y-w-x β + x + y)e - ( y-w-x β +x+y) 2 2t dw = βe µ(y-x) ∞ 0 e -βµw ′ (w ′ + x + y)e -(w ′ +x+y) 2 2t dw ′ = βe µ(y-x) te -(x+y) 2 2t -βµt ∞ 0 e -βµw ′ -(w ′ +x+y) 2 2t dw ′ = βe µ(y-x) te -(x+y) 2 2t - √ 2πβµt 3/2 e β 2 2 µ 2 t e βµ(x+y) N c ( x+y+tβµ √ t ) = β t e µ(y-x) e -(x+y) 2 2t 1 - √ 2πtβµe (x+y+βµt) 2 2t N c ( x+y+tβµ √ t )
which yields the desired result. The cases y < 0 and β < 0 are treated in a similar way.

For the case x < 0, we perform the change of variable x → -x, y → -y, β → -β and µ → -µ.

Exact simulation of bridges of a Skew motion with drift

For 0 < t < T let us denote q β,µ (t, T, a, b, y) the probability density of B β,µ t knowing that B β,µ 0 = a and B β,µ T = b. That is to say

P[B β,µ t ∈ dy | B β,µ 0 = a , B β,µ T = b] = q β,µ (t, T, a, b, y)dy.
Note that with these notations, q 0,0 (t, T, a, b, y) is the probability density of W t knowing that W 0 = a and W T = b. As the law of the Brownian bridge is well known, sampling from q 0,0 (t, T, a, b, y) is easy.

In order to sample along the law given by q β,µ (t, T, a, b, y) with a rejection algorithm using Brownian bridges values as proposals, we will use the two following results.

Lemma 3 Let a, b ∈ R, 0 < t < T . For (β, µ) ∈ (-1, 1) × R, we have ∀y ∈ R, q β,µ (t, T, a, b, y) = p β,µ (t, a, y)p β,µ (T -t, y, b) p β,µ (T, a, b) . ( 13
)
Proof. This comes from the Markov property for solutions of (2) (see Subsubsection 2.2.2 and the references therein).

Lemma 4

We have for all µ ∈ R, ∀0 < t < T, ∀a, b, y ∈ R, q 0,µ (t, T, a, b, y) = q 0,0 (t, T, a, b, y).

Proof. It is a direct consequence of Equations ( 5) and [START_REF] Faugeras | The Inverse EEG and MEG Problems : The Adjoint State Approach I: The Continuous Case[END_REF].

Lemma 5 Let a, b ∈ R, 0 < t < T .
-For (β, µ) ∈ (-1, 1) × R with βµ ≥ 0, we have

∀y ∈ R, q β,µ (t, T, a, b, y) ≤ K β,µ T,a,b q 0,0 (t, T, a, b, y), (14)
where

K β,µ T,a,b := 4α 2 p 0,µ (T, a, b) p β,µ (T, a, b)
,

and α := max( 1+β 2 , 1-β 2 ).
-For (β, µ) ∈ (-1, 1) × R with βµ < 0, set

γ β,µ (t, z) := 1 -βµ √ 2πt exp( (z + tβµ) 2 2t )N c ( βµt + z √ t ).
Then, p β,µ (t, x, y) ≤ 2αγ β,µ (t, |x|)p 0,µ (t, x, y)

and p β,µ (t, x, y) ≤ 2αγ β,µ (t, |y|)p 0,µ (t, x, y).

In particular, ∀y ∈ R, q β,µ (t, T, a, b, y) ≤ K β,µ t,T,a,b q 0,0 (t, T, a, b, y),

where

K β,µ t,T,a,b := 4α 2 γ β,µ (t, |a|)γ β,µ (T -t, |b|) p 0,µ (T, a, b) p β,µ (T, a, b) .
Proof. Case βµ ≥ 0. Let t > 0 and x ≥ 0. Looking at Proposition 2 it is clear that for y < 0, p β,µ (t, x, y) ≤ (1 + β)p 0,µ (t, x, y).

For y ≥ 0 we have

p β,µ (t, x, y) ≤ p 0,µ (t, x, y) + β √ 2πt exp -(x+y) 2 2t + µ(y -x) -1 2 µ 2 t ≤ (1 + β) p 0,µ (t, x, y),
where we have used (yx) 2 ≤ (y + x) 2 (because x, y > 0). We can proceed in a similar way for x < 0 and finally, we get that ∀t > 0, ∀x, y ∈ R, p β,µ (t, x, y) ≤ 2α p 0,µ (t, x, y).

Thus, using the previous inequality gives

q β,µ (t, T, a, b, y) = p β,µ (t, a, y)p β,µ (T -t, y, b) p β,µ (T, a, b) ≤ 4α 2 p 0,µ (T, a, b) p β,µ (T, a, b) p 0,µ (t, a, y)p 0,µ (T -t, y, b) p 0,µ (T, a, b) ≤ 4α 2 p 0,µ (T, a, b) p β,µ (T, a, b) q 0,0 (t, T, a, b, y). ( 19
)
Case βµ < 0. Let us denote Γ β,µ (t, x, y) x,y) is an even function. As we have

:= 1 -βµ √ 2πt exp (|x|+|y|+tβµ) 2 2t N c ( βµt+|x|+|y| √ t ). For fixed x ∈ R, y → Γ β,µ (t,
∀z > 0, z e z 2 2 ∞ z e -u 2 2 du < 1, (20)
the function z → √ 2π exp( z 2 2 )N c (z) has negative first derivative on R + . Therefore y → Γ β,µ (t, x, y) is decreasing on R + and we have max y∈R Γ β,µ (t, x, y) = γ β,µ (t, |x|). Using this and the same kind of computations than in the previous case we get [START_REF] Kulik | On the solution of a one-dimensional stochastic differential equation with a singular drift coefficient[END_REF]. As the roles of x and y are symmetric in Γ β,µ (t, x, y), we get [START_REF] Gall | One-dimensional stochastic differential equations involving the local times of the unknown process[END_REF]. We then obtain [START_REF] Lejay | On the decomposition of excursions measures of processes whose generators have diffusion coefficients discontinuous at one point. Markov Process[END_REF], using the same computations than for [START_REF] Lejay | Simulating a diffusion on a graph. Application to reservoir engineering[END_REF].

Remark 7 Note that (20) also allows to prove that p β,µ (t, x, y) remains strictly positive (see Remark 6).

Lemma 5 suggests to use the following rejection algorithm in order to sample along q β,µ (t, T, a, b, y)dy (see for example Proposition 1 in [START_REF] Beskos | Exact simulation of diffusions[END_REF]).

EXACT SIMULATION ALGORITHM ALONG q β,µ (t, T, a, b, y)dy

1. Set K = K β,µ T,a,b if βµ ≥ 0, K = K β,µ t,T,a,b otherwise. 2. Sample Y along q 0,0 (t, T a, b, y)dy. 3. Evaluate f (Y ) := 1 K q β,µ (t, T, a, b, Y ) q 0,0 (t, T, a, b, Y ) ≤ 1. 4. Simulate U ∼ U ([0, 1]). If U ≤ f (Y ) accept the proposed value Y . Else return to Step 2.

Exact Algorithm Bridges

Acceptance Ratio 0.28 0.18

Numerical experiments

Example 1. We first deal with a toy example. We consider the following SDE

dX t = dW t - π 2 cos( π 5 X t )dt + βdL 0 t (X), X 0 = x 0 , (21)
with β = 0.6, and x 0 = 0.2. Note that, here, the drift b(x) = -π 2 cos( π 5 x) is bounded and of class C ∞ on the whole real line.

For the exact procedure the constant drift involved in Subsection 3.2, equals µ = b(0) = -π 2 . So we will have to sample bridges of SBM with non zero drift µ, using the results of Section 5.

We have first to sample X T from

h(y) = C exp (B(y) -B(x 0 )) p β,µ (T, x 0 , y) = C exp 5 2 sin( π 5 x 0 ) -sin( π 5
y)µ(yx 0 ) p β,µ (T, x 0 , y) (Step 1 of the Exact Algorithm). This can be done by rejecting standard normal random variables with mean x 0 . Indeed, using [START_REF] Kulik | On the solution of a one-dimensional stochastic differential equation with a singular drift coefficient[END_REF], we have here

h(y) C ≤ 2αγ β,µ (T, |x 0 |) exp(5 - µ 2 T 2 ) p 0,0 (T, x 0 , y).
Then we accept or reject the proposed value X T , using Steps 2 to 4 of the Exact Algorithm, with bridges of B β,µ ,

φ(x) = π 2 8 cos 2 ( π 5 x) + π 2 20 sin( π 5 x) + π 2 20 ,
and K = 9π 2 20 as un upper bound for φ. We plot on Figure 1 (top and bottom figures) the approximated density obtained with 10 6 simulations of X T , sampled with our exact procedure. On the top figure we plot the approximated densities obtained with 10 6 simulations of the Euler Scheme used in [START_REF] Martinez | Interprétations probabilistes d'opérateurs sous forme divergence et analyse des méthodes numériques probabilistes associées[END_REF] and [START_REF] Martinez | Discrétisation d'équations différentielles stochastiques unidimensionnelles à générateur sous forme divergence avec coefficient discontinu[END_REF], for decreasing time steps. We can observe the convergence of Euler type simulations to exact ones. Note that to have the Euler scheme fitting the exact procedure we have to take a fine time step (namely ∆t = 10 -4 ). This is because, as shown in [START_REF] Martinez | Interprétations probabilistes d'opérateurs sous forme divergence et analyse des méthodes numériques probabilistes associées[END_REF], the rate of weak convergence of the Euler scheme in this situation is of order (∆t) 1/2-ǫ , for a smooth initial condition.

On the bottom figure the approximated density is compared with the approximated densities obtained with 10 6 simulations of the random walk based method studied in [START_REF] Etoré | A Donsker theorem to simulate one-dimensional processes with measurable coefficients[END_REF], for decreasing space steps. Again we can observe the convergence of the process with discretization error.

In Table 1 we report the empirical acceptance ratios for the rejection step using φ and the Poisson point process in the Exact Algorithm (this corresponds to the column Exact Algorithm in the table), and for the rejection sampling of bridges of the SBM with drift (this is the average acceptance ratio in this case).

In Table 2 we report the CPU times needed to get the 10 6 simulations, with the three different methods (and with the different discretization steps we have used). Programs were written in Clanguage and executed on a personal computer equipped with an Intel Core 2 duo processor, running at 2.23 Ghz.

Exact

Euler Random Walk (∆t = 10 -n , n = 2, 4)

(h = 1 10 , 1 200 ) 239s 17s 3

.52s 1680s 1411s

Table 2: CPU times for 10 6 simulations of XT .

On this example the exact simulation is competitive, compared to schemes with very fine grids.

Example 2. We want now to sample along the law of the continuous Markov process X generated by

L = 1 2 d dx a d dx • (22)
with

a(x) =      x 2 +x+1 (2x+1) 2 if x ≥ 0 3x 2 -x+2 (6x-1) 2 if x < 0.
Note that a(0+) = 1 = 2 = a(0-). The coefficient a(x) is of class C 1 on R * ,-and R * ,+ , and uniformly strictly positive and bounded, which ensures the existence of X; in addition X solves

dX t = a(X t )dW t + a ′ (X t ) 2 dt + a(0+) -a(0-) a(0+) + a(0-) dL 0 t (X), (23)
(see [START_REF] Lejay | A scheme for simulating one-dimensional diffusion processes with discontinuous coefficients[END_REF], [START_REF] Etoré | On random walk simulation of one-dimensional diffusion processes with discontinuous coefficients[END_REF]). We define the Lamperti transformation Φ(x) =

x 0 dz/ a(z) and set Y t := Φ(X t ). Then

dY t = dW t + 1 2 ( √ a) ′ • Φ -1 (Y t )dt + a(0+) -a(0-) a(0+) + a(0-) dL 0 t (Y ), (24)
(this follows from Proposition 3.1 in [START_REF] Etoré | On random walk simulation of one-dimensional diffusion processes with discontinuous coefficients[END_REF]; see also [START_REF] Lejay | A scheme for simulating one-dimensional diffusion processes with discontinuous coefficients[END_REF] and [START_REF] Ouknine | Skew-Brownian motion" et les processus qui en dérivent[END_REF]). Firstly, note that √

a(0+)- √ a(0-) √ a(0+)+ √ a(0-) < 1.
Secondly, we have

( √ a) ′ (x) =      1 2 √ x 2 +x+1 -2 √ x 2 +x+1 (2x+1) 2 if x ≥ 0 - 1 2 √ 3x 2 -x+1 + 6 √ 3x 2 -x+2 (6x-1) 2 if x < 0, Φ(x) =    2 √ x 2 + x + 1 -2 if x ≥ 0 -2 √ 3x 2 -x + 1 + 2 √ 2 if x < 0, and
Φ -1 (y) =        -1+ √ (y+2) 2 -3 2 if y ≥ 0 1- √ 1-12[2-( √ 2-y/2) 2 ] 6 if y < 0.
As ( √ a) ′ (x) is bounded with bounded first derivative on R * ,-and R * ,+ , the explicitly known

coefficients β = √ a(0+)- √ a(0-) √ a(0+)+ √ a(0-)
and b(y) = 1 2 ( √ a) ′ • Φ -1 (y) satisfy the assumptions of Subsubsection 3.1.2. Thus we can perform exact sampling from [START_REF] Martinez | One-dimensional parabolic diffraction equations: pointwise estimates and discretization of related stochastic differential equations with weighted local times[END_REF], and, applying the exact inverse transformation Φ -1 , get samples from [START_REF] Martinez | Discrétisation d'équations différentielles stochastiques unidimensionnelles à générateur sous forme divergence avec coefficient discontinu[END_REF] with absolutely no discretization error.

Here we have,

µ = 1 4 a ′ (0+) -a ′ (0-) a(0+) -a(0-) = - 26 4(1 - √ 2
) . As we have

B(y) =    -µy + 1 2 log( √ a • Φ -1 (y)) if y ≥ 0 -µy + 1 2 [log( √ a • Φ -1 (y)) -log( √ 2)] if y < 0,
we can show (using again [START_REF] Kulik | On the solution of a one-dimensional stochastic differential equation with a singular drift coefficient[END_REF]), that for all y 0 , y ∈ R and T > 0, exp (B(y) -B(y 0 )) p β,µ (T, y 0 , y) ≤ √ 2 √ 24 2αγ β,µ (T, |y 0 |) e -1 2 µ 2 T p 0,0 (T, y 0 , y).

This allows to sample Y T from h(y) = C exp (B(y) -B(y 0 )) p β,µ (T, y 0 , y), by rejecting normal variables with mean y 0 and variance T .

We then accept or reject the proposed value Y T by using bridges of B β,µ and

φ(y) = ((1/2)( √ a) ′ • Φ -1 (y)) 2 + (1/2)(( √ a) ′′ √ a) • Φ -1 (y) 2 , with ( √ a) ′′ (x)        - 2x+1 4(x 2 +x+1) 3/2 - 1 (2x+1) √ x 2 +x+1 + 8 √ x 2 +x+1 (2x+1) 3 if x ≥ 0 6x-1 4(3x 2 -x+2) 3/2 - 3 (2x+1) √ 3x 2 -x+2 + 72 √ 3x 2 -x+2 (6x-1) 3 if x < 0.
We take K = (6

√ 2-1/2) 2 /4+(141-1/8)/2 2
as an upper bound for φ. We plot on Figure 2 the approximated density computed with 10 7 simulations of X T for x 0 = 0.0 and T = 1, obtained from the exact procedure. We plot on the same figure the approximated densities obtained with the Euler scheme and the random walk approximation mentioned in Example 1.

We report in Table 3 the acceptance ratios.

Exact Algorithm Bridges Acceptance Ratio 0.017 0.5 Remark 8 Note that the acceptance ratio for the algorithm in the first example in Table 1 is quite low but decreases to less than 2% in Table 3 in the context of the second example. These figures are closely related to the measurement of the "distance" between the measure of the initial process from the reference measure and so these limitations of the algorithm arise even in the "classical" setting of the reference article [START_REF] Beskos | Retrospective exact simulation of diffusion sample paths with applications[END_REF] (for example with a rapidly varying drift). Nevertheless, in terms of CPU time, the performance of the algorithm seems quite competitive, in comparison with those of discretization schemes.

Remark 9 Note that, at least graphically and contrary to what we can see on Figure 1, the transition density plotted on Figure 2 seems to be continuous at 0 : this matches the well-known theoretical result, which asserts that the transition density of diffusion semigroups corresponding to elliptic divergence form operator of the form ( 22) is always continuous. We refer to Stroock [START_REF] Stroock | Diffusion semigroups corresponding to uniformly elliptic divergence form operators[END_REF] for a proof based on the self-adjoint properties of these semi-groups and Nash's inequality.

7. Discussion and concluding remarks 7.

1. An open problem : the path decomposition of a skew Brownian bridge An important issue for the extension of the initial exact simulation method is to overcome the restraining assumptions made on the drift function b (see Section 2.1) : namely, the assumption of boundedness for b.

For example, it is frustrating that these assumptions do not allow us to simulate exactly what one may call the "Skewed Ornstein-Uhlenbeck" diffusion process. This difficulty appears even in the classical case (solutions of non skewed SDEs) and the fundamental reason is that we do not know how to simulate exactly a Poisson Point Process with σ-finite intensity on the whole space R + × R.

In the classical case, where b is everywhere differentiable and no local time is involved (non skewed SDEs), this problem is solved by decomposing the trajectory of the standard Brownian bridge on [0, T ] w.r.t. the space-time point where it attains its maximum or both its maximum and its minimum : we refer to [START_REF] Beskos | A factorisation of diffusion measure and finite sample path constructions[END_REF] for a detailed presentation of this problem in the classical setting.

Consequently, if one wants to overcome the restraining assumptions in Section 2.1 concerning the drift function b, one has to search for such kind of decompositions for (at least) the Skew Brownian Motion (not to mention the drifted Skew Brownian Motion). Up to our knowledge, no results can be found in the literature concerning this decomposition and this open problem seems difficult to us. However, we give below some insight concerning this problem thanks to an application of a theoretical result stated in [START_REF] Pitman | Decomposition at the maximum for excursions and bridges of onedimensional diffusions[END_REF].

Let τ β,µ Proposition 3 (case µ = 0)

In the simple case where µ = 0, the function u λ is given by

u λ (x; z) =                                                                  sinh √ 2λ(z -x) sinh( √ 2λz) 1 + β 2 cosh( √ 2λz) -(1 -β)e - √ 2λz + sinh( √ 2λx) sinh( √ 2λz) if z ≥ x > 0, e - √ 2λ(x-z) if x ≥ z ≥ 0, e √ 2λx 1 + β 2 cosh( √ 2λz) -(1 -β)e - √ 2λz if x < 0 < z, e √ 2λ(x-z) if 0 ≥ z ≥ x, e - √ 2λx 1 -β 2 cosh( √ 2λz) -(1 + β)e √ 2λz if z < 0 < x, sinh √ 2λ(z -x) sinh( √ 2λz) 1 -β 2 cosh( √ 2λz) -(1 + β)e √ 2λz + sinh( √ 2λx) sinh( √ 2λz) if 0 > x ≥ z. (25)
Remark 10 Note that if β = 0, we retrieve after easy computations the well known result that gives the Laplace Transform of the law of the hitting time of z by a standard Brownian Motion starting from x.

Proof. We only sketch the proof. The different cases may be easily conjectured from the description of the excursion measure for the SBM (B β,0 s ) s≥0 and the known facts concerning the standard Brownian Motion (decomposition of the different cases when a skew Brownian Motion reaches z starting from x). In order to check rigorously the validity of the result, one may verify that the formulas (25) yield a solution of Dynkin's problem associated to the generator of (B β,0 s ) s≥0 namely : (see [START_REF] Lejay | On the constructions of the skew Brownian motion[END_REF]). In particular, the density ℓ β,0 (t, x, y)dy of the SBM w.r.t. the speed measure m(dy) is given by ℓ β,0 (t, x, y) = (of course ℓ β,0 (t, x, y) = ℓ β,0 (t, y, x)). Let M β,0 T := sup 0≤s≤t B β,0 s ; ρ β,0 T := inf{s ≥ 0 : B β,0 s = M β,0 s }.

   1 2 d 2 dx 2 u λ (.; z) = λu λ (.; z) u λ (z; z) = 1, (26)
                               1 (1 + β) √ 2πt exp{- (y -x)
Then, applying the results of Theorem 2 in Pitman-Yor [START_REF] Pitman | Decomposition at the maximum for excursions and bridges of onedimensional diffusions[END_REF], we have the following proposition : An open problem is to find a description of these laws and to give a procedure in order to simulate these laws exactly.

Concluding remarks

In this paper we presented an extension of the exact simulation method of [START_REF] Beskos | Retrospective exact simulation of diffusion sample paths with applications[END_REF] that permits to sample an exact skeleton of a one dimensional diffusion process skewed at 0. This method may be applied to diffusions related to strongly elliptic divergence form operators that possess a discontinuous coefficient at 0. The basic idea of this contribution depends highly on the possibility to perform a Girsanov transformation such that no local time appears in the Girsanov exponential weight and such that the reference measure is tractable.

In our opinion, this first work should be extended in several directions. Firstly, it is necessary to give a complete treatment of the case β = 0 (see Remark 2). In this case, there still exists a way to perform a Girsanov transformation such that no local time appears in the Girsanov exponential weight, but then the reference measure becomes that of a Brownian motion with two-valued drift (see [START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF] for an introduction to these particular types of Brownian motions). As before, the difficulty arises for the simulation of the bridges in the Step 3 of the algorithm.

Further digging shows that, in this particular situation, the solution of the exact simulation problem in the manner of [START_REF] Beskos | Retrospective exact simulation of diffusion sample paths with applications[END_REF] is closely related to the computation of joint laws for the position together with local and occupation times by an arbitrary Brownian bridge (with no drift) but conditioned on its final position and local time at 0. Even if there exists abundant litterature dealing with the Brownian bridge, there is no result for such joint laws.

Consequently, we believe that performing a totally exact simulation algorithm in full generality in the case β = 0 appears to be outside the scope of this paper. Note that a satisfactory treatment of the case β = 0 is crucial if one has the objective to deal with the even more general case, where the discontinuity of b and the local time appear at distinct space points.

Secondly, various questions arise in the treatment of "skewness" : how can we overcome the restraining boundedness assumption on the drift function b ? What about a one dimensional diffusion process skewed at a finite number of points ?

  z→0+ b(z) := b(0+) and lim z→0-b(z) := b(0-) exist and are finite. The value b(0) of the function b at 0 is of no importance and can be fixed arbitrarily to some constant (possibly different from either b(0+) or b(0-)).

SDEs of type ( 2 ) 2 . 2 . 1 .

 2221 Existence and uniqueness Under the assumptions of the previous section, the equation (2) possesses a unique strong solution. In fact, performing the bijective change of variable g(x) → (1β)x½ x≥0 + (1 + β)x½ x<0 allows to consider Y t := g(X t ) solution of a new transformed equation

3. 1 .

 1 Notations and additional assumptions 3.1.1. Notations

Figure 1 :

 1 Figure 1: Approximated densities of the positions at time T = 1.0 of 10 6 paths of the solution of (21) starting from x0 = 0.2: exact versus Euler with time step ∆t = 10 -n , for n = 2, 4 (top) and exact versus random walk with space steps h = 1 10 , 1 200 (bottom).

Figure 2 :

 2 Figure 2: DIVERGENCE FORM OPERATOR: Approximated density of the positions at time T = 1 of 10 7 paths of the solution of (23) starting from x0 = 0.0: exact versus random walk with space step h = 3.10 -3 and Euler scheme with ∆t = 10 -4 .

z:

  = inf(s ≥ 0 : B β,µ s = z). Set u λ (x; z) := E x e -λτ β,µ z which gives the Laplace transform of τ β,µ z at λ > 0 (with B β,µ 0 = x).

  withu λ (.; z) ∈ {g ∈ C 0 (R) ∩ C 2 ((-∞, 0) ∪ (0, ∞)) : (1 + β)g ′ (0+) = (1β)g ′ (0-)}.A scale function s and the corresponding integrated speed measure m of a Skew Brownian Motion are given by s(x) = 2 β+1 x if x ≥ 0

2 1 -

 1 β x if x < 0 ; m(x) = (β + 1)x if x ≥ 0 (1β)x if x < 0.

Proposition 4 1 . 27 ) 2 .

 1272 For any a, b ≤ z < ∞, λ > 0, we have thatE e -λρ β,0 T ½ M β,0 T ∈dz | B β,0 0 = a, B β,0 T = b = u λ (a; z)u λ (z; b) ℓ β,0 (T, a, b) s(dz). (Moreover, under P . | B β,0 0 = a, B β,0 T = b, M β,0 T = z, ρ β,0 T = u ,the path fragments B β,0 s : 0 ≤ s ≤ u B β,0 T -s : 0 ≤ s ≤ Tu are independent, distributed respectively like B β,0 s : 0 ≤ s ≤ τ β,0 z under P . | B β,0 0 = a given τ β,s ≤ τ β,0 z under P . | B β,0 0 = b given τ β,0 z = Tu.

Table 1 :

 1 Acceptance ratios in Example 1.

Table 3 :

 3 Acceptance ratios in Example 2.

Acknowledgements

We would like to thank the anonymous referee for helpful comments and having pointed the reference [31].


Virgile Bocchetti
email: virgile.bocchetti@u-cergy.fr

H T Diep
email: diep@u-cergy.fr

Monte Carlo Simulation of Melting and Lattice Relaxation of the (111) Surface of Silver

Keywords: numbers:05.10.Ln, 64.70.dj, 05.70.Np : Structure

It is experimentally observed and theoretically proved that the distance between topmost layers of a metal surface has a contraction. However, well-known potentials such as Lennard-Jones and Morse potentials lead to an expansion of the surface inter-layer distance. Such simple potentials therefore cannot be used to study metal surface relaxation. In this paper, extensive Monte Carlo simulations are used to study the silver (111) surface with both the Gupta potential (GP) and the Embedded Atom Method (EAM) potential. Our results of the lattice relaxation at the (111) surface of silver show indeed a contraction for both potentials at low temperatures in agreement with experiments and early theories. However at higher temperatures, the EAM potential yields a surface melting at ≃ 700 K very low with respect to the experimental bulk melting at ≃ 1235 K while the GP yields a surface melting at ≃ 1000 K closer to the bulk one. In addition, we observe with the EAM potential an anomalous thermal expansion, i. e. the surface contraction becomes a surface dilatation with respect to the bulk, at ≃ 900 K. The Gupta potential does not show this behavior.We compare our results with different experimental and numerical results.

I. INTRODUCTION

It is well known that surface atoms of a material have a behavior different from that of bulk atoms, mainly because of the lack of neighbors and the surface geometry. Different kinds of surface behavior can occur according to the nature of the material and the surface orientation. The case of metallic materials has been well studied theoretically [1][2][3][4][5][6][7] and by means of different experimental techniques [8][9][10].

The contraction of the lattice spacing is experimentally observed by Medium Energy Ion Backscattering (MEIS) [8], Low-Energy Electron Diffraction (LEED) [9] or X-ray scattering [11]. Other methods such as elastic He scattering and electron energy-loss spectroscopy have also been used. In a theoretical point of view, Gupta has shown analytically that for classical pairwise potentials (Lennard-Jones, Morse, ...), inter-layer distance near the surface have an expansion. He has also shown that the Tight-Binding-Potential (TBP) and the so-called Gupta potential (GP) lead to a contraction of inter-layer distance at metallic surfaces.

Despite the important number of experimental and theoretical techniques used to observe this phenomenon, there exists an unsolved question on how the contraction evolves with increasing temperature. There are two contradictory answers in the literature: X-ray scattering [11] and LEED [9] as well as molecular dynamics (MD) simulations using the Embedded Atom Method (EAM) [7,[START_REF] Ahlam | Rahman[END_REF] show the surface inter-layer distance always smaller than the bulk one at the same temperature, namely surface is contracted, whereas MEIS experiments [8] and ab-initio density-functional theory (DFT) calculations [13] show an anomalous thermal expansion of the surface at some temperature below the bulk melt. Facing this long-standing unsolved question, we wanted to carry out a Monte Carlo (MC) study in an attempt to answer that question. To our knowledge, there is no MC simulations in literature so far about this subject, although some MC simulations have been used to reproduce experimental patterns such as surface blocking pattern during scattering process [see Ref. [10] for the Pb (110) surface]. Given a tremendous number of numerical studies on surface problems, it is surprising that no MC simulation has been performed so far to see the variation of the surface inter-layer distance.

The purpose of this paper is thus to investigate by MC simulation the variation of the lattice spacing between the topmost layers of the (111) silver surface versus temperature.

In order to simulate such a behavior as accurately as possible, we have considered potentials which describe as well as possible the material. The EAM potential is often used in MD simulations and especially for the Ag (111) surface [7,[START_REF] Ahlam | Rahman[END_REF]14]. Working with this potential allows us to compare our results with other numerical studies using the same potential. Furthermore, the EAM potential reproduces accurately the bulk melting temperature of Ag. On the other hand, the Gupta potential describes well surface and cluster behaviors [15][16][17][START_REF] Rogan | Ricardo[END_REF][19]. The melting temperature of bulk Ag is also well reproduced with this potential. That was the reason why the two potentials GP and EAM have been used since many years to simulate silver material and other metallic crystals. However, as will be seen below, the two potentials, although yielding the same result for low-temperature surface contraction, give different results at higher temperature concerning the surface contraction and the surface melting.

In Section II we recall essential properties of the two potentials with their sets of parameters and we briefly describe our algorithm. The results obtained by our computation are shown in Section III. Concluding remarks are given in Section IV.

II. MODEL AND MONTE CARLO METHOD

A. Potentials

Gupta potential

In computer simulation of metals, the Gupta potential is one of the most used semi-empirical potentials. There is a multiple reason for this success but the main ones are the accuracy of its results for metals and its easy and quick implementation. The Gupta potential is derived from the Gupta's expression of the cohesive energy of the bulk material and is based on the second-moment approximation of the electron density of states in the tight-binding theory. It includes implicitly some manybody interactions. The expression of the potential is the following:

V (r) = E n N i   A n i =j e -p r i,j r 0 -1 - N i =j e -2q r i,j r 0 -1   (1)
where A is a constant given in eV, r 0 the equilibrium nearest-neighbor (NN) distance in the bulk metal, p the repulsive interaction range and q the attractive one. r i,j is the distance between the atoms i and j and E n an energy constant which depends on the size of the system. In Eq. ( 1), the first part gives is a Born-Mayer pairwise repulsion energy term and the second part is the manybody attractive contribution.

The parameters used in this work are reported in Table I.

EAM potential

Several authors [21][22][23] have proposed a method based on density-functional ideas called EAM. We used a version of the EAM potential given in Ref. [23]. The parameters of the EAM potential used in this work for silver are reported in Table II. In the EAM potential, the total potential energy is given by:

E p = i   F i (ρ i ) + 1 2 j =i φ (r ij )   (2)
where φ ij represents the pair energy between two atoms i and j at the distance r ij . F i is the embedding energy function which represents the energy to embed an atom into a local site with electron density ρ i . The electron density ρ i has the following expression:

ρ i = j =i f (r ij ) (3)
with

f (r ij ) = f e exp -β rij re -1 1 + rij re -λ 20 (4)
where r e , f e , β and λ are constant parameters which are given in Table II. The embedding function has the following form:

F (ρ) =          3 i=0 F ni ρ ρn -1 i ρ < 0.85ρ e , 3 i=0 F i ρ ρe -1 i 0.85ρ e ≤ ρ < 1.15ρ e , F e 1 -ln ρ ρe η ρ ρe η 1.15ρ e ≤ ρ          (5)
The pair energy expression between atoms i and j is :

φ (r ij ) = A exp -α rij re -1 1 + rij re -κ 20 - B exp -β rij re -1 1 + rij re -λ 20 (6)
All constant parameters in the above expressions such as F ni etc. are listed in Table II. Note that MD simulations using EAM potential yield the bulk melting temperature at 1170 K for Ag [22] which is to be compared to the experimental value of 1235 K [START_REF] Touloukian | Thermophysical Properties of Matter[END_REF].

B. The algorithm

The complete description of the algorithm can be found in Ref. [START_REF] Bocchetti | [END_REF]. For the present surface problem, we have a Reference [23] used the the periodic boundary conditions (PBC) in two directions x and y. Let us briefly recall the main steps of the algorithm. The step that requires a long CPU time is the construction of the PBC table which contains neighbors of each atom of the simulation sample. The reader is referred to Ref. [START_REF] Bocchetti | [END_REF] for details. Note that to apply correctly the PBC we have to take into account instantaneous fluctuations of the sample size in each direction at every MC step.

Using the least square approximation for experimental data of the NN distance measured for a bulk sample taken from Ref. [START_REF] Touloukian | Thermophysical Properties of Matter[END_REF] we obtain the linear thermal expansion d N N ( Å) = α L T + 2.87053 with the coefficient α L = 6.27617 × 10 -5 . In Fig. 1 we have plotted experimental data taken from Ref. [START_REF] Touloukian | Thermophysical Properties of Matter[END_REF] and the least square approximation used in this paper. The lattice constant at a given temperature is now an input data. In our calculations, we begin at T = 200 K in order to stay in the classical limit. Indeed, as we can see in Ref. [26] that at lower T , the thermal expansion coefficient does not have the same linear behavior below and above 200 K. At each temperature, we start the simulation with the lattice constant taken from Fig. 1. In doing so, we favor a fast convergence to equilibrium. This was what we observed in test runs. Only the surface atoms and an adjustable number of layers beneath the surface are moving. It allows us to choose a number of moving layers to produce precise results while saving CPU time. Of course, we have to make sure that final results do not depend on the choice of the number of moving atoms (see an example in Fig. 2).

Note that by using the correct value of the lattice constant at each temperature for interior layers far from the surface, we take into account some temperature effects but we neglect local fluctuations in those far layers. The effect of such approximation has been tested: the result confirms that this does not affect our conclusion. As seen below, in most simulations which have been made, atoms are allowed to move in the first three layers. We have observed that the contraction of the surface lattice spacing disappears rapidly at a few layers from the sur- face. This approximation is justified by tests we have made and by experimental observations and theoretical calculations mentioned above. Note that we use the Metropolis updating criterion. At a given temperature, we equilibrate the system before computing the averages of quantities needed for analyzing the system behavior as functions of temperature.

III. RESULTS AND DISCUSSIONS

In this section we describe our MC results for the (111) surface of silver. We have studied a system of 256 atoms per layer in which atoms in the first three layers are allowed to relax at each MC step. From the fourth layer inward, the atoms dilate uniformly with temperature using data of Fig. 1. For comparison, we have also considered a system of 100 atoms per layer with 3 moving layers. As it turned out, the thermal relaxation of the first moving layers is not significantly different for these two surface sizes. Therefore, we show below our results obtained with the surface size of 256 atoms on the top of five beneath layers (6 layers in all) with the PBC in the x and y directions. Usually, our simulations are carried out over 5 to 6 millions of MC steps per atoms. We discarded the two first millions during equilibrating time and averaged physical quantities over the following 2-3 millions of MC steps.

A. Surface melting

The first step is the monitoring of the surface stability versus temperature. Of course, our calculations have been carried out for temperatures below the bulk melting temperature (T m ) because we use the temperaturedependent bulk lattice spacing below T m shown in Fig. 1 for bulk atoms. Note however that very close to T m , the thermal coefficient does not follow the linear behavior. We have therefore heated the system from the ground state to temperatures between 200 K and 1200 K, just below T m , with a step of 50 K. The mean displacement amplitude allows us to establish the validity domain of temperatures. This quantity is computed as follows : at each MC step, we compute a spatial average of all separation distances between atoms with their ground-state node positions and then we compute its thermal average through MC steps. As we can see in Fig. 3, with the Gupta potential, after 1000 K the atoms undergo a sudden jump indicating a phase change: the surface becomes disordered (melted) between 1000 K and 1050 K. This value from our MC simulation using Gupta potential is very close to the value 1100 K obtained from MD simulation using the EAM potential [START_REF] Ahlam | Rahman[END_REF]. A study is needed to determine the precise surface melting temperature but it is not the aim to search for such a precision of this paper.

For the EAM potential, we have obtained a result different from that of MD simulation: the surface melting temperature takes place at about 700 K, while MD simulation fonds it at 1100 K as we mentioned above. In order to see clearly this surface melting, we have com-puted with the EAM potential the structure factor for the topmost layers. The structure factor S K is computed as follow:

S K = 1 N l N l j=1 e i - → K . - → dj (7)
where -→ d j is the position vector of an atom in the layer, N l the number of atoms in a layer and -→ K the reciprocal lattice vector which has the following coordinate (in reduced units): 2π 1; -√ 3; 0 . The angular brackets < ... > indicate thermal average taken over MC run time. The above "order parameter", which allows us to monitor the long-range surface order, is plotted for the surface layer in Fig. 4. As we can see, the long-range order is lost at ≃ 700 K. In order to investigate in more details the surface melting, we have also computed the O 6 order parameter introduced in Ref. [27]. This order parameter describes the short-range hexatic orientational order of the surface which is defined as follows:

O 6 = jk W jk e i6Θ jk jk W jk (8)
with

W jk = e -(z j -z k ) 2 2δ 2 (9)
where the sum runs over the NN pairs and Θ jk is the angle which the j -k bond, when projected on the xy plane, forms with the x axis. The δ parameter is taken as onehalf the average inter-layer spacing. The weighting function, W jk , allows us to differentiate the "non coplanar" and the "coplanar" neighbors. With a coplanar neighbor, the weighting function takes a maximum value. We have to calculate the spatial average of O 6 taken over all atoms of the surface layer and then calculate its thermal average over MC run time. We plot the so-averaged O 6 parameter and the S K structure factor versus temperature in Fig. 4. We observe that the long-range and short-range orders are lost at the same temperature: the melting of the surface occurs at around 700 K where both order parameter change their curvature. Snapshots of the system at two temperatures close to the surface melting are shown in Figs. 5 and6, for the two potentials.

The difference in the surface melting temperature for the two potentials is surprising since both give the same bulk melting temperature. One possible explanation comes from the construction of their set of parameters: the Gupta potential was fitted in Ref. [20] with the bulk melting temperature, while the EAM potential is fitted to reproduce physical quantities such as the lattice constant, the elastic constants C 11 , C 12 and C 44 and sublimation energy, etc. but not with the bulk melting temperature.

A consequence is that surface atoms in the Gupta case are more strongly attached to the interior part, so that we need a much higher temperature to make the surface melt. We believe that the high "stability" of surface atoms is the reason why the Gupta potential was very suitable for the study of clusters of very few number of atoms [15].

For both potentials, a contraction of the lattice spacing between two topmost layers occurs but its temperature dependence is different in the two cases, as will be shown below.

B. Surface contraction

In order to see the variation of the inter-layer distance at the surface, we have plotted in Fig. 7 the z-position distribution obtained for the two potentials EAM and GP. As we can see, when the temperature increases, the three peaks which represent the z coordinates of the three moving layers, are shifted to the left, indicating a contraction toward the bulk. Note that the two potentials lead to the same behavior.

In these two z-axis distribution functions, we can see that at a given T , the deviation of the topmost layer from its perfect position is largest. This deviation becomes smaller for the next two layers and it should disappear at a few layers from the surface. This observation justifies our approximation to let only the three topmost layers to relax. Previous theoretical results and experimental data also give support to our hypothesis. In order to examine closely the inter-layer contraction, we have plotted in Fig. 8 the inter-layer distance ∆ 12 between layer 1 and 2, at different temperatures. This quantity is determined by averaging over the corresponding peaks during MC simulation time. We only focus our attention on ∆ 12 since for this quantity there is a general agreement between experiments and between experiments and theories. For inner layers, there is no such agreement between experiments and theories, and even between experiments, depending on the technique used: LEED [9], MEIS [8], ... Let us comment the results shown in Fig. 8. Several remarks are in order:

(i) At low temperatures, for both potentials, ∆ 12 has a contraction of about 2.5% at 300 K, compared to a contraction of 0.5% with MD simulations using EAM potential [7,[START_REF] Ahlam | Rahman[END_REF]. Experimental contraction of 2.5% was observed in Ref. [8].

(ii) The distance ∆ 12 increases with increasing temperature: in the case of EAM potential, ∆ 12 crosses the bulk limit at ≃ 900 K, indicating a surface expansion, in agreement with experimental results in Ref. [8]. This "expansion" was not found by MD simulations with EAM potential (see Refs. [7,[START_REF] Ahlam | Rahman[END_REF]) up to 1100 K and in LEED [9]. A surface expansion was observed in an ab-initio DFT calculation [13].

(iii) As we can see in Fig. 8, MC results with GP does not show the above-mentioned anomalous expansion up to 1150 K slightly above the surface melting temperature observed above with this potential. We will give below an explanation about this point at least from our MC results. As said above, the contraction of the first inter-layer distance ∆ 12 has been observed by both MEIS and LEED at low temperatures. For ∆ 23 , the distance between the second and third layers, these two methods are not in agreement: at low T , MEIS finds a little expansion of ∆ 23 while LEED finds a contraction. Here, our computations are in agreement with LEED results.

Let us discuss about the controversial point on the anomalous expansion of ∆ 12 , namely the transition from contraction to expansion when ∆ 12 crosses the bulk line (see Fig. 8). As said earlier, this has been experimentally observed by MEIS at about 750 K [8]. We have found this with our simulations using EAM potential at ≃ 900 K. Our expansion is about 5% at 1150 K, with respect to the bulk value. The disagreement between our MC result and MD result using the same potential resides not only in the difference of surface melting temperatures (ours is 700 K, the MD one is 1100 K) but also in the existence of the anomalous expansion.

For GP, we do not observe the anomalous expansion up to 11500 K as shown in Fig. 8. So the existence or not of an anomalous expansion depends on the potential at least with our MC results. Let us explain this as follows:

(i) With the EAM potential: The surface melts at a temperature much lower than that of the anomalous expansion. As we can see in Fig. 8, at around 700 K, the melting of the first layer causes the apparition of a little peak but ∆ 12 is still smaller than the bulk inter-layer distance. The fact that the anomalous expansion does not occur at the melting temperature of the first layer means that the melted surface is in a "two-dimensional" liquid state in the temperature region between 700 K and 900 K. This liquid surface layer is "detached" from the remaining crystal only at temperatures higher than 900 K. More accurate experiments of surface melting are desirable to allow a better understanding of this point.

(ii) With the GP: Since the surface melts at ≃ 1000 -1050 K (see Fig. 3), if there is an anomalous expansion, this should happen at a higher temperature. Given the fact that with EAM, our MC result indicates the anomalous expansion occurring 200 K above the surface melting temperature, we can imagine the same scenario with the GP: an expansion may be possible at ≃ 1200 -1235 K, namely at the bulk melting.

To close this section, let us summarize here that (i) the surface melting occurs below the bulk one, but the distance to the bulk melting temperature depends on the potential: the EAM and Gupta potentials give the same bulk melting temperature but different surface melting temperatures (ii) an anomalous dilatation of surface is possible only at a temperature much higher than the surface melting temperature.

IV. CONCLUDING REMARKS

Monte Carlo simulations with both EAM and Gupta potentials show two different surface melting temperatures for the (111) surface of a silver sample: while both EAM and Gupta potentials give the same bulk melting temperature (≃ 1200 K), the EAM potential yields a surface melting at ≃ 700 K and the Gupta potential shows a surface melting at ≃ 1000 -1050 K close to the bulk melting.

However, both potentials show a contraction of the topmost inter-layer distance of that surface at low temperatures, in good agreement with experiments and theories. Our results of the temperature-dependence of the inter-layer contraction indicate a strong potentialdependence: the variation of the contraction with increasing temperature shows a difference of the two potentials. Using the EAM potential, our MC simulations show an anomalous thermal expansion, namely the distance between the topmost layers is larger than that between two adjacent bulk layers. This can be explained by the fact that the anomalous thermal expansion occurs only at a temperature much higher than the surface melting temperature: the already-disordered surface layer is loosely attached to the crystal. This surface "dilatation" is in agreement with MEIS experiments and ab-initio DFT calculations. On the other hand, using the Gupta potential, our results show that the surface inter-layer distance varies with the temperature but it is always smaller than the bulk distance at least up to temperatures close to the bulk melting. We attribute the non observation of a surface dilatation in this case by the fact that surface melting occurs too close to the bulk melting. The high surface melting temperature has been found in MD simulations with EAM potential [START_REF] Ahlam | Rahman[END_REF]. The existence of surface dilatation has been observed in MEIS experiments [8] but not in LEED experiments [9] and in X-ray scattering [11]. We think that this difference comes in part from the difference of surface flatness, surface cleanness (contaminated or not) and experimental conditions. We hope that experimentalists could resolve this important question on anomalous high-temperature behavior of Ag(111). Experiments to determine the melting temperature of the Ag (111) surface is also needed to choose a potential appropriate for Ag crystal.

Figure 1 .

 1 Figure 1. (Color online) Experimental data by X-rays diffraction measurements(green circles) from Ref. [24] and the least square linear fit (red line).

Figure 2 .

 2 Figure 2. (Color online) Side view of a silver sample, with periodic boundary conditions, after several millions of MC steps at 300 K. Only atoms of the three topmost layers are allowed to move. See text for explanation.

Figure 3 .

 3 Figure 3. (Color online) Mean displacement amplitude (in Å) between atoms and their perfect ground-state lattice positions computed with Gupta potential.

Figure 4 .

 4 Figure 4. (Color online) Structure factor S K (green triangles) and O6 order parameter (red circles) of the first layer versus temperature for the EAM potential.

Figure 5 .

 5 Figure 5. Pictures of the system for two different temperatures computed with the EAM potential: 550 K (left) and 650 K (right). Only three layers are allowed to relax.

Figure 6 .

 6 Figure 6. Pictures of the system for two different temperatures computed with the Gupta potential : 950 K (left) and 1050 K (right).Only three layers are allowed to relax.

Figure 7 .

 7 Figure 7. (Color online) z-position distribution at different temperatures calculated with EAM and Gupta potentials. The surface at its non relaxed position corresponds to z=0.

Figure 8 .

 8 Figure 8. (Color online) Inter-layer ∆12 computed with Gupta and EAM potentials. Dashed line represents the bulk inter-layer distance. See text for comments.

Table I .

 I Gupta parameters for silver

	Parameter	Value a
	A (eV)	0.09944
	p	10.12
	q	3.37
	r0( Å)	2.892
	En	2.52 b
	a Reference [16]
	b Fitted with experimental bulk melting temperature in [20]

Table II .

 II EAM parameters for silver

	Parameter	Value a	Parameter	Value
	re	2.891814	Fn 1	-0.221025
	fe	1.106232	Fn 2	0.541558
	ρe	15.539255	Fn 3	-0.967036
	α	7.944536	F0	-1.75
	β	4.237086	F1	0
	A	0.266074	F2	0.983967
	B	0.386272	F3	0.520904
	κ	0.425351	η	1.149461
	λ	0.850703	Fe	-1.751274
	Fn 0	-1.729619


Bo Xiang

Chaohui Wang

Jean-Francois Deux

Alain Rahmouni

Nikos Paragios

3D CARDIAC SEGMENTATION WITH POSE-INVARIANT HIGHER-ORDER MRFS

Keywords: Left ventricle segmentation, pose-invariant shape prior, higher-order MRF

This paper proposes a novel pose-invariant segmentation approach for left ventricle in 3D CT images. The proposed formulation is modular with respect to the image support (i.e. landmarks, edges and regional statistics). The prior is represented as a third-order Markov Random Field (MRF) where triplets of points result to a low-rank statistical prior while inheriting invariance to global transformations. The ventricle surface is determined through triangulation where image discontinuities can be easily evaluated and the Divergence theorem provides an exact calculation of regional statistics acting on the image or a derived feature space. Promising results using boosting along with the learned prior demonstrate the potential of our method.

INTRODUCTION

Left ventricle (LV) segmentation is a prerequisite for quantitative analysis on diagnosing cardiovascular diseases. It is a well studied problem in biomedical imaging and most of the existing methods combine anatomical prior knowledge with observed visual support towards automatic delineation of LV.

Early approaches have considered slice-by-slice segmentation [START_REF] Jolly | Automatic segmentation of the left ventricle in cardiac MR and CT images[END_REF] using active contours.The use of active shape models [START_REF] Cootes | Active shape models-their training and application[END_REF] as well as active appearance models was a step forward where geometric and photometric variations were modeled using linear sub-spaces either at the frame level or during the entire cardiac cycle. Level set methods [START_REF] Huang | Metamorphs: Deformable shape and appearance models[END_REF] were also considered for LV segmentation towards providing accurate integration of boundary and visual support. Non-parametric statistics was an alternative to linear subspaces towards capturing more complex variations. Graph-based methods [START_REF] Grady | Random walks for image segmentation[END_REF] have also emerged due to their computational efficiency combined with certain user interaction. Spatio-temporal models involving poly-parametric class of transformations were recently investigated. Last but not least, the biomechanical prior models [START_REF] Sermesant | An electromechanical model of the heart for image analysis and simulation[END_REF] modeled the complete cardiac behavior which was calibrated through the image observations. All the above methods are modeled in a common subspace where statistics are learned in order to describe the variability among individuals. Such a process introduces strong statistical bias and at the same time makes the applicability of the model problematic when referring to diseased subjects. Furthermore, the construction of prior models is problematic due to the lack of sufficient number of training examples given the dimensionality of the model. Last but not least, these methods account partially for visual support with some of them acting on the boundaries, some of them on the regional statistics and some on anatomical landmarks. The aim of this paper is to propose a unified formulation that can integrate all existing sources of information, while being pose-invariant and able to cope with linear and non-linear statistics.

To this end, we adopt a Point Distribution Model associated with higher-order MRF. Triangulation on the surface provides third-order terms where an exact estimation of the boundary support and regional statistics is feasible through the Divergence theorem [START_REF] Xiang | Tagged cardiac mr image segmentation using boundary & regional-support and graph-based deformable priors[END_REF]. Furthermore, third-order cliques are also considered to encode prior knowledge [START_REF] Wang | 3D knowledge-based segmentation using pose-invariant higher-order graphs[END_REF]. Incremental deformation of the model coupled with landmark constraints is considered within a dual-decomposition approach to recover the optimal solution into a new image. Our approach is most closely related to the work [START_REF] Xiang | Tagged cardiac mr image segmentation using boundary & regional-support and graph-based deformable priors[END_REF][START_REF] Wang | 3D knowledge-based segmentation using pose-invariant higher-order graphs[END_REF]. The problem solved in [START_REF] Wang | 3D knowledge-based segmentation using pose-invariant higher-order graphs[END_REF] is simply graph matching (corresponding to a certain extend to the landmark part of the proposed formulation), while the 2D applied method in [START_REF] Xiang | Tagged cardiac mr image segmentation using boundary & regional-support and graph-based deformable priors[END_REF] neither inherits global pose invariance nor encompasses the use of landmarks.

The reminder of the paper is organized as follows. Section 2 presents the pose-invariant model and its associated statistical shape prior. In Section 3 we discuss the integration of visual and prior terms in a sound theoretical framework. Experimental results are presented in Section 4 and discussion in Section 5 concludes the paper.

POSE-INVARIANT POINT DISTRIBUTION MODEL

An essential task in knowledge-based image segmentation is to determine a statistical shape model as a prior knowledge which is then combined with visual support in an inference process. We represent the shape as a Point Distribution Model X = {x 1 , • • • , x n } (Fig. 1), consisting of a set V = {1, • • • , n} of n control points (marked as blue dots in Fig. 1(a)) lying on the boundary. Let x i (i ∈ V) denote the 3D position of the control point i. Meanwhile a subset V l ⊂ V of points are considered as anatomical landmarks (marked as green crosses in Fig. 1(a)). Furthermore, we associate this model with a set of triangles

T = {t ijk = (x i , x j , x k )|i, j, k ∈ V, i = j, j = k, k = i}
connecting any three points i, j, k. A subset S ⊂ T of triangles are defined as the surface triangles (Fig. 1(b)) which compose a closed mesh (2-manifold).

Based on the shape representation, the next task is to build a statistic model from a training set in order to construct the prior. The existing global approaches, such as ASM, often suffer from their lack of flexibility, while they also require a significant number of training samples to capture statistics due to their high-dimensional representation. To overcome these limitations, we propose a similarity-invariant model based on local interactions.

Let us consider a triplet of points i, j and k, the geometrical constraint between their positions t ijk = (x i , x j , x k ) can be charactered in a pose-invariant manner using two inner angles {α = (x i x j x k ), β = (x j x k x i )} of the triangle t ijk , while the third angle is a linear combination of the two. The inner angle can be computed as follows:

α = arccos --→ x i x j • ---→ x j x k x i x j x j x k (1)
An important property of the angle measurement is that it leads to a shape model invariant with respect to the global pose (i.e. translation, rotation and scale). Given a training set X = {X 1 , • • • , X m } of m instances of the shape, each example consists of n control points which are located in correspondences. For ∀i ∈ V, we have the set

X i = {x 1 i , • • • , x m i }
representing the instances of the i th control point in the shape. Then using a statistical model, we learn from the training set the probability density distributions p ijk (α, β) of the angles which are determined by a triplet of control points i, j and k. Let C = {(i, j, k)|i, j, k ∈ V, i = j, j = k, k = i} denote the triplet set corresponding to the triangle set T . Thus ∀c ∈ C, a distribution p c (α, β) is learned to describe the point position dependencies in the triplet c. With the accumulation of all these local distribution densities, the prior manifold is constructed through all triplet interactions. In this manner, our shape model encodes global consistency as well as local variations. Without loss of generality, both simple and complex distribution models can be used to learn the statistics. We use a Gaussian distribution model in practice.

So far, we build our pose-invariant point distribution model M = (X , C, P). It consists of a set of control points X , a triplet set C, and the statistical shape prior P = {p c (α, β)|c ∈ C}. Based on this representation,

KNOWLEDGE-BASED SEGMENTATION

Cardiac computed tomography (CT) image segmentation is challenging due to its complex background and low contrast. We aim to segment the myocardium structure which is the muscle between endocardium and epicardium. However, the boundaries of myocardium are not clearly visible because: (1) endocardium detection can be disturbed by the papillary muscles in the blood pool; (2) epicardium is not distinct between the myocardium and the right ventricle. To deal with these problems, we combine the anatomical landmarks, regional and boundary support to achieve a good segmentation.

We formulate the segmentation problem as an energy minimization problem. Given an initial position of the shape model, we aim to search for the optimal control point displacements to best compromise between the visual support and shape-fitness in an observed image. To this end, we use a higher-order Markov Random Field (MRF) formulation, in which a node models a control point and a third-order clique models the interaction by a triplet.

Given a graph G = (V, C), V the node set and C denotes the three-order clique set. We also denote a node set V l ⊂ V corresponding to the anatomical landmarks and a clique set C s ⊂ C corresponding to the surface triangles. Let X i (i ∈ V) denote the latent variable (i.e. the 3D position) of node i, and U i denote the candidate space for the configuration x i . To this end, the segmentation problem is formulated as the inference of the configuration Θ = (x i ) i∈V of all the nodes over the candidate space U = i∈V U i :

Θ opt = arg min Θ∈U E(Θ) (2)
The MRF energy E(Θ) is defined as a combination of landmark energy φ, region and boundary energy ψ (1) and shape prior energy ψ (2) while the first is a singleton term and the other two are higher-order (third-order) terms:

E(Θ) = i∈V l φ i (x i ) + c∈Cs ψ (1) c (x c ) + c∈C ψ (2) c (x c ) (3)
The definitions of each energy are described in following.

Landmark Correspondences

In the proposed graphical-model framework, we consider a set of salient points/anatomical landmarks (e.g. an apical control point) among the vertices of the surface model. They can be localized by a detector based on machine learning techniques such as Randomized Forests [START_REF] Breiman | Random forests[END_REF]. A singleton term is then introduced to enforce the consistency between the detected position and the estimated one for each landmark i:

φ i (x i ) = C 0 • |x i -xi | 2 (4)
where xi denotes the detected position for landmark i.

Regional Statistics & Boundary Support

Region-driven and boundary-based terms are based on the hypothesis that the populations in the image can be separated by their statistical properties. To represent the statistics, voxel intensity as a low-level representation has the limitation on the capability to recognize different objects, for example when the intensities of objects are overlapped. In this context, we use a high-dimensional feature vector for each voxel instead of an intensity value. The feature vector consists of patches of intensities and Gabor features [START_REF] Bernardino | Fast iir isotropic 2-d complex gabor filters with boundary initialization[END_REF] in view of their abilities to capture texture features in different scales and orientations.

In our case, a 5 × 5 × 5 patch and Gabor features with 3 scales and 6 orientations are used.

Based on the feature images, we learn an Adaboost classifier [START_REF] Collins | Logistic regression, AdaBoost and Bregman distances[END_REF] for the object and background and we apply the classifier responses to obtain a likelihood image of the testing image (Fig. 2). Each voxel is represented by a signed score in [-1, 1], whose sign symbol indicates the pixel label (i.e. negative for object and positive for background) and absolute value indicates the possibility to be object or background (distinguished in blue and red in Fig. 2(b)). The likelihood image is used for both region and boundary related potentials.

With regard to the regional potential, the exact calculation of regional likelihood on the volume can be factorized by using Divergence theorem, which states that volume V integral can be transferred into surface S integral.

V f dv = S (F • n) ds = c∈Cs γ region c (x c ) (5)
where f = ∇ • F is the likelihood image, n is the outward pointing unit normal of the boundary surface S. For surface triplet c ∈ C s , the regional potential is defined as follows:

γ region c (x c ) = C 1 • sc (F • n) ds (6)
where s c denotes the surface area determined by triplet c. Furthermore, we define the boundary related energy as a surface integral on surface triplet c ∈ C s : where M (x) is a boundary map which is the output of an boundary detector. For simplicity, we measure it by the dot product of the unit facet normal and the likelihood gradient at the position x, i.e. M (x) = -n • ∇f . This weak boundary detector is sufficient since we combine other image cues and shape prior in our framework. To this end, the data related energy defined on surface triplet c ∈ C s is:

γ boundary c (x c ) = C 2 • sc M (x)ds (7)
ψ (1) c (x c ) = γ region c (x c ) + γ boundary c (x c ) (8)

Prior Knowledge Constraints

The prior term encodes the spatial constraints of the control points configuration with respect to the shape manifold. The energy of local interaction on triplet c ∈ C is defined as:

ψ (2) c (x c ) = -log p c (α, β) (9)
where p c (α, β) denotes the possibility density of angles in the triplet c, which is obtained in the statistical shape model.

EXPERIMENTAL VALIDATION

We validate the proposed method on a dataset of 20 3D CT cardiac volumes. The volumes from different subjects have an approximate mean size of 512 × 512 × 250 voxels, where the voxel size is about 0.35 × 0.35 × 0.5mm 3 . Manual segmentation on the dataset is available and is considered as ground truth using for both learning and validation procedures. We perform a leave-one-out cross-validation on the dataset. A dual-decomposition optimization framework [START_REF] Komodakis | MRF optimization via dual decomposition: Message-passing revisited[END_REF] is adopted to perform the Maximum-a-Posteriori (MAP) inference for the proposed higher-order MRF. The weight parameters C 0 , C 1 , C 2 in MRF energy potentials are once adjusted manually to balance the corresponding terms for a good performance and then are kept for all the testing volumes. For quantitative evaluation, we compute the Dice coefficient between the segmented volume and the ground truth, which is one of the most common similarity measures between two volumes. We compare the proposed method with a Random Walks algorithm [START_REF] Grady | Random walks for image segmentation[END_REF] and a baseline ASM [START_REF] Cootes | Active shape models-their training and application[END_REF] using the same likelihood image. In Fig. 3 we present the distributions of the Dice coefficients obtained by the three methods evaluating on the whole dataset. Noted that a higher Dice coefficient implies a more accurate segmentation performance, our method outperforms the related work. Visual examples of experimental results are presented in Fig. 4, showing that the proposed method deals well with the shape variations and poor image quality. Both quantitative and qualitative analyses demonstrate the potential of our method.

CONCLUSION

In this paper, we propose a novel segmentation method based on graph theory with its explicit representations and have tested it using cardiac images. To the best of our knowledge, this is the first formulation that is pose invariant and can integrate common visual supports of any nature. It is achieved through the use of higher-order MRF where thirdorder cliques encode pose-invariant prior knowledge as well as boundary and regional statistics.

The main limitation of the method is the inherited computational complexity that is proportional to the number of higher-order cliques. We are investigating a more compact manner to encode prior knowledge that requires the smallest possible number of third-order cliques without altering the ability to express the shape manifold. Combined segmentation with landmark selection is another promising direction since geometric methods are known to be robust to the initial conditions. Last but not least, the extension of the method to account for the temporal nature of the cardiac motion through higher-order spatio-temporal priors is under investigation.

  (a) Distribution of the control points (b) Triangulated mesh

Fig. 1 .

 1 Fig. 1. 3D myocardium Point Distribution Model.

Fig. 2 .

 2 Fig. 2. Image Likelihood.

Fig. 3 .

 3 Fig. 3. Boxplot of the Dice coefficient statistics.