Jean-Dominique Polack

Hugo Dujourdy

Baptiste Pialot

Thomas Toulemonde

Energetic Wave Equation for Diffuse Sound Fields

Keywords: , 43, 55, Dt

The diffusion equation for modelling diffuse sound fields was proposed some fifty years ago on heuristic principles as an extension to Sabine's diffuse field model, and still receives much attention. Recently, the authors developed the model one step further to provide the missing relations between sound intensity and sound energy. This introduces some extra terms that, in case of non-Sabine spaces (narrow or flat rooms), can be defined with the help of the boundary conditions in terms of absorption and scattering coefficients on the walls. Integrating the divergence of the stress-energy tensor across the shortest dimensions of the space leads to a propagation equation of the Telegraph type, which can be solved using finite difference time domain simulation. For the two-dimensional case (open-space), numerical results are compared to measurements in real spaces. The comparison makes it possible to evaluate the absorption and scattering coefficients by an adjustment procedure.

Introduction

Since the first empirical treatment of reverberation by W.C. Sabine in 1894 [START_REF] Sabine | Collected papers acoustics[END_REF], many authors have noted discrepancies between theory and measurements. But few have questioned the very foundation of Sabine's model, the so-called diffuse field approximation. Among those, ISO 14257 [START_REF]Acoustics -Measurement and parametric description of spatial sound distribution curves in workrooms for evaluation of their acoustical performance[END_REF] recommends measuring the decay of the sound field with distance, not reverberation time, to characterize industrial halls Ollendorff [START_REF] Ollendorff | Statistical room acoustics as a problem of diffusion: a proposal[END_REF] was the first to come up with the idea of a diffusive model for non-uniform distribution of sound energy in disproportionate rooms. His model escaped attention before being rediscovered independently 20 years later by Picaut et al. [START_REF] Picaut | A mathematical model of diffuse sound field based on a diffusion equation[END_REF].

Here we present a new approach based on the stressenergy tensor. Like the diffusion equation [START_REF] Morse | Methods of Theoretical Physics[END_REF], this approach was initiated by P.M. Morse [START_REF] Morse | Theoretical Acoustics[END_REF]. It allows us to write the conservation of sound energy in the form of a Telegraph equation, which is solved by a finite difference method. We then compare the resulting simu-(c) European Acoustics Association lation to actual measurements in an open-space office. This approach has constituted the doctoral research project of the second author [START_REF] Dujourdy | Diffusion acoustique sur les lieux de travail[END_REF], and the topic of two publications [START_REF] Dujourdy | An energetic wave equation for modelling diffuse sound fields -application to corridors[END_REF][START_REF] Dujourdy | An energetic wave equation for modelling diffuse sound fields -application to open offices[END_REF].

State of the art

Sabine's model

Sabine's model considers the conservation of sound energy at any position in space:

∂ t w + ∇ • I = 0 (1)
where ∂ t is the partial derivation relative to time, ∇• is the divergence operator, w the total instantaneous density of acoustic energy, and I the acoustic intensity vector. It adds a heuristic relationship between intensity and sound energy, based on the diffuse field hypothesis:

• sound energy is equi-distributed in space;

• it propagates isotropically in space.

The relationship that binds energy and intensity is then:

| I| = w 4π c ( 2
)
where c is the speed of and the total energy incident on the walls with total surface S is given by:

S I • ndS = w 4 cS (3)
Taking into account the fact that the acoustic energy is constant in the volume, the integration of the equation over the whole volume V gives the well-known equation:

V ∂ t w + ᾱSc 4 w = 0 ( 4
)
whose solution is an exponential function of time with decrement ᾱSc 4V , where ᾱ is the average absorption coefficient on the walls.

Ollendorff-Picaut's model

Following a suggestion by Morse and Feshback [START_REF] Morse | Methods of Theoretical Physics[END_REF], Ollendorf and Picaut proposed another heuristic relationship between acoustic intensity and energy:

I = -δ ∇w (5)
where the total energy is no longer equi-distributed in space. Here, ∇ is the gradient operator and δ is the diffusion coefficient. They then get an equation of diffusion:

∂ t w -δ∆w = 0 (6)
or by "renormalizing" the absorption within the volume:

∂ t w + σw -δ∆w = 0 (7)
This equation accepts an analytical solution in infinite space:

w( r, t) = 1 (4πδt) d 2 exp{- r 2 4δt -σt} (8)
where d is the space dimension. In flat spaces (d = 2) or narrow spaces (d = 1), it can be shown that σ = ᾱSc 4V , and solution 8 generalizes Sabine's exponential decay.

The stress-energy tensor

General equation

Following the example of Morse and Ingard [START_REF] Morse | Theoretical Acoustics[END_REF], Dujourdy et al. [START_REF] Dujourdy | An energetic wave equation for modelling diffuse sound fields -application to corridors[END_REF][START_REF] Dujourdy | An energetic wave equation for modelling diffuse sound fields -application to open offices[END_REF] propose to introduce the symmetric stress-energy tensor:

T =     E E tx E ty E tz E tx E xx E xy E xz E ty E xy E yy E yz E tz E xz E yz E zz     = E J T J E (9)
where J = I c and J T is the transposed matrix of J; E is the wave-stress tensor, with diagonal elements satisfying relation E = E xx + E yy + E zz since kinetic and potential energy are equal on average. The conservation of the tensor T is simply expressed by a null covariant derivative:

∇ 4 • T = 0 ( 10
)
where

∇ 4 = ( 1 c ∂ t , ∇)
is the quadri-divergence operator.

Dimensional reduction

The dimensional reduction of the system is achieved by integrating on the small dimensions of the volume. Thus, for a flat space such as an open-space office, one integrates on the height of the space (coordinated z), by introducing the average value of the quantities, except those for which the derivation is carried out according to z and which give the flows entering the walls. For the conservation of energy (first line of the stress-energy tensor), we obtain:

1 c ∂ t El z + ∂ x J x l z + ∂ y J y l z + J + z -J - z = 0 (11)
where l z is the height of the space, and where J + z , resp. -J - z , is the energy flow through the ceiling (side +), resp. through the floor (side -), i.e. the sound intensity absorbed by each of them. Similarly, for the following two lines of the stressenergy tensor, the following conservation relationships are obtained:

1 c ∂ t J x l z + ∂ x E xx l z + ∂ y E xy l z (12) + E + xz -E - xz = 0 1 c ∂ t J y l z + ∂ x E yx l z + ∂ y E yy l z (13)
+ E + yz -E - yz = 0
where E + xz and E + yz , resp. -E - xz and -E - yz , are the wave stresses that are applied on the ceiling (side +), resp. on the floor (side -). It should be noted that the last line of the stress-energy tensor does not provide additional information. The dimensional reduction is also accompanied by a simplification of the residual wave-stress tensor E, which remains symmetrical. We postulate:

E xx = E yy = E 2 (14)
E xy = 0 (15)
that is, an isotropic distribution of energy in the twodimensional space.

Boundary conditions

There remains to write the energy balance on the ceiling and the floor of the space, and more generally on the walls. To do this, we introduce the energy distribution function f ( r, v, t), which gives the probability that a wave is localized at position r heading in the direction of velocity v at instant t. Following the example of [START_REF] Morse | Methods of Theoretical Physics[END_REF][START_REF] Picaut | A mathematical model of diffuse sound field based on a diffusion equation[END_REF], we define f by:

f ( r, v, t) = E 4π + 3 J 4π • v c (16)
The energy incident on a wall is obtained by integrating on the solid angle

Ω = (θ, φ) with θ ∈ [0, π 2 ], φ ∈ [0, 2π],
and the energy reflected by integration on the solid angle

Ω = (θ, φ) with θ ∈ [ π 2 , π], φ ∈ [0, 2π].
All calculations done, we get:

J z,inc = E 4 + Jz 2 J z,ref = E 4 -Jz 2 (17)
and we recover Jing and Xiang's definition of the absorption coefficient A [10]:

J z,inc -J z,ref = J z,abs = αJ z,inc (18)
= α 2(2 -α) E = A 4 E
where α is a proportionality factor between absorbed and incident intensities, equivalent to an absorption coefficient.

For the stress balance on the walls, it is necessary to add terms to the energy distribution function:

f ( r, v, t) = E 4π + 3 J 4π • v c (19)
+ 15 4π {E xz cos θ sin θ cos φ + E yz cos θ sin θ sin φ} For the second line of the stress-energy tensor, the incoming and outgoing stresses are obtained by integration in the form of:

M xz,in = 3 16 J x + 1 2 E xz (20) M xz,out = 3 16 J x - 1 2 E xz (21)
This allows the introduction of a scattering coefficient D:

M xz,in -M xz,out = M xz,scat = βM xz,in (22) = 3 4 β 2(2 -β) J x = D 4 J x
where β is a proportionality factor between scattered and incoming stresses. Similarly, the third line of the stress-energy tensor becomes:

M yz,in -M yz,out = M yz,scat (23)
= 3 4 β 2(2 -β) J y = D 4 J y

Solution

The substitution of equations 19, 23, and 24 in equations 11-14, as well as the use of the mean free path (λ = 2l z in a 2-dimensional space), leads to the equation system:

1 c ∂ t E + ∂ x J x + ∂ y J y + A λ E = 0 1 c ∂ t J x + ∂ x E xx + ∂ y E xy + D λ J x = 0 (24) 1 c ∂ t J y + ∂ x E yx + ∂ y E yy + D λ J y = 0
which reduces to the Telegraph equation:

1 c 2 ∂ tt E -∆ E 2 + A + D λc ∂ t E + AD λ 2 E = 0 (25)

Finite difference time domain simulation

We choose to solve numerically equation 25 by a numerical Finite Difference Time Domain (FDTD) method. For this, we choose a non staggered grid and a centred-time centred-space scheme. The derivatives are approximated by:

∂ 2 E ∂t 2 t x,y = E n+1 i,j -2E n i,j + E n-1 i,j ∆t 2 + O(∆t) 2 ∂ 2 E ∂x 2 t x,y = E n i+1,j -2E n i,j + E n i-1,j ∆x 2 + O(∆x) 2 ∂ 2 E ∂y 2 t x,y = E n i,j+1 -2E n i,j + E n i,j-1 ∆y 2 + O(∆y) 2 ∂E ∂t t x,y = E n+1 i,j -E n-1 i,j 2∆t + O(∆t) 2 (26) ∂E ∂x t x,y = E n i+1,j -E n i-1,j 2∆x + O(∆x) ∂E ∂y t x,y = E n i,j+1 -E n i,j-1 2∆y + O(∆y)
where O(•) is the truncation error and the ∆x, ∆y and ∆t the discretization steps in space and time. The discrete equation then takes the form:

E n+1 i,j (a + 1) = E n-1 i,j (a -1) + E n i,j (2 -C 2 rx -C 2 ry -b) + 1 2 C 2 rx (E n i+1,j + E n i-1,j ) (27)
+ 1 2 C 2 ry (E n i,j+1 + E n i,j-1 ) + O[(∆t) 2 , (∆x) 2 , (∆y) 2 ]
completed by 4 boundary equations on the 4 lateral walls of the space [START_REF] Dujourdy | An energetic wave equation for modelling diffuse sound fields -application to open offices[END_REF]. C rx = c∆t ∆x and C ry = c∆t ∆y are the Courant-Friedrichs-Lewy coefficients for each dimension, and a = (A + D) c∆t 2λ , b = AD( c∆t λ ) 2 . The scheme is explicit, and the accuracy is second order in time and space.

Stability

We study the stability of the scheme with the Von Neumann procedure, which looks for a solution of equation 28 in the form E n i,j = Z n e j(θi+ωj) . All calculations made, we obtain two conditions of stability [START_REF] Dujourdy | An energetic wave equation for modelling diffuse sound fields -application to open offices[END_REF]:

C 2 rx (1 + AD( ∆x 2λ ) 2 ) < 1 (28) C 2 ry (1 + AD( ∆y 2λ ) 2 ) < 1 (29)

Results

The discrete equation 28 was applied to a quasiinfinite two-dimensional space of 200m by 200m. The source is located in the centre of the space. The absorption and scattering coefficients range from α = 0.01 to α = 0.8 and β = 0.01 to β = 1.5. The mean free path is λ = 5.4m, and the discretization steps are set at ∆x = ∆y = 1m and ∆t = 1.10 -3 s. The speed of sound is c = 344m/s.

Spatial decay

Figures 1-2 show spatial sound decays, and compare them to the analytical solution in the steady-state case:

E = E 0 exp(- √ 2AD λ r) √ r (30)
for large values of r, where r is the distance to the source and E 0 an integration constant.

It is clear that the numerical solution corresponds to the analytical solution, except for the very low values of the coefficients, or near the walls where reflections on the latter exert influence.

Temporal decay

Figures 34show temporal sound decays at 15m from the source, and compare them to the two analytical approximations corresponding resp. to a large value of one of the coefficients and to equal coefficients. In the first case, equation 25 is reduced to an equation of diffusion, with solution:

E = E 1 e (-(A+D)r 2 2λct -ADc (A+D)λ t) t (31)
In the second case (A = D), the solution is a travelling wave packet [START_REF] Bruneau | Introduction aux théories de l'acoustique[END_REF]:

E + = E 2 e (-A λ ct) U ( c √ 2 t -r) 2π ( c √ 2 t) 2 -r 2 (32)
where U (•) is the Heaviside function. It is clear that the numerical solution corresponds to the analytical solutions with appropriate choice of the integration constants E 1 and E 2 , especially for the large values of t when the algorithm has stabilized. In Figures 34, E 1 = 74dB for analytical solution 31; for analytical solution 32, E 2 = 101dB, but reduces to 95dB for large values of α and β. Similar results are obtained at the other distances.

Measurements

We were able to conduct measurements in three openspace offices in Paris. Here we present the results obtained for Ovalie (Figure 5), an office space of 75m long by 35m wide, oval-shaped with a central islet. Its surface is 1260m 2 , and its height 2.7m, yielding a mean free path of 5.4m. The office was furnished with desks, chairs and computers during measurements, but without the presence of users. Figure 6 shows the spatial decay in Ovalie along the green axis and Figure 7 the temporal decay at 10m from the S3 source. Both Figures compare measured decays to the decays calculated after adjusting the parameters α and β. The measurement frequency is the 1kHz octave, and the optimum values of α and β are respectively 0.21 and 0.7. The same optimum values are obtained at the other distances. In Figure 7, the levels of the temporal decay are slightly overestimated, an observation that is general to all analysed positions.

Adjustment procedure

The absorption and scattering coefficients are obtained by comparison with the measurements. Now we see in Figures 34that the smaller of the α or β values drives the temporal decay. In other words, temporal decay allows the evaluation of the smaller of the coefficients, in general α, and the second coefficient is obtained from spatial decay. Thus an iterative procedure was used to estimate α and β, assuming first β large for a first estimation of the absorption coefficient by equation 31, which reduces to:

E ≈ E 1 e (-(D)r 2 2λct -Ac λ t) t (33)
Then β is inferred from the spatial decay by equation 30. The procedure however exaggerates the slopes of the two decays, so that the coefficients are slightly adjusted to the margin for better agreements with the measures.

Conclusions

We have shown that energy conservation can be broken down into two equations: the conservation of total energy; and the conservation of the acoustic intensity. The two equations combine into a single tensor equation, the conservation of the stress-energy tensor. The resolution of these equations involves two coefficients on the walls: an absorption coefficient and a scattering coefficient. The smaller of these two coefficients, usually the absorption coefficient, pilots the temporal decay, and the second coefficient is then obtained from spatial decay. This adjustment procedure has been validated in real cases of open-space offices.

Figure 1 .

 1 Figure 1. Spatial decays for different α with β = 0.2; lozenges correspond to steady-state analytical solution, with E0 = 74dB.

Figure 2 .Figure 3 .Figure 4 .

 234 Figure 2. Spatial decays for different β with α = 0.2; lozenges correspond to steady-state analytical solution, with E0 = 74dB.

5 .Figure 6 .

 56 Figure 6. Model adjustment by comparison with spatial decay measured at Ovalie: source S3, green axis, α = 0.21 and β = 0.7.

Figure 7 .

 7 Figure 7. Model adjustment by comparison with temporal decay at 10m from the source measured at Ovalie: source S3, green axis, α = 0.21 and β = 0.7.

Acknowledgement

This work was carried out within the framework of a CIFRE convention (nr 2012/1184) between Impedance S.A.S. and Université Pierre et Marie Curie. The authors thank Impedance for providing access to the open-space offices for the measurements.


Damien Lamberton
email: damien.lamberton@u-pem.fr

On the binomial approximation of the American put

We consider the binomial approximation of the American put price in the Black-Scholes model (with continuous dividend yield). Our main result is that the error of approximation is O((ln n) α /n), where n is the number of time periods and the exponent α is a positive number, the value of which may differ according to the respective levels of the interest rate and the dividend yield.

The binomial approximation

Consider the Black-Scholes model, in which the stock price at time t is given by

where, under the risk-neutral probability measure, (B t ) t≥0 is a standard Brownian motion. Here, r is the instantaneous interest rate, and d is the dividend rate (or the foreign interest rate in the case of forex options). We assume r > 0 and d ≥ 0.

Denote by P the price function of the American put with maturity T and strike price K, so that

P (t, x) = sup τ ∈T 0,T -t E x e -rτ f (S τ ) , 0 ≤ t ≤ T, x ∈ [0, +∞),
with f (x) = (K -x) + , and E x = E (• | S 0 = x). Here, T 0,t denotes the set of all stopping times with respect to the Brownian filtration, with values in the interval [0, t].

For technical reasons (especially for the derivation of regularity estimates for the second time derivative of the price function), it is more convenient to use the log-stock price. So, we introduce

X x t = x + µt + σB t , with µ = r -d - σ 2 2 ,
and U(T, x) = sup τ ∈T 0,T E e -rτ ϕ(X x τ ) , with ϕ(x) = (K -e x ) + . We then have P (t, x) = U(T -t, ln(x)), t > 0, x > 0.

Note that U(t, x) satisfies the following parabolic variational inequality max -∂U ∂t + (A -r)U, ϕ -U = 0, with the initial condition U(0, .) = ϕ.

Here, A is the infinitesimal generator of X, namely

A = σ 2 2 ∂ 2 ∂x 2 + µ ∂ ∂x .
Recall that, for each T > 0, there is a real number b(T ) ≤ ln(K) such that U(T, x) > ϕ(x) ⇔ x > b(T ).

In fact, if (b(t), 0 ≤ t ≤ T ) is the exercise boundary of the American put with maturity T , we have b(t) = ln(b(T -t)). We will also need the European value function, defined by Ū (T, x) = E e -rT ϕ(X x T ) .

Note that Ū (0, .) = ϕ and

- ∂ Ū ∂t + (A -r) Ū = 0.
Note that, in Section 3, the function Ū will be denoted by u ϕ . We now introduce the random walk approximation of Brownian motion. To be more precise, assume (X n ) n≥1 is a sequence of i.i.d. real random variables satisfying EX 2 n = 1 and EX n = 0, and define, for any positive integer n, the process B (n) by

B (n) t = T /n [nt/T ] k=1 X k , 0 ≤ t ≤ T,
where [nt/T ] denotes the greatest integer in nt/T .

We will assume the following about the common distribution of the X n 's (cf. hypothesis (H4) of [START_REF] Lamberton | Brownian optimal stopping and random walks[END_REF]). Note that, in the binomial case, X 1 takes its values in {-1, +1}.

(H4) The random variable X 1 is bounded and satisfies EX 2 1 = 1 and EX 1 = EX 3 1 = 0.

In the following, we fix S 0 and set P 0 = P (0, S 0 ) = U(T, ln S 0 ).

Note that, if we introduce the notation g(x) = (K -S 0 e σx ) + , we have

P 0 = sup τ ∈T 0,T E e -rτ g(µ 0 τ + B τ ) ,
with µ 0 = µ/σ. We now have a natural approximation of P 0 , given by

P (n) 0 = sup τ ∈T (n) 0,T E e -rτ g(µ 0 τ + B (n) τ ) ,
where T

(n) 0,T denotes the set of all stopping times (with respect to the natural filtration of B (n) ), with values in [0, T ] ∩ {0, T /n, 2T /n, . . . , (n -1)T /n, T }. Our main result is the following.

Theorem 1.1. There exists a positive constant C such that, for all positive integers n,

-C (ln n) ᾱ n ≤ P (n) 0 -P 0 ≤ C (ln n) α n ,
where ᾱ = α = 1 if d > r, and ᾱ = 3/2, α = 5/4 if d ≤ r.

The above estimates improve our previous results (see [START_REF] Lamberton | Brownian optimal stopping and random walks[END_REF], Theorem 5.6) which gave an upper bound of the form C √ ln n n 4/5

. Note that, for European options, the error estimate is O(1/n) (see [START_REF] Diener | Asymptotics of the price oscillations of vanilla option in a tree model[END_REF], [START_REF] Walsh | The Rate of Convergence of the Binomial Tree Scheme[END_REF]). We also mention the results of [START_REF] Liang | Optimal Convergence Rate of the Binomial Tree Scheme forAmerican Options with Jump Diffusion and Their Free Boundaries[END_REF] about finite difference schemes, which give the rate O(1/ √ n), but their estimate is uniform over the time interval, while we concentrate on the error estimate for a fixed time. The paper [START_REF] Liang | Optimal Convergence Rate of the Binomial Tree Scheme forAmerican Options with Jump Diffusion and Their Free Boundaries[END_REF] also has results about the approximation of the exercise boundary. We also refer to [START_REF] Silvestrov | American Type Options. Stochastic Approximation Methods[END_REF] and its references for a review of recent results on the approximation of American option prices.

Our approach remains the same as in [START_REF] Lamberton | Brownian optimal stopping and random walks[END_REF]: we relate the error estimates to the regularity of the value function. The improvement comes from a refinement of the quadratic estimates for the second order time derivative, in the spirit of Friedman and Kinderlehrer (see [START_REF] Friedman | Parabolic variational inequalities in one space dimension and smoothness of the free boundary[END_REF] and [START_REF] Kinderlehrer | An introduction to variational inequalities and their applications[END_REF]). We also exploit the smoothness of the exercise boundary and its asymptotic properties close to maturity.

The constant C in Theorem 1.1 is related to the Berry-Esseen estimate and to the regularity of the value function. Although it is hard to keep track of the constants in the regularity estimates, it may be worth mentioning that they remain uniform with respect to µ and σ as long as (µ, σ) remains in a compact subset of R × (0, ∞). A consequence of this observation is that the bounds in Theorem 1.1 are also valid for variants of the approximation in which the process approximating ln(S t /S 0 ), instead of being µt + σB

(n) t , is given by µ n t+ σ n B (n) t at discrete times t, with µ n = µ + O(1/n) and σ 2 n = σ 2 + O(1/n),
as occurs in the classical risk-neutral approximation. Indeed, standard arguments show that the value function is locally Lipschitz-continuous with respect to σ 2 (away from 0) and µ.

The paper is organized as follows. In the next Section we recall some results of [START_REF] Lamberton | Brownian optimal stopping and random walks[END_REF]. Section 3 is devoted to estimates for the derivatives of the value function. The estimates are then used in Sections 4 and 5 to prove Theorem 1.1: in Section 4, we give an upper bound for P (n) 0 -P 0 and in Section 5, we derive the lower bound.

Acknowledgement: The research on this paper has been stimulated by fruitful discussions on the approximation of American options with Martijn Pistorius, to whom the author is very grateful.

The value function and the approximating process

As in [START_REF] Lamberton | Brownian optimal stopping and random walks[END_REF], we introduce the modified value function

u(t, x) = e -rt U(T -t, ln(S 0 ) + µt + σx), t ≥ 0, x ∈ R.
We have P 0 = u(0, 0) and u(T, x) = e -rT U(0, ln(S 0 ) + µT + σx) = e -rT (K -S 0 e µT +σx ) + and, for t ∈ [0, T ], u(t, x) ≥ e -rt (K -S 0 e µt+σx ) + = e -rt g(µ 0 t + x).

(1)

We will need the European analogue of u, namely

ū(t, x) = e -rt Ū(T -t, ln(S 0 ) + µt + σx) = e -rT E (g(µ 0 T + x + B T -t )) , t ≥ 0, x ∈ R.
We will also use the notation:

h = T n .
With this notation, we have

B (n) t = √ h [t/h] k=1 X k , 0 ≤ t ≤ T.
We have, for all t ∈ {0, h, 2h, . . . , (n -1)h, nh = T } (cf. Proposition 3.1 of [START_REF] Lamberton | Brownian optimal stopping and random walks[END_REF]),

u(t, B (n) t ) = u(0, 0) + M t + t/h j=1 Du((j -1)h, B (n) (j-1)h ),
where (M t ) 0≤t≤T is a martingale (with respect to the natural filtration of B (n) ), such that M 0 = 0, and

Du(t, x) = E u t + h, x + √ hX 1 -u(t, x), 0 ≤ t ≤ T -h, x ∈ R.

The above decomposition of u(t, B

(n) t ) (which is in fact Doob's decomposition) can be viewed as a discrete version of Itô's formula, which, for a smooth function v :

[0, T ]×R → R, implies that v(t, B t ) -t 0 δv(s, B s )ds is a (local) martingale, where δv = ∂v ∂t + 1 2 ∂ 2 v ∂x 2 .
It is also easy to check that, if v is smooth and

Dv(t, x) = E v t + h, x + √ hX 1 -v(t, x), we have (1/h) × Dv(t, x) = δv(t, x) + O(h).
The main technical difficulty that we have to deal with is the lack of smoothness of the modfied value function u.

Remark 2.1. The derivatives of u are related to those of U by the following formulas. We have

∂u ∂t (t, x) = e -rt - ∂U ∂t + µ ∂U ∂x -rU (T -t, ln(S 0 ) + µt + σx)
and

∂ 2 u ∂t 2 (t, x) = e -rt ∂ 2 U ∂t 2 -2µ ∂ 2 U ∂t∂x + µ 2 ∂ 2 U ∂x 2 +2r ∂U ∂t -2rµ ∂U ∂x + r 2 U (T -t, ln(S 0 ) + µt + σx).
We also have

δu(t, x) = ∂u ∂t (t, x) + 1 2 ∂ 2 u ∂x 2 (t, x) = e -rt - ∂U ∂t + (A -r)U (T -t, ln(S 0 ) + µt + σx) = e -rt (A -r)ϕ(ln(S 0 ) + µt + σx)1 {ln(S 0 )+µt+σx≤ b(T -t)} ,
where the last equality follows from regularity results (see, for instance, [START_REF] Jaillet | Variational inequalities and the pricing of American options[END_REF]).

We will need a more precise description of the operator D, given by the following proposition (see Proposition 3.4 of [START_REF] Lamberton | Brownian optimal stopping and random walks[END_REF]). For convenience, we denote by X a random variable with the same distribution as X 1 , which is independent of the sequence

(X n ) n≥1 . Proposition 2.1. Assume that (H4) is satisfied and that v is a function of class C 3 on [0, T ] × R. For 0 ≤ t ≤ T -h and x ∈ R, define Dv(t, x) = 2 √ h 0 dξ ξ 0 dzE X ξ -X 2 (ξ -z) ∂ 2 v ∂t∂x (t + ξ 2 , x + zX) .
We have

Dv(t, x) = Dv(t, x) + 2 √ h 0 dξ ξ 0 dzE X 2 δv(t + ξ 2 , x + zX) ,
with the notation δv = ∂v ∂t + 1 2

∂ 2 v ∂x 2 ,

and

Dv(t, x) = 2 √ h 0 dξ ξ 0 dz(ξ -z)E X 2 ξ -X 2 (ξ -z) 2 ∂ 3 v ∂t∂x 2 (t + ξ 2 , x + zX) . Remark 2.2. Note that, if δv(s, x + zX) = 0 for all s ∈ [t, t + h] and z ∈ [0, √ h], we have v(t, x) = Dv(t, x).
From the last equality in Proposition 2.1, we derive the following estimates.

Dv(t, x) ≤ 2 √ h 0 ξ 2 dξ ξ 0 dzE X 2 + X 4 2
∂ 3 v ∂t∂x 2 (t + ξ 2 , x + zX) ≤ √ h √ h 0 2ξdξE ξ 0 dz X 2 + X 4 2 ∂ 3 v ∂t∂x 2 (t + ξ 2 , x + zX) ≤ √ h t+h t dsE dy1 {|y-x|≤ √ h|X|} |X| + |X| 3 2 ∂ 3 v ∂t∂x 2 (s, y) = √ h t+h t ds dyE 1 {|y-x|≤ √ h|X|} |X| + |X| 3 2 ∂ 3 v ∂t∂x 2 (s, y)
We know from Proposition 3.2 of [START_REF] Lamberton | Brownian optimal stopping and random walks[END_REF] (based on Berry-Esseen estimates) that, for every k ∈ (1, 3], there exists a positive constant C k (which does not depend on X), such that, for all y ∈ R, n ≥ 1 and j ∈ {1, 2, . . . , n},

E   |X| + |X| 3 2 1 B (n) jh -y ≤ √ h|X|   ≤ C k √ j E (|X| 3 ) 1 + E |X| 3+k 1 + |y| k .
Hence, for j = 1, . . . , n -1,

E Dv(jh, B (n) jh ) ≤ C k,X √ j √ h jh+h jh ds dy 1 + |y| k ∂ 3 v ∂t∂x 2 (s, y) ≤ C k,X h √ 2 jh+h jh ds √ s dy 1 + |y| k ∂ 3 v ∂t∂x 2 (s, y) , ( 2
)
where, for the last inequality, we used the inequality jh ≥ (j + 1)h/2.

3 Estimates for the second order time derivative

In this section, we refine the regularity results that we used in [START_REF] Lamberton | Brownian optimal stopping and random walks[END_REF]. We first establish some elementary L 1 -estimates. Then, we obtain a quadratic estimate for the second order time derivative of the difference Ũ = U -Ū . For the definition of the relevant weighted Sobolev spaces, we will use the notation

ν j (dx) = dx (1 + x 2 ) j/2 , j > 1.
3.1 Some elementary L 1 -estimates Proposition 3.1. Assume that the function ϕ is continuous and satisfies ϕ ∈ L 1 (ν j ), ϕ ′ ∈ L 1 (ν j ) and the second derivative ϕ ′′ is a Radon measure on R, with R

|ϕ ′′ (dz)| (1+z 2 ) j/2 < ∞. Let u ϕ (t, x) = e -rt E(ϕ(X x t )), t ≥ 0, x ∈ R. Then, for all T > 0, there exists a constant C T > 0, such that ∀t ∈ (0, T ], ∂ 2 u ϕ ∂t 2 (t, .) L 1 (ν j ) ≤ C T t .
We will easily deduce this proposition from the following lemma.

Lemma 3.1. If ρ is a Radon measure on R and q a nonnegative integrable function on R, we have

||ρ * q|| L 1 (ν j ) ≤ 2 j/2 R |ρ(dz)| (1 + z 2 ) j/2 ∞ -∞ q(x)(1 + x 2 ) j/2 dx.
We also have, for any measurable function f on R,

∀y ∈ R, ||f (. -y)|| L 1 (ν j ) ≤ 2 j/2 (1 + y 2 ) j/2 ||f || L 1 (ν j ) . Proof: We have ||ρ * q|| L 1 (ν j ) ≤ ∞ -∞ dx (1 + x 2 ) j/2 R |ρ(dz)|q(x -z) = R |ρ(dz)| (1 + z 2 ) j/2 ∞ -∞ q(x -z) (1 + z 2 ) j/2 (1 + x 2 ) j/2 dx = R |ρ(dz)| (1 + z 2 ) j/2 ∞ -∞ q(x) (1 + z 2 ) j/2 (1 + (x + z) 2 ) j/2 dx. Note that z 2 ≤ 2((x + z) 2 + x 2 ),
so that we deduce

1 + z 2 1 + (x + z) 2 ≤ 1 + 2(x + z) 2 + 2x 2 1 + (x + z) 2 ≤ 2(1 + x 2 ). Hence ||ρ * q|| L 1 (ν j ) ≤ 2 j/2 R |ρ(dz)| (1 + z 2 ) j/2 ∞ -∞ q(x)(1 + x 2 ) j/2 dx.
Similarly, we have, for any measurable function f and y ∈ R,

||f (. -y)|| L 1 (ν j ) = |f (x -y)| dx (1 + x 2 ) j/2 = |f (x)| dx (1 + (x + y) 2 ) j/2 = |f (x)| 1 + x 2 1 + (x + y) 2 j/2 dx (1 + x 2 ) j/2 ≤ |f (x)| 1 + 2(x + y) 2 + 2y 2 1 + (x + y) 2 j/2 dx (1 + x 2 ) j/2 ≤ 2 j/2 (1 + y 2 ) j/2 ||f || L 1 (ν j ) . ⋄ Proof of Proposition 3.1: We have u ϕ (t, x) = e -rt ∞ -∞ ϕ(x + y) exp - (y -µt) 2 2σ 2 t dy σ √ 2πt = e -rt p t * ϕ(x), with p t (x) = 1 σ √ 2πt exp - (x + µt) 2 2σ 2 t = 1 σ √ t n x + µt σ √ t .
Here, n denotes the standard normal density function.

On the other hand, we know that u ϕ satisfies the equation

∂u ϕ ∂t = (A -r)u ϕ , (3)
so that

∂u ϕ ∂t (t, .) = e -rt (A -r)p t * ϕ = e -rt p t * [(A -r)ϕ].
It follows from our assumptions that

(A -r)ϕ is a Radon measure satisfying R |(A -r)ϕ(dz)| 1 (1 + z 2 ) j/2 < ∞.
So that, using Lemma 3.1,

∂u ϕ ∂t (t, .) L 1 (ν j ) ≤ C j ∞ -∞ p t (x)(1 + |x| j )dx = C j ∞ -∞ 1 σ √ t n x + µt σ √ t (1 + |x| j )dx = C j ∞ -∞ n (y) (1 + |yσ √ t -µt| j )dy ≤ C j 1 + t j .
On the other hand, by differentiating (3), we have

∂ 2 u ϕ ∂t 2 = (A -r) ∂u ϕ ∂t = e -rt ((A -r)p t ) * (A -r)ϕ.
Hence, using Lemma 3.1, and the definition of p t ,

∂ 2 u ϕ ∂t 2 (t, .) L 1 (ν j ) ≤ C j ∞ -∞ |(A -r)p t (x)|(1 + |x| j )dx ≤ C j t 1 + t j . ⋄

Quadratic estimates

Recall the notation:

U(t, x) = sup τ ∈T 0,t E e -rτ ϕ(X x τ ) , u ϕ (t, x) = e -rt E(ϕ(X x t )), t ≥ 0, x ∈ R,
with ϕ(x) = (K -e x ) + . We now introduce the difference Ũ = U -u ϕ (which corresponds to the early exercise premium). We have the following L 2 -estimate for the second time derivative of Ũ = U -u ϕ .

Theorem 3.1. Fix T > 0 and j > 1. There exists a constant C > 0 such that, for all ξ ∈ (0, T ],

T ξ (t -ξ) ∂ 2 Ũ ∂t 2 (t, .) 2 L 2 (ν j ) dt ≤ C 1 + | ln ξ| β , with β =      3/2, if d ≤ r, 1, if d > r.
This estimate is closely related to Theorem 2.4 of [START_REF] Lamberton | Brownian optimal stopping and random walks[END_REF], a variant of results due to Friedman and Kinderlehrer (see [START_REF] Friedman | Parabolic variational inequalities in one space dimension and smoothness of the free boundary[END_REF], Lemma 4.1, and [START_REF] Kinderlehrer | An introduction to variational inequalities and their applications[END_REF], Chapter VIII). Note that by considering the difference Ũ = U -u ϕ , we are able to derive a logarithmic upper bound, instead of a power of ξ, which would come up by considering U (see Theorem 2.4 of [START_REF] Lamberton | Brownian optimal stopping and random walks[END_REF]). For the proof of Theorem 3.1, we need some preliminary estimates on the derivatives

∂ 2 Ũ ∂x 2 and ∂ 2 Ũ ∂t∂x . Lemma 3.2. Fix T > 0 and j > 1.
For any ε ∈ (0, 1/4), there exists a constant C > 0 such that, for all t ∈ (0, T ],

∂ Ũ ∂x (t, .) L 2 (ν j ) ≤ C √ t and ∂ 2 Ũ ∂x 2 (t, .) L 2 (ν j ) ≤ Ct ε .
Proof: We know that Ũ solves the equation

- ∂ Ũ ∂t + (A -r) Ũ = h,
with initial condition Ũ (0, .) = 0, where the function h is given by

h(t, x) = (A -r)ϕ(x)1 {x≤ b(t)} , t > 0, x ∈ R.
We have the following identity (which can be viewed as a form of the early exercise premium formula). Ũ (t, .) = -t 0 e -r(t-s) p t-s * h(s, .)ds, where

p t (x) = 1 σ √ 2πt exp - (x + µt) 2 2σ 2 t = 1 σ √ t n x + µt σ √ t ,
with n denoting the standard normal density function. It is straightforward to check that

∂ Ũ ∂x (t, .) = - t 0 e -r(t-s) p t-s * ∂ h ∂x (s, .)ds,
and, with the notation δ z for the Dirac measure at a point z,

∂ h ∂x (t, x) = (A -r)ϕ ′ (x)1 {x≤ b(t)} -(A -r)ϕ(x)δb (t) (x) = -κ(t, x) + γ(t)δb (t) (x), (4)
with

κ(t, x) = -(A -r)ϕ ′ (x)1 {x≤ b(t)} and γ(t) = -(A -r)ϕ( b(t))
. Note that κ is a bounded function on (0, ∞) × R and γ is a continuous, nonnegative and bounded function on (0, +∞). At this stage, it is clear that

||p t-s * ∂ h ∂x (s, .)|| ∞ ≤ C/ √ t -s, so that ∂ Ũ ∂x (t, .) L 2 (ν j ) ≤ C √ t.
On the other hand, we have

∂ 2 Ũ ∂x 2 (t, .) L 2 (ν j ) ≤ t 0 e -r(t-s) p ′ t-s * κ(s, .) L 2 (ν j ) ds + ||ζ(t, .)|| L 2 (ν j ) , with ζ(t, .) = t 0 e -r(t-s) γ(s)p ′ t-s * δb (s) ds = t 0 e -r(t-s) γ(s)p ′ t-s (. -b(s))ds.
We have, using Lemma 3.1,

p ′ t-s * κ(s, .) L 2 (ν j ) = p ′ t-s (y)κ(s, . -y)dy L 2 (ν j ) ≤ |p ′ t-s (y)| ||κ(s, . -y)|| L 2 (ν j ) dy ≤ 2 j/4 ||κ(s, .)|| L 2 (ν j ) |p ′ t-s (y)|(1 + y 2 ) j/4 dy.
Note that, since κ is bounded and j > 1, sup s>0 ||κ(s, .)|| L 2 (ν j ) < ∞, so that, for some constant C > 0 (which may vary from line to line)

p ′ t-s * κ(s, .) L 2 (ν j ) ≤ C |p ′ t-s (y)|(1 + y 2 ) j/4 dy = C 1 σ 2 (t -s) n ′ y + µ(t -s) σ √ t -s (1 + y 2 ) j/4 dy = C 1 σ √ t -s |n ′ (z)| (1 + (-µ(t -s) + σ √ t -sz) 2 ) j/4 dz ≤ C √ t -s 1 + (t -s) j/2 . Hence, if 0 < t < T , t 0 e -r(t-s) p ′ t-s * κ(s, .) L 2 (ν j ) ds ≤ C t 0 ds √ t -s = 2C √ t.
We now estimate ||ζ(t, .)|| L 2 (ν j ) . We have, using the boundedness of γ,

|ζ(t, x)| = t 0 e -r(t-s) γ(s)p ′ t-s (x -b(s))ds ≤ C t 0 1 σ 2 (t -s) n ′ x -b(s) + µ(t -s) σ √ t -s ds Recall that n ′ (x) = -xn(x). Therefore |ζ(t, x)| ≤ C t 0 |x -b(s) + µ(t -s)| (t -s) 3/2 n x -b(s) + µ(t -s) σ √ t -s ds ≤ C t 0 ds √ t -s + C t 0 |x -b(s)| (t -s) 3/2 n x -b(s) + µ(t -s) σ √ t -s ds. Note that n(x 1 + x 2 ) = n(x 1 ) exp - x 2 2 2 -x 1 x 2 ≤ n(x 1 ) exp(-x 1 x 2 ) ≤ n(x 1 ) exp x 2 1 4 + x 2 2 = n(x 1 / √ 2)e x 2 2 .
Hence, for t ∈ (0, T ),

|ζ(t, x)| ≤ C t 0 ds √ t -s + C T t 0 |x -b(s)| (t -s) 3/2 n x -b(s) √ 2σ √ t -s ds.
Note that, for all α > 0, there exists

C α > 0, such that, for all y ∈ R, n(y/ √ 2) ≤ C α /|y| 2α . Hence, for t ∈ (0, T ), |ζ(t, x)| ≤ C √ t + C α t 0 |x -b(s)| (t -s) 3/2 (t -s) α |x -b(s)| 2α ds = C √ t + C α t 0 (t -s) α-3 2 |x -b(s)| 2α-1 ds = C √ t + C α t α-1 2 1 0 1 (1 -u) 3 2 -α |x -b(tu)| 2α-1 du. Now, take α = 1 2 + ε (with 0 < ε < 1/4) and put β(t, x) = |x -b(t)| 1-2α = |x -b(t)| -2ε . We get ||ζ(t, .)|| L 2 (ν j ) ≤ C √ t + Ct ε 1 0 1 (1 -u) 1-ε ||β(tu, .)|| L 2 (ν j ) du.
Using Lemma 3.1, we have

||β(tu, .)|| 2 L 2 (ν j ) ≤ 2 j/2 1 + b(tu) 2 j/2 1 |x| 4ε dx (1 + x 2 ) j/2 .
Since ε < 1/4, the integral on the righthand side is finite, and the lemma easily follows. ⋄

We now turn to the study of ∂ 2 Ũ ∂t∂x . Recall that ∂U/∂t solves the parabolic equation

-∂v/∂t + (A -r)v = 0 in the set {(t, x) | t > 0, x > b(t)}.
Since the exercise boundary is differentiable and ∂U/∂t is continuous and vanishes on the exercise boundary, it follows that ∂ 2 U ∂t∂x is continuous "up to the boundary", i.e. on the set {(t, x) | t > 0, x ≥ b(t)} (see [START_REF] Friedman | Parabolic variational inequalities in one space dimension and smoothness of the free boundary[END_REF], Lemma 4.5). We first show that ∂ 2 U ∂t∂x is nonnegative along the exercise boundary. Lemma 3.3. We have, for any t > 0,

∂ 2 U ∂t∂x (t, b(t)) ≥ 0.
Proof: We have, for all t > 0, due to the smooth fit property,

∂U ∂x (t, b(t)) = ϕ ′ ( b(t)),
so that, by differentiating with respect to t,

∂ 2 U ∂t∂x (t, b(t)) + ∂ 2 U ∂x 2 (t, b(t)) b′ (t) = ϕ ′′ ( b(t)) b′ (t) and ∂ 2 U ∂t∂x (t, b(t)) = - ∂ 2 U ∂x 2 (t, b(t)) -ϕ ′′ ( b(t)) b′ (t).
Observe that, for each t > 0, the function

x → U(t, x) -ϕ(x) is C 2 on the interval [ b(t), ∞)
and has a minimum at b(t). Therefore, its second derivative must be nonnegative at this point. Since b′ (t) ≤ 0, the lemma is proved. ⋄ Lemma 3.4. Fix T > 0 and j > 1. There exists a constant C > 0 such that, for all

t 1 ∈ (0, T ∧ 1], T t 1 ∂ 2 Ũ ∂t∂x (t, .) 2 L 2 (ν j ) dt ≤ C ln(1/t 1 ).
For the proof of Lemma 3.4, we will need the bilinear form associated with the operator A -r.

We first introduce the relevant weighted Sobolev spaces. For j > 1, let

H j = L 2 (R, ν j ) and V j = {f ∈ H j | f ′ ∈ H j }.
The inner product on H j will be denoted by (•, •) j and the associated norm by | • | j . The natural norm on V j will be denoted by || • || j . Thus, we have

|f | 2 j = +∞ -∞ f 2 (x) dx (1 + x 2 ) j/2 ,
and ||f || 2 j = |f | 2 j + |f ′ | 2 j .
Recall that the partial differential operator A is defined by

A = σ 2 2 ∂ 2 ∂x 2 + µ ∂ ∂x .
We associate with the operator A -r a bilinear functional on V j , defined by

a j (f, g) = σ 2 2 ∞ -∞ f ′ (x)g ′ (x) dx (1 + x 2 ) j/2 - jσ 2 2 ∞ -∞ f ′ (x)g(x) x (1 + x 2 ) (j/2)+1 dx -µ ∞ -∞ f ′ (x)g(x) dx (1 + x 2 ) j/2 + r ∞ -∞ f (x)g(x) dx (1 + x 2 ) j/2 , so that, if f ′ ∈ V j , a j (f, g) = -((A -r)f, g) j .
It will be convenient to write a j (f, g) as a j (f, g) = ãj (f, g) + āj (f, g), with

ãj (f, g) = σ 2 2 [(f ′ , g ′ ) j + (f, g) j ] and āj (f, g) = a j (f, g) -ãj (f, g). ( 5
)
With these notations, it is easy to check that |ā j (f, g)| ≤ C||f || j |g| j and |ā j (f, g)| ≤ C||g|| j |f | j , for some constant C which does not depend on f nor g.

- ∂W m ∂t + (A -r)W m = ∂h m ∂x .
Multiply by ∂W m /∂t and integrate with respect to ν j to get, for any fixed t > 0,

- ∂W m ∂t (t, .), ∂W m ∂t (t, .) j -a j W m (t, .), ∂W m ∂t (t, .) = ∂h m ∂x (t, x) ∂W m ∂t (t, x)ν j (dx).
Note that

a j W m (t, .), ∂W m ∂t (t, .) = ãj W m (t, .), ∂W m ∂t (t, .) + āj W m (t, .), ∂W m ∂t (t, .) = 1 2 d dt (ã j (W m (t, .), W m (t, .))) + āj W m (t, .), ∂W m ∂t (t, .) .
By integrating with respect to time, we get, if 0 < t 1 < T ,

- T t 1 ∂W m ∂t (t, .) 2 j dt + 1 2 (ã j (W m (t 1 , .), W m (t 1 , .)) -ãj (W m (T, .), W m (T, .))) = T t 1 āj W m (t,
∂h m ∂x (t, x) = -κ m (t, x) + γ m (t, x),
where κ m = κ * ρ m , and κ is a bounded function, and

γ m (t, x) = ρ m (t -τ, x -b(τ ))γ(τ )dτ.

Hence

∂t = ∂ 2 U m ∂t∂x - ∂ 2 u m ∂t∂x ,
where

U m = U * ρ m and u m = u ϕ * ρ m , so that - T t 1 γ m (t, .), ∂W m ∂t (t, .) j dt = J (1) m + J (2) m ,
with

J (1) m = - T t 1 γ m (t, .), ∂ 2 U m ∂t∂x (t, .) j dt and J (2) m = T t 1 γ m (t, .), ∂ 2 u m ∂t∂x (t, .) j dt.
We have

J (1) m = - T t 1 dt dx (1 + x 2 ) j/2 γ m (t, x) ∂ 2 U m ∂t∂x (t, x) = - T t 1 dt dx (1 + x 2 ) j/2 dτ γ(τ )ρ m (t -τ, x -b(τ )) ∂ 2 U m ∂t∂x (t, x) = -dτ ds dy1 {t 1 <τ +s<T } γ(τ )ρ m (s, y) ∂ 2 U m ∂t∂x (τ + s, b(τ ) + y) 1 (1 + (y + b(τ )) 2 ) j/2 = - T + 1 m t 1 γ(τ )η m (τ )dτ, where η m (τ ) = T -τ t 1 -τ ds dy (1 + (y + b(τ )) 2 ) j/2 ρ m (s, y) ∂ 2 U m ∂t∂x (τ + s, b(τ ) + y) = T -τ t 1 -τ ds ρ m (s, y)dy (1 + (y + b(τ )) 2 ) j/2 ds ′ dy ′ ρ m (s ′ , y ′ ) ∂ 2 U ∂t∂x (τ + s -s ′ , b(τ ) + y -y ′ ). Note that ∂ 2 U ∂t∂x = 0 on the open set S = {(t, x) | t > 0,
x < b(t)} (which is the interior set of the stopping region), so that

η m (τ ) = T -τ t 1 -τ ds ρ m (s, y)dy (1 + (y + b(τ )) 2 ) j/2 ds ′ dy ′ ρ m (s ′ , y ′ ) W (τ, s -s ′ , y -y ′ ), where W (τ, θ, z) = ∂ 2 U ∂t∂x (τ + θ, b(τ ) + z)1 { b(τ )+z≥ b(τ +θ)} . Since ∂ 2 U ∂t∂x (τ, b(τ )) ≥ 0, for τ > 0, we have η m (τ ) ≥ T -τ t 1 -τ ds ρ m (s, y)dyρ m (s ′ , y ′ )ds ′ dy ′ 1 + (y + b(τ )) 2 j/2 D(τ, s -s ′ , y -y ′ ),
where

D(τ, θ, z) = ∂ 2 U ∂t∂x (τ + θ, b(τ ) + z) - ∂ 2 U ∂t∂x (τ, b(τ )) 1 { b(τ )+z≥ b(τ +θ)} .
Hence (since γ ≥ 0)

J (1) m ≤ - T + 1 m t 1 γ(τ )ε m (τ )dτ, with ε m (τ ) = T -τ t 1 -τ ds ρ m (s, y)dyρ m (s ′ , y ′ )ds ′ dy ′ (1 + (y + b(τ )) 2 ) j/2 D(τ, s -s ′ , y -y ′ ). We have |ε m (τ )| ≤ T -τ t 1 -τ D m (τ )ds = (T -t 1 )D m (τ ),
where

D m (τ ) = sup |θ|≤1/m,|z|≤2/m |D(τ, θ, z)| 1 { b(τ )+z≥ b(τ +θ)} .
Due to the continuity properties of ∂ 2 U ∂t∂x , as m → ∞, the function D m converges to 0, uniformly on the interval [t 1 , T + 1]. Therefore, we have

lim sup m→∞ J (1) m ≤ 0. (7)
We now examine J (2) m . We have, using the boundedness of γ,

|J (2) m | ≤ T t 1 dt dx (1 + x 2 ) j/2 |γ m (t, x)| ∂ 2 u m ∂t∂x (t, x) ≤ C T t 1 dt dx (1 + x 2 ) j/2 dτ ρ m (t -τ, x -b(τ )) ∂ 2 u m ∂t∂x (t, x) .
Note that, since ϕ is Lipschitz, we have

∂ 2 uϕ ∂t∂x (t, .) ∞ ≤ C t and, since supp ρ ⊂ [-1, 0] × [-1, +1], ∂ 2 u m ∂t∂x (t, .) ∞ ≤ dτ dyρ m (τ, y) ∂ 2 u ϕ ∂t∂x (t -τ, .) ∞ ≤ C t .
Hence

|J (2) m | ≤ C T t 1 dt t dx dτ ρ m (t -τ, x -b(τ )) = C ln T t 1 . ( 8
)
It follows from ( 6), ( 7) and ( 8) that

lim sup m→∞ T t 1 ∂W m ∂t (t, .) 2 j dt ≤ C 1 + t 2ε 1 + ln T t 1 ,
which proves the lemma. ⋄

Proof of Theorem 3.1

For the proof of Theorem 3.1, we will work on the equation satisfied by ∂ Ũ /∂t. Let

V = ∂ Ũ ∂t .
We have

- ∂V ∂t + (A -r)V = ∂ h ∂t , where h(t, x) = (A -r)ϕ(x)1 {x≤ b(t)} , t > 0, x ∈ R.
The following lemma will clarify the computation of the derivative ∂ h/∂t in the sense of distributions.

Lemma 3.5. Define the function I on (0, +∞) × R by

I(t, x) = 1 {x≤ b(t)} , t > 0, x ∈ R.
The distribution ∂I/∂t applied to a compactly supported

C ∞ function ρ on (0, +∞) × R is given by ∂I ∂t , ρ = b′ (t)ρ(t, b(t))dt.
This can be written (less precisely): ∂I ∂t (t, .) = b′ (t)δb (t) .

Proof: We have

∂I ∂t , ρ = -I, ∂ρ ∂t = -dt dxI(t, x) ∂ρ ∂t (t, x)
Let J be the range of b. We have J = ( b(∞), b(0)). Note that, if x ≤ b(∞), I(t, x) = 1 for all t > 0 and if x ≥ b(0) I(t, x) = 0 for all t > 0, so that, in both cases, I(t, x) ∂ρ ∂t (t, x)dt = 0. Therefore

∂I ∂t , ρ = - J dx dt1 {x≤ b(t)} ∂ρ ∂t (t, x) = - J dx dt1 {t≤ b-1 (x)} ∂ρ ∂t (t, x) = - J dxρ( b-1 (x), x) = b′ (t)ρ(t, b(t))dt.
Here, we have used the fact that b is strictly decreasing (which is proved in [START_REF] Villeneuve | Options américaines dans un modèle de Black-Scholes multidimensionnel[END_REF]), but we can also approximate b by the strictly decreasing functions bε (t) = -εt + b(t) to derive the formula. In fact, we only need b to be C 1 : indeed, we can replace b(t) by bµ (t) = -µt + b(t)

and choose µ so that bµ is strictly increasing in a neighborhood of the time projection of the support of ρ. ⋄

We now proceed with the proof of Theorem 3.1. As in the proof of Lemma 3.4, we introduce a regularizing sequence ρ m , and set

V m = V * ρ m and χ m = ∂ h ∂t * ρ m , so that - ∂V m ∂t + (A -r)V m = χ m
Note that the functions V m , χ m are C ∞ , with bounded derivatives on any subset [t 1 , T ]×R, with 0 < t 1 < T . This is due to the fact that V is bounded on such subsets. For any fixed t > 0, multiply by ∂V m /∂t and integrate with respect to ν j to get

- ∂V m ∂t (t, .) 2 j -a j V m (t, .), ∂V m ∂t (t, .) = χ m (t, x) ∂V m ∂t (t, x)ν j (dx).
We have

a j V m (t, .), ∂V m ∂t (t, .) = 1 2 d dt (ã j (V m (t, .), V m (t, .))) + āj V m (t, .), ∂V m ∂t (t, .) .
By integrating with respect to time, we get, if 0 < t 1 < T ,

- T t 1
∂V m ∂t (t, .)

2 j dt + 1 2 [ã j (V m (t 1 , .), V m (t 1 , .)) -ãj (V m (T, .), V m (T, .))] = T t 1 āj V m (t, .), ∂V m ∂t (t, .) dt + T t 1 χ m (t, .), ∂V m ∂t (t, .) j dt. Hence T t 1 ∂V m ∂t (t, .) 2 j dt ≤ C ||V m (t 1 , .)|| 2 j - T t 1 āj V m (t, .), ∂V m ∂t (t, .) dt - T t 1 χ m (t, .), ∂V m ∂t (t, .) j dt ≤ C   ||V m (t 1 , .)|| 2 j + T t 1 ||V m (t, .)|| j ∂V m ∂t (t, .) j dt   + J m (t 1 , T ), with J m (t 1 , T ) = - T t 1 χ m (t, .), ∂V m ∂t (t, .) j dt.
Using the inequality

2 ||V m (t, .)|| j ∂V m ∂t (t, .) j ≤ ε ∂V m ∂t (t, .) 2 j + 1 ε ||V m (t, .)|| 2 j ,
we derive

1 2 T t 1 ∂V m ∂t (t, .) 2 j dt ≤ C ||V m (t 1 , .)|| 2 j + T t 1 ||V m (t, .)|| 2 j dt + J m (t 1 , T ). (9)
We now study J m (t 1 , T ). Note that, for any fixed t > 0,

χ m (t, .), ∂V m ∂t (t, .) j = ν j (dx) ∂V m ∂t (t, x) ∂ h ∂t * ρ m (t, x) We have ∂ h ∂t (t, .) = (A -r)ϕ ∂I ∂t (t, .),
so that, using Lemma 3.5, and the notation

γ(t) = -(A -r)ϕ( b(t)) ∂ h ∂t * ρ m (t, x) = -dτ ρ m (t -τ, x -b(τ )) b′ (τ )γ(τ ), Recall that γ(τ ) ≥ 0. Hence χ m (t, .), ∂V m ∂t (t, .) j = -ν j (dx) ∂V m ∂t (t, x) dτ ρ m (t -τ, x -b(τ )) b′ (τ )γ(τ ) = -dτ dx (1 + x 2 ) j/2 ∂V m ∂t (t, x)ρ m (t -τ, x -b(τ )) b′ (τ )γ(τ ) = -dτ dy (1 + (y + b(τ )) 2 ) j/2 ∂V m ∂t (t, y + b(τ ))ρ m (t -τ, y) b′ (τ )γ(τ ).
Going back to J m (t 1 , T ), we have

J m (t 1 , T ) = - T t 1 dt dτ dy ∂V m ∂t (t, y + b(τ ))ρ m (t -τ, y)γ j (τ, y), with γj (τ, y) = - 1 (1 + (y + b(τ )) 2 ) j/2 b′ (τ )γ(τ ).
Note that γj (τ, y) ≥ 0. We have

J m (t 1 , T ) = -dτ dt dy1 {t 1 <t<T } ∂V m ∂t (t, y + b(τ ))ρ m (t -τ, y)γ j (τ, y) = -dτ ds dy1 {t 1 <τ +s<T } ∂V m ∂t (τ + s, y + b(τ ))ρ m (s, y)γ j (τ, y).
Observe that

d dτ V m (τ + s, y + b(τ )) = ∂V m ∂t (τ + s, y + b(τ )) + ∂V m ∂x (τ + s, y + b(τ )) b′ (τ ), so that ∂V m ∂t (τ + s, y + b(τ )) = d dτ V m (τ + s, y + b(τ )) - ∂V m ∂x (τ + s, y + b(τ )) b′ (τ ).
Hence

J m (t 1 , T ) = Ĵm (t 1 , T ) + Jm (t 1 , T ), with Ĵm (t 1 , T ) = -dτ ds dy1 {t 1 <τ +s<T } d dτ V m (τ + s, y + b(τ )) ρ m (s, y)γ j (τ, y)
and

Jm (t 1 , T ) = + dτ ds dy1 {t 1 <τ +s<T } ∂V m ∂x (τ + s, y + b(τ )) b′ (τ )ρ m (s, y)γ j (τ, y).
We have, using integration by parts,

Ĵm (t 1 , T ) = -ds dyρ m (s, y) T -s t 1 -s d dτ V m (τ + s, y + b(τ )) γj (τ, y)dτ = -ds dyρ m (s, y)V m (T, y + b(T -s))γ j (T -s, y) + ds dyρ m (s, y)V m (t 1 , y + b(t 1 -s))γ j (t 1 -s, y) + ds dyρ m (s, y) T -s t 1 -s V m (s + τ, y + b(τ )) ∂γ j ∂τ (τ, y)dτ.
Note that, due to the continuity of V (= ∂ Ũ /∂t) on (0, ∞) × R, the sequence V m converges uniformly to V on compact sets. We also have the continuity of γj and ∂γ j /∂τ (due to the fact that b is C 2 ). We easily deduce thereof that

lim m→∞ Ĵm (t 1 , T ) = -V (T, b(T ))γ j (T, 0) + V (t 1 , b(t 1 ))γ j (t 1 , 0) + T t 1 V (τ, b(τ )) ∂γ j ∂τ (τ, 0)dτ,( 10
)
and the convergence is uniform with respect to t 1 , as long as t 1 remains in a compact set of the form [ξ, T ], where 0 < ξ < T . For Jm (t 1 , T ), we have

Jm (t 1 , T ) = J(1) m (t 1 , T ) + J(2) m (t 1 , T ), with J(1) m (t 1 , T ) = + dτ ds dy1 {t 1 <τ +s<T } ∂ 2 U m ∂t∂x (τ + s, y + b(τ )) b′ (τ )ρ m (s, y)γ j (τ, y)
and

J(2) m (t 1 , T ) = -dτ ds dy1 {t 1 <τ +s<T } ∂ 2 u m ∂t∂x (τ + s, y + b(τ )) b′ (τ )ρ m (s, y)γ j (τ, y).
We deal with J(1) m (t 1 , T ) in the same way as for the proof of [START_REF] Lamberton | Critical price near maturity for an American option on a dividend paying stock (avec S. Villeneuve)[END_REF]. Using the fact that b′ (τ )γ j (τ, y) ≤ 0, we have J( 1)

m (t 1 , T ) ≤ J(1) m (t 1 , T ), with J(1) m (t 1 , T ) = dτ ds dy1 {t 1 <τ +s<T } ds ′ dy ′ ρ m (s ′ , y ′ )D(τ, s -s ′ , y -y ′ ) b′ (τ )ρ m (s, y)γ j (τ, y),
where

D(τ, θ, z) = ∂ 2 U ∂t∂x (τ + θ, b(τ ) + z) - ∂ 2 U ∂t∂x (τ, b(τ )) 1 { b(τ )+z≥ b(τ +θ)} .
Due to the continuity properties of

Jm (t 1 , T ) -J(t 1 , T ) = 0,
where

J(t 1 , T ) = -V (T, b(T ))γ j (T, 0) + V (t 1 , b(t 1 ))γ j (t 1 , 0) + T t 1 V (τ, b(τ )) ∂γ j ∂τ (τ, 0)dτ - T t 1 dτ ∂ 2 u ∂t∂x (τ, b(τ )) b′ (τ )γ j (τ, 0).
Since ∂U/∂t vanishes along the exercise boundary, we have

V (t, b(t)) = -∂uϕ ∂t (t, b(t)), so that -V (T, b(T ))γ j (T, 0) + V (t 1 , b(t 1 ))γ j (t 1 , 0) + T t 1 V (τ, b(τ )) ∂γ j ∂τ (τ, 0)dτ = ∂u ϕ ∂t (T, b(T ))γ j (T, 0) - ∂u ϕ ∂t (t 1 , b(t 1 ))γ j (t 1 , 0) - T t 1 ∂u ϕ ∂t (τ, b(τ )) ∂γ j ∂τ (τ, 0)dτ = T t 1 d dτ ∂u ϕ ∂t (τ, b(τ )) γj (τ, 0)dτ, so that J (t 1 , T ) = T t 1 d dτ ∂u ϕ ∂t (τ, b(τ )) - ∂ 2 u ϕ ∂t∂x (τ, b(τ )) b′ (τ ) γj (τ, 0)dτ = T t 1 ∂ 2 u ϕ ∂t 2 (τ, b(τ ))γ j (τ, 0)dτ.
We now go back to [START_REF] Silvestrov | American Type Options. Stochastic Approximation Methods[END_REF] and integrate with respect to t 1 to derive where the last inequality follows from Lemma 3.4. Moreover,

1 2 T ξ dt 1 T t 1 ∂ 2 Ũm ∂t 2 (t, .) 2 j dt ≤ C T ξ ||V m (t 1 , .)||
lim m→∞ T ξ dt 1 Jm (t 1 , T ) = T ξ dt 1 J (t 1 , T ) = T ξ (t -ξ) ∂ 2 u ϕ ∂t 2 (t, b(t))γ j (t, 0)dt.
Hence

1 2 T ξ (t -ξ) ∂ 2 Ũ ∂t 2 (t, .) 2 j dt ≤ C 1 + ln T ξ + T ξ (t -ξ) ∂ 2 u ϕ ∂t 2 (t, b(t))γ j (t, 0)dt.
Theorem 3.1 now follows from the following lemma, which relies on the asymptotic behavior of the exercice boundary near maturity (see [START_REF] Barles | Critical Stock Price near expiration[END_REF], [START_REF] Lamberton | Critical price near maturity for an American option on a dividend paying stock (avec S. Villeneuve)[END_REF]).

Lemma 3.6. We have

T ξ (t -ξ) ∂ 2 u ϕ ∂t 2 (t, b(t)) | b′ (t)|dt ≤ C 1 + | ln ξ| β , with β =      3/2, if d ≤ r, 1, if d > r.
Proof: We first note that, since ϕ is bounded and Lipschitz continuous, we have

∂ 2 u ϕ ∂t 2 (t, .) ∞ ≤ C t 3/2 .
This can be seen by arguing, as in the proof of Proposition 3.1, that u ϕ (t, .) = e -rt p t * ϕ, so that and ∂ 2 uϕ ∂t 2 (t, .) = (A -r) 2 u ϕ . In order to estimate the x-derivatives of u ϕ up to the order 4, we may differentiate p t three times and use the boundedness of ϕ ′ .

We then have -P 0 , we relate this quantity to the modified value function u using (1) as follows:

T ξ (t -ξ) ∂ 2 u ϕ ∂t 2 (t, b(t)) | b′ (t)|dt ≤ C T ξ t -ξ t 3/2 | b′ (t)|dt ≤ C T ξ 1 √ t | b′ (t)|dt. Now, since b′ (t) ≤ 0, we have T ξ 1 √ t | b′ (t)|dt = - T ξ 1 √ t b′ (t)dt = - b(T ) -b(0) √ T - b(ξ) -b(0) √ ξ - 1 2 T ξ 1 t 3/2 ( b(t) -b(0))dt ≤ b(0) -b(T ) √ T + 1 2 T ξ 1 t 3/2 ( b(0) -b(t))dt 22 If d ≤ r, we have b(0) -b(t) ≤ C t| ln t| for t close to 0, so that T ξ 1 t 3/2 ( b(0) -b(t))dt ≤ C T ξ 1 t | ln t|dt ≤ C(1 + | ln ξ| 3/2 ). If d > r, we have b(0) -b(t) ≤ C √ t for t close to 0, so that T ξ 1 t 3/2 ( b(0) -b(t))dt ≤ C T ξ 1 t dt = C ln(T /ξ). ⋄ 4
P (n) 0 -P 0 = sup τ ∈T (n) 0,T E e -rτ g(µ 0 τ + B (n) τ ) -u(0, 0) ≤ sup τ ∈T (n) 0,T E u(τ, B (n) τ ) -u(0, 0) = sup τ ∈T (n) 0,T E   τ /h j=1 Du((j -1)h, B (n) (j-1)h )   .
We observe that Du ≤ Du, and recall from [START_REF] Lamberton | Brownian optimal stopping and random walks[END_REF] (Lemma 4.1) that sup 0≤j≤n-1 E Du(jh, B

(n) jh ) ≤ Ch, so that

P (n) 0 -P 0 ≤ E   n-2 j=1 Du(jh, B (n) jh )   + O(h) ≤ C k,X h √ 2 T -h h ds √ s dy 1 + |y| k ∂ 3 u ∂t∂x 2 (s, y) + O(h). (11)
Here, we have a regularity problem, since u is not C 3 . This problem can be fixed as follows. By convolution, one can approximate u by a sequence u m which is smooth, uniformly bounded and satisfies δu m ≤ 0, and Du m ≤ Du m . We need the following variant of Lemma 3.1. Lemma 4.1. If ρ is a Radon measure on (0, T )×R and q a nonnegative integrable function on (0, T ) × R, with q(t, x) = 0 for t / ∈ (0, a), where a satisfies 0 < a < h, we have

T -h h ds √ s |ρ * q(s, y)| (1 + |y| 2 ) k/2 dy ≤ 2 k/2 T -h h-a R |ρ(dt, dz)| √ t(1 + z 2 ) k/2 a 0 ds ∞ -∞ q(s, x)(1 + x 2 ) k/2 dx. Proof: We have T -h h ds √ s dy (1 + |y| 2 ) k/2 |ρ * q(s, y)| ≤ T -h h ds √ s dy (1 + |y| 2 ) k/2 |ρ(dt, dz)| q(s -t, y -z) = |ρ(dt, dz)| T -h h ds √ s dy (1 + |y| 2 ) k/2 q(s -t, y -z) ≤ T -h h-a |ρ(dt, dz)| a 0 dθ √ t + θ dy (1 + |y| 2 ) k/2 q(θ, y -z) ≤ T -h h-a |ρ(dt, dz)| √ t a 0 dθ dy (1 + |y| 2 ) k/2 q(θ, y -z) = T -h h-a |ρ(dt, dz)| √ t(1 + z 2 ) k/2 a 0 dθ (1 + z 2 ) k/2 dx (1 + |x + z| 2 ) k/2 q(θ, x) ≤ 2 k/2 T -h h-a |ρ(dt, dz)| √ t(1 + z 2 ) k/2 a 0 dθ dx(1 + x 2 ) k/2 q(θ, x),
where the last inequality follows from 1+z 2 1+(x+z) 2 ≤ 2(1 + x 2 ). ⋄ Using Lemma 4.1 and the fact that ∂ 3 u/(∂t∂x 2 ) is a Radon measure (see (13) and the comment below), we derive the correct version of (11), namely

P (n) 0 -P 0 ≤ C k,X h √ 2 T -h h 1 √ s 1 1 + |y| k ∂ 3 u ∂t∂x 2 (ds, dy) + O(h). ( 12
)
If we introduce the function ũ := u -ū, we have, using the fact that δū = 0,

T -h h 1 √ s 1 1 + |y| k ∂ 3 u ∂t∂x 2 (ds, dy) ≤ T -h h 1 √ s 1 1 + |y| k ∂ 3 ũ ∂t∂x 2 (ds, dy) +2 T -h h ds √ s dy 1 + |y| k ∂ 2 ū ∂t 2 (s, y) ≤ T -h h 1 √ s 1 1 + |y| k ∂ 3 ũ ∂t∂x 2 (ds, dy) +2C T T -h h ds √ s(T -s) ≤ T -h h 1 √ s 1 1 + |y| k ∂ 3 ũ ∂t∂x 2 (ds, dy) + C T | ln h|,
where we have used Proposition 3.1.

We now need to estimate T -h The modification of τ 1 into τ is motivated by the unboundedness of ∂u/∂t near T . We have, due to the definition of τ 1 , P (n) -P ≥ E e -rτ g(µ 0 τ + B (n) τ ) -u(0, 0) = E e -rτ g(µ 0 τ + B (n) τ ) -u(τ, B (n) τ ) + u(τ, B (n) τ ) -u(0, 0) .

h 1 √ s 1 1+|y| k ∂ 3 ũ ∂t∂x 2 (ds, dy) . Recall from Remark 2.1 that ∂ ũ ∂t (t, x) + 1 2 ∂ 2 ũ ∂x 2 (t, x) = ζ(t,
We have

E

Proof of Lemma 3 . 4 :

 34 In order to rule out regularity issues, we introduce a C ∞ , nonnegative function ρ on R × R, with ρ(t, x)dtdx = 1 and supp ρ ⊂ [-1, 0] × [-1, +1] and set, for any positive integer m, ρ m (t, x) = m 2 ρ(mt, mx).Now, let W = ∂ Ũ ∂x , W m = W * ρ m , and h m = h * ρ m .For each m > 0, the functions W m , h m are C ∞ with bounded derivatives and we have

  In order to derive an upper bound for P (n) 0

  u(τ, B (n) τ ) -u(0, 0) = E equality follows from Lemma 4.1 of[START_REF] Lamberton | Brownian optimal stopping and random walks[END_REF]. Now, if j < (τ /h) ∧ (n -2), we haved(B (n) jh , I jh+h ) > √ h||X|| ∞ + |µ 0 |h, so that B (n) jh > b(jh + h) + √ h||X|| ∞ + |µ 0 |h.We then have, for s ∈ [jh, jh + h] and z ∈ [0,jh + zX > b(jh + h) + zX + √ h||X|| ∞ + |µ 0 |h ≥ b(jh + h) + |µ 0 |h = b(s) + b(jh + h) + µ 0 (jh + h) -( b(s) + µ 0 s) -µ 0 (jh + h -s) + |µ 0 |h ≥ b(s),the last inequality coming from the fact that t → b(t) + µ 0 t is increasing. We can now assert that, for j < (τ /h) ∧ (n -2), Du(B(n) jh , jh) = Du(B (n) jh , jh), so that E u(τ, B (n) τ ) -u(0, 0) ≤ E α = 1 if d >r and α = 5/4 if d ≤ r, as follows from the discussion in the previous section.

  ||W (t, .)|| j ≤ Ct ε . Hence ãj (W m (t 1 , .), W m (t 1 , .)) ≤ C ρ m (t 1 -t, y)2 j/2 (1 + y 2 ) j/2 t 2ε dtdy = C ρ m (t, y)2 j/2 (1 + y 2 ) j/2 |t 1 -t| 2ε dtdy

	We have
	ãj (W Using Lemma 3.2, we have (for ε ∈ (0, 1/4)) = C ρ(t, y)2 j/2 1 +	y 2 m 2	j/2	t 1 -	t m	2ε	dtdy
								≤ C(1 + t 2ε 1 ).
	We also have
	T t 1	||W m (t, .)|| 2 j dt =	T t 1				j dsdyρ 2	dt
	.), ãj (W T ∂W m ∂t (t, .) dt + T t 1 t 1 -1 2 T t 1 ∂h m ∂x (t, .), ∂W m (t, .) dt āj W m (t, .), ∂h m ∂x (t, .), ∂W m ∂t ∂W m (t, .) ∂t (t, .) dt ≤ T t 1 dsdyρ Since T +1 Hence T t 1 ∂W m ∂t (t, .) 2 j dt ≤ 0 ||W (s, .)|| 2 j ds < ∞, we deduce that ∂t j ≤ 1 2 ãj (W T t 1 ||W m (t, .)|| j ∂W m j ∂t (t, .) dt 1 2 T t 1 ∂W m ∂t (t, .) 2 j dt ≤ C 1 + t 2ε 1 -T t 1 ∂h m ∂x (t, .), ∂W m ∂t (t, .)
						-	T t 1		∂h m ∂x	(t, .),	∂W m ∂t	(t, .)
	Using the inequality
			2 ||W m (t, .)|| j	∂W m ∂t	(t, .)	j	≤ ε	∂W m ∂t	(t, .)	2 j	+	1 ε	||W m (t, .)|| 2 j ,
	we get
	1 2	T t 1	∂W m ∂t	(t, .)	2 j	dt ≤	1 2	ãj (W m (t 1 , .), W m (t 1 , .)) + C	T t 1	||W m (t, .)|| 2 j dt
								-	T t 1	∂h m ∂x	(t, .),	∂W m ∂t	(t, .)

j dt. m (t 1 , .), W m (t 1 , .))m (t 1 , .), W m (t 1 , .)) + C j dt. j dt. m (t 1 , .), W m (t 1 , .)) ≤ C ||W m (t 1 , .)||

2 j and, using Lemma 3.1, ||W m (t 1 , .)|| j ≤ ρ m (t 1 -t, y) ||W (t, . -y)|| j dtdy ≤ ρ m (t 1 -t, y)2 j/4 (1 + y 2 ) j/4 ||W (t, .)|| j dtdy. m (t -s, y)W (s, . -y) m (t -s, y) ||W (s, . -y)|| 2 j dt ≤ T t 1 dt dsdyρ m (t -s, y)2 j/2 (1 + y 2 ) j/2 ||W (s, .)|| 2 j ≤ T + 1 m t 1 ds ||W (s, .)|| 2 j dtdyρ m (t -s, y)2 j/2 (1 + y 2 ) j/2 . j dt. It follows from the proof of Lemma 3.2 (see (4)) that

  and the convergence is uniform with respect to t 1 , as long as t 1 remains in [ξ, T ]. On the other hand, due to the continuity of ∂ 2 uϕ ∂t∂x , we have uniformly with respect to t 1 ∈ [ξ, T ]. At this stage, we can state that J m (t 1 , T ) ≤ Jm (t 1 , T ),

			∂ 2 U ∂t∂x , we have
			lim m→∞	J(1) m (t 1 , T ) = 0,
	lim m→∞	J(2) m (t 1 , T ) = -	T t 1	dτ	∂ 2 u ϕ ∂t∂x	(τ, b(τ )) b′ (τ )γ j (τ, 0),
	with Jm (t 1 , T ) = Ĵm (t 1 , T ) + J(1) m (t 1 , T ) + J(2) m (t 1 , T ), and
		lim m→∞	sup t 1 ∈[ξ,T ]

  Moreover, | b′ (t)| ≤ | b′ (T -t)| + |µ/σ|, so thatTherefore, the righthand side of (13) is a Radon measure and since, due to Theorem 3.1, ∂ 2 ũ/∂t 2 is locally integrable, it follows that ∂ 3 ũ/∂t∂x 2 is a Radon measure.Moreover, we have∂ 3 ũ ∂t∂x 2 (ds, dy) ≤ C T + 2Now, using the Cauchy-Schwarz inequality and Theorem 3.1, we have The last inequality follows from Theorem 3.1 and Lemma 3.4, and the connection between the derivatives of the functions Ũ and ũ (see Remark 2.1; we also use the classical bounds ||∂U/∂t(t, .)|| ∞ + ||∂ 2 U/∂x 2 (t, .)|| ∞ ≤ C/ √ t).For the derivation of the lower bound, we use the stopping time introduced in[START_REF] Lamberton | Brownian optimal stopping and random walks[END_REF] (see the proof of Theorem 5.6). Namelyτ = τ 1 1 {τ 1 <T -h} + T 1 {τ 1 =T -h} ,whereτ 1 = inf t ∈ [0, T -h] | t/h ∈ N and d(B , I t+h ) ≤ √ h||X|| ∞ + |µ 0 |h .Here, I t = {x ∈ R | u(t, x) = g(t, x + µ 0 t)}. Note that, if t < T , I t = (-∞, b(t)].

	T We conclude that 0 dt √ t | b′ (t)| ≤ ≤ = T -h h 1 √ s 1 + |y| k T /2 0 sup 0≤t≤T /2 dt √ t | b′ (T -t)| | b′ (T -t)| + T /2 0 sup 0≤t≤T /2 | b′ (T -t)| T /2 0 1 T -h h ds √ s dy 1 + |y| k ∂ 2 ũ ∂t 2 (s, y) ≤ C T -h h ds s(T -s -h 2 ) 1/2   T -h T /2 T dt √ t dt √ t h ds(T -s -dt √ t | b′ (T -t)| + 2|µ/σ| + + 2 T T T /2 | b′ (T -t)|dt + 2|µ/σ| √ T + 2 T b(0) -b(T /2) + 2|µ/σ| √ T < ∞. √ T T -h h ds √ s 1 1 + |y| k ∂ 2 ũ ∂t 2 (s, y) . h 2 ) ∂ 2 ũ ∂t 2 (s, .) 2 k   1/2 ≤ C | ln h|   T -h 2 h ds(T -s -h 2 ) ∂ 2 ũ ∂t 2 (s, .) 2 k   1/2 = C | ln h|   T -h h/2 dt(t -h 2 ) ∂ 2 ũ ∂t 2 (T -t, .) 2 k   1/2 P (n) 0 -P 0 ≤ C (ln n) α n , with α = 1 if d > r and α = 5/4 if d ≤ r. 5 Lower bound for P (n) 0 -P 0 (13) where we have used Lemma 3.5. Note that sup x≤ b(t) ∂ζ ∂t (t, x) < ∞ and sup (n)

x)1 {x≤ b(t)} ,

where ζ(t, x) = e -rt (A -r)ϕ(ln(S 0 ) + µt + σx) and b(t) = b(T -t) -µt -ln(S 0 ) /σ. By differentiating wit respect to t, we derive the following expression

∂ 2 ũ ∂t 2 (t, x) + 1 2 ∂ 3 ũ ∂t∂x 2 (t, x) = ∂ζ ∂t (t, x)1 {x≤ b(t)} + ζ(t, b(t)) b′ (t)δ b(t) , 0<t<T ζ(t, b(t)) < ∞. ≤ C | ln h| 1+β , with β = 1 if d > r and β = 3/2 if d ≤ r. t

We now want a lower bound for E e -rτ g(µ 0 τ + B (n) τ ) -u(τ, B (n) τ ) .

We have, using the equality {τ ≥ T -h} = {τ = T },

Using the estimate ∂u ∂t (t, .)

The estimate

is now an easy consequence of Lemma 5.7 and Remark 5.8 of [START_REF] Lamberton | Brownian optimal stopping and random walks[END_REF], which can be summarized in the following statement.


Alejandro Fernandez-Martinez

Gabriel J Cuello
email: cuello@ill.fr

Gabriela Roman-Ross

Mark R Johnson

Fabrizio Bardelli

Xavier Turrillas

Laurent Charlet

A Fernández-Martínez

G Román-Ross

Arsenic uptake by gypsum and calcite: Modeling and

Keywords: arsenic, minerals, simulation, diffraction, EXAFS

come

Introduction

Arsenic is recognized as a dangerous pollutant of the environment [START_REF] Nriagu | Arsenic in the Environment[END_REF]. Arsenic present in groundwaters may be trapped in the solid phase of minerals like calcite or gypsum, either by adsorption or by co-precipitation. When a contaminant is incorporated in the bulk rather than simply adsorbed at the surface, it is less available and it can be considered "immobilized" in the environment at least until the host phase dissolution. The aim of this study is to elucidate whether the incorporation of As(III) and As(V) into the bulk of calcite and gypsum, respectively, occurs or not, and to what extent.

Surface-sensitive X-ray Standing Wave (XSW) studies by Cheng et al. [START_REF] Cheng | [END_REF] show that the As atom replaces the C atom in the carbonate molecules of calcite. The geometry of the carbonate group is not preserved, showing a displacement of the As atom of 0.76 Å in the [001] direction. Density Functional Theory (DFT) based simulations have proved that this replacement drives to a similar displacement of the As atom when C atoms are replaced by As atoms in the bulk of calcite (see below). This fact leads us to keep the same hypothesis of arsenite/carbonate replacement in our study of As incorporation into the bulk.

The stability diagram of aqueous solutions of HnAsO n-3

4

shows that the arsenate (AsO 3- 4 ) is the most stable specie under oxidizing conditions. This anion is a tetrahedron with the As atom at the centre, surrounded by four O atoms. The fact that both, arsenate and sulphate groups have the same geometry supports the hypothesis of a possible replacement of sulphate by arsenate groups when gypsum is precipitated in presence of As(V). The charge is compensated by bonding to an extra proton.

Materials and methods

In order to test the possible mechanisms for As immobilization by calcite and gypsum, samples of both minerals were synthesized in the presence of As(III) and As(V), respectively. Calcite precipitation was conducted at pH = 7.5 by addition of CaCl2 and Na2CO3 solutions. Gypsum was precipitated from Na2SO4 and CaCl2 solutions at three different pH values: 4, 7.5 and 9. Arsenic concentrations incorporated in the solids range between 30 mM/kg and 1200 mM/kg for calcite and 100 mM/kg to 1000 mM/kg for gypsum.

Powder samples were analysed by neutron diffraction at the high flux powder diffractometer D20 (ILL). Experiments were carried out using a Cu(200) monochromator which gives a wavelength of λ = 1.30 Å and at ambient conditions of pressure and temperature [START_REF] Román-Ross | [END_REF]. Diffraction patterns were taken for the samples in their container and for the empty cell, in the range of 10 • to 130 • . Also powder diffraction experiments were performed with x-ray at ID11 (ESRF), reproducing the same experimental conditions but using a wavelength of λ = 0.5 Å. Both diffraction data sets were analysed using Fullprof [START_REF] Rodríguez-Carvajal | Collected Abstracts of Powder Diffraction Meeting[END_REF].

Geometrical optimisations of the unit cell and the supercells of pure and As-doped calcite and gypsum were done with the Vienna Ab-initio Simulation Package (VASP) [START_REF] Kresse | Software VASP[END_REF]. The PBE functional and PAW pseudopotentials were used. The goal was to reproduce the expansion of the unit cell produced by the incorporation of As atoms into the structure of both minerals. Unit cells of pure calcite and gypsum obtained from Rietveld refinements were used as starting point for all the models. Geometrical optimizations of single unit cells and of supercells of gypsum were done replacing the sulphate molecules SO 2- 4 by arsenates AsO 3- 4 . Similar simulations were performed with supercells of calcite replacing C by As atoms [START_REF] Román-Ross | [END_REF].

EXAFS data were collected on a diffraction and absorption beamline (GILDA-BM8) at the ESRF of Grenoble and extracted using standard procedures [7]. The theoretical photoelectron paths were generated using the FEFF8 [8] code and the fit performed using the MINUIT library from CERN [START_REF] James | CERN Program Library[END_REF]. A monochromator of Si(311) was used to set the incident energy at the K-edge of As (11867 eV).

Results and discussion

Our diffraction data show an expansion of the unit cell due to As incorporation into calcite crystallites (Fig. 1). By modelling, As concentration in the bulk of the samples can be extrapolated by comparison of the relative volume changes between the experimental and simulated data. Calcite unit cell and two supercells (2x2x1 and 3x2x1) were geometrically optimised replacing one and two units of AsO 3- 3 (150 mM/kg and 290 mM/kg) by CO 2- 3 units. The simulations showed a volume expansion linearly dependent on the replacement of As, as expected by Vegard's law (Fig. 2). This augmentation is due to the lattice expansion along the a axis as As concentration in solids increases. The observed As ion displacement of 0.57 Å over the O base along the (0001) is compatible with Cheng's results [START_REF] Cheng | [END_REF]. The lower value of this displacement (0.57 Å vs. 0.76 Å) can be due to the fact that atoms near the surface are less attached to the solid and can move more freely.

The experimental value of the unit cell volume was interpolated using the linear fit of the simulated volume expansion, giving values of 9, 10 and 16 mM/kg of As in calcite for the three synthesised samples. Simulations of 3x2x1 supercells with one CO 2- 3 unit replaced by one AsO 3- 3 unit were done to check whether the replacement is more likely in sites at the same or at different crystallographic planes. The higher enthalpy of formation (∆H = 211.36 meV) for the calcite structure with two As atoms lying on the same plane shows that replacement is more likely to happen in different planes, leading to a more stable structure. We found from EXAFS data analysis a nearest neighbour distance of ≈ 1.77(3) Å corresponding to the As-O bond distance. This value lies in between the reported one [START_REF] Loewenschuss | [END_REF] for the arsenite molecule (d L As-O = 1.84 Å) and that obtained from the simulations (d sims

As-O = 1.75 Å). The coordination number was kept fixed to its theoretical value (N = 3) in order to reduce the correlation between free parameters in the fitting procedure. This result supports the hypothesis of the incorporation of the As atoms into the C crystallographic sites.

Figure 1 shows the gypsum unit cell volume obtained from combined refinements of neutron and x-ray data. The expansion of the unit cell is proportional to the As concentration in solids and strongly dependent on the pH value: the biggest expansion is found in samples synthesized at pH 9. This result is in good agreement with the hypothesis of replacement of sulphate (SO 2- 4 ) by arsenate groups (AsO 3- 4 ). This replacement is more likely at higher pH values, according to available thermodynamical data regarding speciation of As [START_REF] Pourbaix | Atlas of Electrochemical Equilibria in Aqueous Solutions[END_REF]. The expansion of the unit cell parameters is due to the different lengths for As-O (≈ 1.70 Å) [12] and S-O (≈ 1.50 Å) [11] bonds.

Simulations show an increasing of the unit cell volume proportional to the number of atoms of S replaced by As (Fig. 2). Four single cells (with 940, 1809, 3357 and 4696 mM/kg of As) and four supercells were simu-lated: two 2x1x2 supercells with one and two As atoms, a 2x1x3 and a 3x1x3 with one As atom each, giving As concentrations of 358, 705, 240 and 160 mM/kg, respectively. The simulations allow us to extrapolate the As concentration in the bulk of the samples by comparing the relative volume variations between the experimental and simulated data (Table 1). 1 Extrapolated values of the As concentration in samples of gypsum. The result for pH 4 is due to the near zero relative expansion of the unit cell.

Our results support the hypothesis of As immobilisation by incorporation into the bulk of these minerals. This improves the knowledge on the long term stability of contaminated sludges and it has important consequences for site remediation actions. This work is an example of a direct link between fundamental research and environmental issues. The understanding of the As compounds behaviour in sedimentary environments is essential to estimate and predict possible consequences of forecast or accidental events.

Fig. 1 .

 1 Fig. 1. Experimental unit cell volume of As-doped calcite and gypsum samples in function of As concentration.

Fig. 2 .

 2 Fig. 2. Simulated relative expansion of unit cell volume of gypsum and calcite. The insets show a detail of the interpolation of the experimental data using the simulated expansion.