The snapping out Brownian motion

The Annals of Applied Probability 2016, Vol. 26, No. 3, 1727–1742 DOI: 10.1214/15-AAP1131 © Institute of Mathematical Statistics, 2016 THE SNAPPING OUT BROWNIAN MOTION B Y A NTOINE L EJAY1 Inria Nancy Grand-Est We give a probabilistic representation of a one-dimensional diffusion equation where the solution is discontinuous at 0 with a jump proportional to its flux. This kind of interface condition is usually seen as a semi-permeable barrier. For this, we use a process called here the snapping out Brownian motion, whose properties are studied. As this construction is motivated by applications, for example, in brain imaging or in chemistry, a simulation scheme is also provided. 1. Introduction. Many diffusion phenomena have to deal with interface conditions. Let D be a diffusivity coefficient which is smooth away from a regular surface S, but presents some discontinuity there. In this case, the solution to the diffusion equation (1)   ∂t u(t, x) = 21 ∇ D(x)∇u(t, x) = 0 with u(0, x) = f (x) has to be understood as a weak solution. However, u is smooth away from S and satisfies (2) u(t, x+) = u(t, x−) and D(x+)n+ (x) · ∇u(t, x+) = D(x−)n− (x) · ∇u(t, x−), for x ∈ S, when S is assumed to separate locally Rd into a “+” and a “−” part and where n± is a vector normal to S at x pointing to the “±” side. The second condition is called the continuity of the flux. Now, let us assume that D takes scalar values, and is constant away from a thin layer of width 2ℓ enclosed between two parallel surfaces S + and S − . When the width ℓ of the layer tends to 0, S + and S − merge into a single interface located on a surface S. When the diffusivity D0 decreases to 0 with ℓ and D0 /ℓ → λ > 0, then the solution to (1) converges to a function v satisfying (1) away from S with the interface condition for x ∈ S: (3) ∇v(t, x+) = ∇v(t, x−) and  λ v(t, x+) − v(t, x−) = D(x±)∇v(t, x±). 2 Received January 2013; revised July 2015. 1 Supported by the ANR SIMUDMRI (ANR-10-COSI-SIMUDMRI). MSC2010 subject classifications. Primary 60J60; secondary 60G20, 60J35, 60J55. Key words and phrases. Interface condition, elastic Brownian motion, semi-permeable barrier, thin layer, piecing out a Markov process. 1727 1728 A. LEJAY F IG . 1. The thin layer problem. The solution has a continuous flux on S but is discontinuous on S (see, e.g., [33], Chapter 13). A heuristic explanation is given Figure 1. If D is smooth on Rd , it is well known that   u(t, x) = Ex f (Xt ) , (4) where X is the diffusion process generated by 21 ∇(D∇) which is solution under Px to the stochastic differential equation (SDE) (5) Xt = x +  t 0 σ (Xs ) dBt +  t  1 d Di,· 0 2 i=1 ∂xi (Xs ) ds with σ σ T = D for a Brownian motion B. When D presents some discontinuities, (5) has no longer a meaning. However, a Feller processes (X, (Ft )t≥0 , (Px )x∈R ) is associated to 21 ∇(D∇·) for which (4) holds. In particular, the marginal distributions Xt have a density p(t, x, ·) under Px , where p(t, x, y) is the fundamental solution to (1) (see, e.g., [36]). Let us now assume that the dimension of the space is equal to 1 and that D is discontinuous at some separated points {xi } with left and right limit there, and smooth elsewhere. The process X is solution to a SDE with local time. The Itô– Tanaka formula is the key tool to manipulate it, and several simulation algorithms have been proposed (see the references in [25], e.g.). The process called the Skew Brownian motion is the main tool for this construction [22, 24]. Coming back to the thin layer problem, we assume that D is constant and equal to D1 on (−∞, −ℓ) and (ℓ, ∞), and to D0 on (−ℓ, ℓ). The associated stochastic process is solution to Xt = x +  t 0 D(Xs ) dBs + D0 − D1 −ℓ D1 − D0 ℓ Lt (X) + L (X), D1 + D0 D1 + D0 t where L±ℓ t (X) is the local time of X at ±ℓ [24]. Letting D0 /ℓ converging to 2κ with ℓ → 0, one may expect that X converges in distribution to a stochastic process Y such that the solution to (1) with the interface condition (3) is given by v(t, x) = Ex [f (Yt )]. THE SNAPPING OUT BROWNIAN MOTION 1729 The article then aims at constructing and giving several properties related to the process Y which we call a snapping out Brownian motion (SNOB). This process is Feller on G = (−∞, 0−] ∪ [0+, +∞) but not on R. The intervals in the definition of G are disjoint so that 0 corresponds either to 0+ or 0− seen as distinct points. The behavior of this process is the following: Assume that its starting point is x ≥ 0. It behaves as a positively reflected Brownian motion until its local time is greater than an independent exponential random variable of parameter 2κ. Then its decides its sign with probability 1/2 and starts afresh as a new reflected Brownian motion, until its local time is greater than a new exponential random variable, and so on. Using the properties of the exponential random variable, it is equivalent to assert that the particle changes its sign when its local time is greater than an exponential random variable with parameter κ, and behaves like a positively or negatively reflected Brownian motion between these switching times. Its name is justified by the following fact: As the time at which the particle possibly changes it signs is the same as for the elastic Brownian motion [10, 15, 18, 19] (also called the partially reflected Brownian motion), it could also be seen as some elastic Brownian motion which is reborn once killed. The elastic Brownian motion, also called a partially reflected Brownian motion, is associated to the Robin boundary condition and has then many applications [8, 15, 35]. This process is the “basic brick” for constructing the SNOB. The behavior of the SNOB justifies also the old heuristic that the interface condition (3) corresponds to a semi-permeable barrier, which arises, for example, in diffusion magnetic resonance imaging [11] or in chemistry [1, 8]. The interface condition (3) is different from (2), to which is associated a Skew Brownian motion and where the particle crosses the interface when it reaches it, and which corresponds to a permeable barrier (see references in [22, 25]). Here, we work under the condition of a single interface at 0. In short time, it is sufficient to describe the behavior of the process even in a more complex media, since other interface or boundary conditions far enough have “exponentially small” influence on the distribution of the process. This is sufficient for simulation purposes, where particles positions are represented by the stochastic process and move according to its dynamic. Using similar computations, one may generalize our work to the case where D(x) = D + if x ≥ 0, D − if x ≤ 0 and an interface condition ∇u(t, 0+) = β∇u(t, 0−) and λu(t, 0+) − μu(t, x−) = ∇u(t, x+) with λ, μ > 0. Diffusions on graphs specified by a condition at each vertex could also be considered, which could be of interest in several applications. This process has been described without proof by Bobrowski in [6], which have studied its limit behavior when the diffusion coefficients increase. Although the SNOB may be seen as a diffusion on a graph, it is not a diffusion on a metric graph, where the edges are joined by vertices. Such diffusions have 1730 A. LEJAY been classified by Freidlin and Wentzell in [12, 13]. The conditions that are required at the vertices of the graphs are some extension of the possible boundary conditions for a Markov process studied by Feller [10]. See also [21], for example, for the related problem of pasting diffusions.2 Our interface condition does not fall in these categories. Our process is best thought as a kind of random evolution process which switches back and forth randomly among a collection of processes (see, e.g., [16, 34]). Outline. In Sections 2 and 3, we present quickly the main results related to the elastic Brownian motions and the piecing out procedure. The SNOB is constructed in Section 4 through its resolvent. In Section 5, we show the relationship between the SNOB and the thin layer problem. Finally, in Section 6, we show how to simulate this process. 2. Elastic Brownian motion. Let (Rt )t≥0 a reflected Brownian motion, and denote by (Lt )t≥0 its symmetric local time at 0. We add a cemetery point † to R+ . For a constant κ > 0, we consider an exponential random variable ξ with parameter κ independent from B. Set Zt = Rt , †, if Lt ≤ ξ , if Lt > ξ . Thanks to the properties of the local time, this process, called the elastic Brownian motion (EBM), is still a strong Markov process. Its semi-group is   Pte f (x) = Ex exp(−κLt )f (Xt ) for f in the set C0 (R+ , R) of continuous functions that vanishes at infinity. Closed form expressions of the density transition function are given in [14, 35]. Let k be the time at which the EBM is killed, which means k = inf{t > 0|Lt ≥ ξ }. This is a stopping time. Since the local time increases only on the closure of Z = {t > 0|Xt = 0}, it holds that Zk = 0 almost surely. Using standard computations in the inverse of the local time of the Brownian motion, √   κ exp(− 2αx). (6) ψ(x, α) = Ex exp(−αk) = √ 2α + κ Using the Itô formula, it is easily shown that u(t, x) = Pte f (x) is solution to the heat equation with Robin (or third kind) boundary condition [3, 15, 31] ⎧ ∂u(t, x) 1 ⎪ ⎪ = △u(t, x), ⎨ ∂t 2 ⎪ ∂u(t, 0) ⎪ ⎩ = κu(t, 0). ∂x on (0, +∞)2 , 2 The article [32] defines a notion of semipermeable membrane which is different from ours, where the solution is continuous with a discontinuous gradient. 1731 THE SNAPPING OUT BROWNIAN MOTION For a Markov process X, let us recall that its resolvent (Gα )α>0 is a family of  operators defined by Gα f (x) = Ex [ 0+∞ e−αs f (Xs ) ds] for any f ∈ C0 and any α > 0. It has a density gα when Gα f (x) = gα (x, y)f (y) dy. Using standard computations on the Green functions, the density gαe (x, y) of the resolvent of the EBM is for x, y ≥ 0, ⎧√ √ ⎪ 2α − κ −√2α(y+x) ⎪ ⎪ √ ⎪ + e− 2α(x−y) , for y ∈ [0, x], e ⎨ 1 2α + κ √ gαe (x, y) = √ √ 2α − κ −√2α(x+y) 2α ⎪ ⎪ ⎪ 2α(x−y) ⎪ e +√ , for y ≥ x. ⎩e 2α + κ We extend the EBM to a process on G by symmetry, so that its resolvent becomes (7) Geα f (x) := Ex  k 0  e−αs f (Xs ) ds =  +∞ 0     gαe |x|, y f sgn(x)y dy for x ∈ G. This process evolves either on R− or R+ but never crosses 0 and is naturally identified with a process on G. 3. Piecing out Markov processes. The procedure of piecing out is a way to construct a Markov process from a killed one. We present in this section a result due to Ikeda, Nagasawa and Watanabe [17] (similar considerations are given in [29]). On a probability space ( , F , P) and a state space S, let ((Xt )t≥0 , (Ft )t≥0 , (Px )x∈S ) be a right continuous strong Markov process living in the extended state space S† = S ∪ {†} with a death point †. The lifetime of X is denoted by k. The shift operator associated to X is denoted by (θt )t≥0 . We also consider a family μ defined on × S† such that μ(ω, ·) is a probability measure on S† and for any fixed Borel subset A, μ(·, A) is σ (Xt , t ≥ 0)measurable. We assume additionally that μ(ω, dy) = δ† (dy) when k(ω) = 0 and   Px μ(ω, dy) = μ(θt(ω) ω, dy), t(ω) < k(ω) = Px [t < k] for any stopping time t. The family μ, called an instantaneous distribution, describes the way the process is reborn once killed. Let  be the product of an infinite, countable, number of copies of × S† . We define X on  by ) = Xt (ω ⎧ xt (ω1 ), ⎪ ⎪ ⎪ ⎪ ⎪ y1 , ⎪ ⎪ ⎪ ⎨ 2 ), xt−k(ω1 ) (ω ⎪ y2 , ⎪ ⎪ ⎪ ⎪ ⎪··· ⎪ ⎪ ⎩ †,   if t ∈ 0, k(ω1 ) , if t = k(ω1 ),   if t ∈ k(ω1 ), k(ω1 ) + k(ω2 ) , if t = k(ω2 ), if t ≥ k(ω1 ) + · · · + kN (ωN ) 1732 A. LEJAY  = (ω1 , y1 , ω2 , y2 , . . .) ∈  and N = inf{k ≥ 0; k(ωk ) = 0}. with ω We consider the probability measure    Px dω1 , dx 1 , . . . , dωn , dx n  1  1 1         = Px dω μ ω , dx Px 1 dω2 μ ω1 , dx 2 · · · Px n dω2 μ ωn , dx n . Px , when the path X(ω) is killed, we let it reborn by placing Under this measure  it at the point x1 with probability μ(ω, dx1 ) and then start again. We left the technical details about the construction of the probability space and the filtration and presents the main result on piecing out Markov process. T HEOREM 1 ([17]). Using the above defined notation, there exists a probabil,  t )t≥0 on which (X, (B t )t≥0 , ( ity space (  , B P) and a filtration (B Px )x∈S† ) is a † strong Markov process on S with P† [Xt = †, ∀t ≥ 0] = 1. 4. The snapping out Brownian motion. D EFINITION 1. A snapping out Brownian motion (SNOB) X is a strong Markov stochastic process living on G constructed by making EBM reborn on 0+ or 0− with probability 1/2 using the piecing-out procedure. The sign of X changes with probability 1/2 when its local time Lt at 0 is greater than uk with u0 = 0, uk − uk−1 ∼ exp(κ) is independent from (ui )i≤k−1 . From the properties of the exponential and binomial distributions, the sign of X changes when its local time is greater than sk with s0 = 0, sk − sk−1 ∼ exp(κ/2) is independent from (si )i≤k−1 . It is also immediate that |X| is a reflected Brownian motion, where | · | is the canonical projection of G onto [0, +∞). P ROPOSITION 1.  The resolvent family (Gα )α>0 of the SNOB is solution to  1 α − △ Gα f (x) = f (x) for x ∈ G 2  κ with ∇Gα f (0+) = ∇Gα f (0−) and Gα f (0+) − Gα f (0−) = ∇Gα f (0) 2 for any bounded, continuous function f on G that vanishes at infinity. This proposition identifies the infinitesimal generator of the process X. The points 0+ and 0− are then interpreted as the sides of a semi-permeable barrier. P ROOF OF P ROPOSITION 1. From this very construction and the strong Markov property, for any continuous function f on G which vanishes at infinity, (8) Gα f (x) = Geα f (x) +  ψ(|x|, α)  Gα f (0+) + Gα f (0−) , 2 1733 THE SNAPPING OUT BROWNIAN MOTION where Geα is defined by (7). Using x = 0+ and x = 0− in (8) and summing the two resulting equations leads to √ Gα f (x) = Geα f (x) + (9) κe− 2α|x| √ β(f ) 2 2α with β(f ) = Geα f (0+) + Geα f (0−). Then √ κ (10) Gα f (x) + Gα f (−x) = Geα f (x) + Geα f (−x) + √ e− 2α|x| β(f ), 2α (11) Gα f (x) − Gα f (−x) = Geα f (x) − Geα f (−x). Derivating (10) and setting x = 0+, since ∇Geα f (0±) = ±κGeα f (0±), ∇Gα f (0+) − ∇Gα f (0−) = 0. Derivating (11), 2∇Gα f (0±) = ∇Gα f (0+) + ∇Gα f (0−) = ∇Geα f (0+) + ∇Geα f (0−)     = κ Geα f (0+) − Geα f (0−) = κ Gα f (0+) − Gα f (0−) . In addition, it is easily seen that (α − 21 △)Gα f = f since ψ(x, α) is solution to (α − 12 △)ψ(x, α) = 0. The resolvent is then identified.  P ROPOSITION 2. resentation: The semi-group (Pt )t≥0 of the SNOB has the following rep- Pt f (x) = Ex (12)  + Ex for a Brownian motion B. P ROOF.    1 + e−κLt f sgn(x)|Bt | 2       1 − e−κLt f − sgn(x)|Bt | 2 Let us decompose a function f as its even and odd parts: fˆ(x) =  1 2 f (x) + f (−x) and fˇ(x) =  1 2 f (x) − f (−x) . Then Geα fˆ(−x) = Geα fˆ(x) and Geα fˇ(−x) = −Geα fˇ(x), so that β(fˇ) = 0 for β defined by (9). Thus Gα fˇ(x) = Geα fˇ(x). In addition, since fˆ(|x|) = fˆ(x) and the SNOB has the same distribution as the reflected Brownian motion |B|, Gα fˆ(x) = Grα fˆ(x) := Ex  +∞ 0    e−αs fˆ |Bs | ds . 1734 A. LEJAY This gives an alternative representation for the resolvent of the SNOB: Gα f (x) = Grα fˆ(x) + Geα fˇ(x). Inverting the resolvent to recover the semi-group (Pt )t≥0 ,       Pt f (x) = Ptr fˆ(x) + Pte fˇ(x) = Ex fˆ |Bt | + Ex exp(−κLt )fˇ sgn(x)|Bt | . This expression could be arranged as (12).  5. The thin layer problem. We now fix ε > 0 and we consider the process Xε generated by (see, e.g., [36] for general considerations on this process) Lε :=  ∂ 1 ∂ a ε (x) 2 ∂x ∂x  with a ε (x) := when x ∈ / [−ε, ε], when x ∈ [−ε, ε] 1, κε, whose domain Dom(Lε ) = {f ∈ L2 (R)|Lε f ∈ L2 (R)} is a subset of the Sobolev space H1 (R) [hence, any function in Dom(Lε ) is identified with a continuous function], where L2 (R) is the set of square integrable functions on R with scalar product f, g = R f (x)g(x) dx. Let us set [h](x) := h(x−) − h(x+) and (13) ⎫   f, f ′′ ∈ L2 (R), ⎬   . D ε := f ∈ C (−∞, −ε) ∪ (−ε, ε) ∪ (ε, ∞)  [f ](±ε) = 0, ⎩  a ε ∇f (±ε) = 0 ⎭ ⎧ ⎨ 2 For k ≥ 0, we write Cck (R) the set of functions with compact support and continuous derivatives up to order k. With an integration by parts, for f ∈ D ε and g ∈ Cc2 (R),   (α − L)f, g = α f, g +    R a ε (x)∇f (x)∇g(x) dx   + a ε ∇f (−ε)g(−ε) − a ε ∇f (ε)g(ε). Using this formula and the regularity of the solution to (α − L)f = g when g ∈ C ∞ (I, R) with −ε, ε ∈ / I , we easily get that D ε contains (α − Lε )−1 (Cc∞ (R)) and is then dense in Dom(Lε ) for the operator norm ( f, f  + Lf, Lf )1/2 . A fundamental solution may be associated to Lε , as well as a resolvent density ε gα , which we will compute explicitly. This operator is self-adjoint with respect to ·, ·, so that its resolvent density satisfies gαε (x, y) = gαε (y, x). This process is a Feller process, and is a strong solution to the SDE with local time  t      ε 1 − κε , a ε Xsε dBs + ηε Lεt X ε − ηε L−ε X with ηε = Xtε = x + t 1 + κε 0 where B is a Brownian motion and Lxt (X ε ) is the symmetric local time at x of Xε (see, e.g., [24], and [4, 28] among others for general results on SDEs with local time). In [10], Section 11, the elastic Brownian motion is constructed as the limit of a process which either jumps at ε or is killed with probability κε when it arrives at 0. 1735 THE SNAPPING OUT BROWNIAN MOTION Using the piecing out procedure, we construct a strong Markov process Z ε by considering the process X ε which is instantaneously replaced at −ε or ε with probability 1/2 when it reaches 0, and then behaving again as Xε until it reaches 0, and so on. This process Z ε could be identified as a process living in G by defining P0+ as Pε and P0− as P−ε , since the process is instantaneously killed when at 0. T HEOREM 2. The process Z ε with Z0ε = x converges in distribution to the SNOB starting from x in the Skorohod topology. The proof relies on the next two results. P ROPOSITION 3. Let gαε be the resolvent density of X ε . Then gαε (x, y) converges to g(x, y) for any x, y = 0 and any α > 0 as ε → 0. R EMARK 1. This result follows from classical results in deterministic homogenization theory (see, e.g., [33]) where the convergence holds in Sobolev spaces. Here, we consider a direct computational proof for the convergence of the Green kernel, which we use later. P ROOF OF P ROPOSITION 3. We assume that x > 0 and we set μ := some α > 0. The resolvent density gαε of X ε has the form, for x > ε, gαε (x, y) = ⎧ Cε (x)e−μy , ⎪ ⎪ ⎪ ⎨ −μy + Bε Aε (x)e √ (x)eμy , √ ⎪ Hε (x)eμy/ κε + Eε (x)e−μy/ κε , ⎪ ⎪ ⎩ μy Fε (x)e , √ 2α for for y > x, for y ∈ [ε, x], for y ∈ [−ε, ε], for y < −ε. By this, we mean that for any bounded, measurable function f , Ex The kernel gεα  +∞ 0 e −αt    ε f Xs ds = gαε (x, y)f (y) dy. R satisfies the conditions gαε (x, ε+) = gαε (x, ε−), gαε (x, ε−) = gαε (x, ε+), ∇y gαε (x, −ε−) = κε∇y gαε (x, −ε+), κε∇y gαε (x, ε−) = ∇y gαε (x, ε+), ∇y gαε (x, x+) − ∇y gαε (x, x−) = 2. √ With μ = 2α, the coefficients Aε , Bε , Cε , Hε and Fε are then expressed with the help of √ √ √ √ √   Gε := 2e4με/ κε κε + e4με/ κε κε + 2 κε − κε + e4με/ κε − 1 μ. 1736 A. LEJAY √ Since ε → 0, Gε = 4 κε(1 + μ κ + O(κε))μ. After tedious computations, √ √   Aε (x) = − e4με/ κε κε − κε − e4με/ κε + 1 eμ(2ε−x) /Gε −→ A0 (x) ε→0 := − e−μx κ +μ , e−μx , μ √   Cε (x) = −2 sinh(2ε/ κε) de−μx+2με − e−μx+2με + deμx + e−μx √ √ √ √ × e2με/ κε /Gε + 4 κεeμx cosh(2ε/ κε)e2με/ κε /Gε Bε (x) = B0 (x) := − −→ C0 (x) := ε→0 Hε (x) = −2eμ(3ε+ε := − κeμx , μ(κ + μ) √ √ √ κε−x κε)/ κε κe−μx 2μ(κ + μ) Eε (x) = −2eμ(ε+ε √ Fε (x) = −4 κεe √ (1 + √ √ κε) κε/Gε −→ H0 (x) ε→0 , √ √ κε−x κε)/ κε (1 − ε→0 √ μ(2ε+2ε κε−x κε)/ κε √ √ √ √ κε) κε/Gε −→ H0 (x), /Gε −→ F0 (x) := ε→0 = −C0 (−x). Let gα be the function −κe−μx μ(κ + μ) ⎧ −μy , ⎪ ⎨ C0 (x)e if y > x, A0 (x)e−μy + B0 (x)eμy , if y ∈ [0, x], ⎪ ⎩ μy F0 (x)e , if y < 0. A similar work may be performed for x < 0. Thus, we easily obtain that gαε (x, y) −→ε→0 gα (x, y) converges to gα and that gα is the density resolvent of the SNOB by checking it satisfies the appropriate conditions at the interface.  gα (x, y) := P ROPOSITION 4. Let hε0 be the first hitting time of 0 for Xε . Under Px , hε0 converges in distribution to a random variable k distributed as the lifetime of the EBM of parameter κ. P ROOF. defined by As in [9, 24], we introduce ε (x) as the piecewise linear function √ dε 1/ κε, (x) = dx 1, if x ∈ [−ε, ε], otherwise. 1737 THE SNAPPING OUT BROWNIAN MOTION Set Y ε = ε (X ε ) so that Y ε is solution to the SDE [9, 24]  y   −y  Ytε = ε (x) + Bt + θ ε Lt ε Y ε − θ ε Lt ε Y ε √  ε 1 − κε √ and y ε := ε (ε) = . with θ ε = 1 + κε κ The infinitesimal generator of Y ε is Lε := 21 △ whose domain contains as a dense subset [it is similar to the discussion on D ε in (13)] ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎫   f, f ′′ ∈ L2 (R), ⎪  ⎪   ⎪  [f ] ±y ε = 0, ⎪ ⎪  ⎪ ⎪     ⎬ ε ′ ε   y − f 1 − θ  2     f ∈ C (−∞, −yε ) ∪ (−yε , yε ) ∪ (yε , ∞)  .  ⎪ = 1 + θ ε f ′ y ε + , ⎪ ⎪ ⎪    ⎪ ⎪ ⎪ ⎪  1 + θ ε f ′ −y ε − ⎪ ⎪ ⎪  ⎪ ⎪   ′ ε  ⎪ ⎩ ⎭ ε  = 1 − θ f −y + , From now, we assume for the sake of simplicity that x > 0. The hitting time hε0 is also the first hitting time of zero by Y ε . Since by symmetry ψ(−x, α) = ψ(x, α) for any x ≥ 0, we consider only that x ≥ 0. Since the Feynman–Kac formula is valid for the process Y ε , ψ ε (x, α) := ε Ex [e−αh0 ] is solution to ⎧1 △ψ ε (x, α) = αψ ε (x, α), ⎪ ⎪ ⎪2ε ⎨ ψ (0, α) = 1,     ε y ε −, α = ψ ε y ε +, α , ⎪ ψ ⎪ ⎪     ⎩ for x = y ε ,    1 − θ ε ∇x ψ ε y ε −, α = 1 + θ ε ∇x ψ ε y ε +, α . Hence, ψ ε (x, α) is sought as  √ γ ε exp(− 2αx), √ √ ψ (x, α) = cos( 2αx) + β ε sin( 2αx), ε After some computations, √ √ √ − cos( 2αy ε ) + κε sin( 2αy ε ) ε √ √ β = √ sin( 2αy ε ) + κε cos( 2αy ε ) and if x > y ε ,   if x ∈ 0, y ε . √ √ − κ √ . εβ ∼ ε→0 κ + 2α ε Besides, γε =e √ 2αy ε √ Hence, for any x > 0, (14)  κε β ε cos √  2αy ε − β ε sin ψ ε (x, α) −→ ψ(x, α) := with ψ defined by (6). ε→0 √ 2αy ε  κ √ . ε→0 κ + 2α ∼ √ κ √ e− 2αx κ + 2α 1738 A. LEJAY This proves that under Px , hε0 converges to a random variable k whose Laplace transform is ψ(x, α) under Px . This random variable k is then the lifetime of an EBM.  P ROOF OF T HEOREM 2. Using the properties of the resolvent, for α > 0 and a bounded, measurable function f , Gεα f (x) := Ex  +∞ e 0 −αt    f Xsε ds ε = Rαε f (x) + Ex e−αh0 with Rαε f (x) := Ex Since ψ ε (x, α) = ψ ε (−x, α), Gεα f (x) = Rεα f (x) +  hε 0 0 e   1 ε Gα f (ε) + Gεα f (−ε) 2 −αt   ε f Xs ds . ψ ε (x, α) Rεα f (ε) + Rεα f (−ε) . 1 − ψ ε (ε, α) 2 For the sake of simplicity, we assume that x > 0. Using the symmetry properties of Lε , Rαε f (x) = But  +∞  0  gαε (x, y) − gαε (x, −y) f (y) dy. gαε (x, y) − gαε (x, −y) −→ gα (x, y) − gα (x, −y) = gαe (x, y), ε→0 gαe (x, y) is the resolvent density of the EBM. Thus, Rεα f (x) −→ε→0 where e Gα f (x) for any x > 0. It is also easily obtained that Rεα f (ε) −→ Geα f (0+) ε→0 and Rεα f (−ε) −→ Geα f (0−). ε→0 Using (9) and (14), Gεα f (x) −→ε→0 Gα f (x). The Trotter–Kato theorem (see, e.g., [20], Theorem IX.2.16, page 504) and the Markov property imply the convergence in finite-dimensional distributions of Z ε to X under Px for x ≥ 0. By symmetry, this could be extended to x ≤ 0. The only remaining point of the tightness. When away from [−ε, ε], X ε behaves like a Brownian motion. Hence, for 0 ≤ s ≤ t ≤ T , let us set f(s, t) := inf{u > s; |Xuε | = ε} with possibly f(s, t) = +∞ and l(s, t) := sup{u < t; |Xuε | = ε} with possibly l(s, t) = −∞. If f(s, t) ≥ t and l(s, t) ≤ s, then for δ < 1/2, there exists an integrable random variable C(ω) such that |Xtε (ω) − Xsε (ω)| ≤ C(ω)(t − s)δ for any 0 ≤ s ≤ t ≤ T . 1739 THE SNAPPING OUT BROWNIAN MOTION If f(s, t) ≤ t and l(s, t) ≤ s, then    ε X − X ε  ≤ X ε t s f(s,t)    ε  ≤ C(t − s)β + 2ε − Xsε  + Xtε − Xf(s,t) since Xtε belongs to [−ε, ε]. A similar analysis could be carried for the other cases, which means that for some integrable random variable C, |t−s|<δ (Z ε )ε>0   sup Xtε − Xsε  ≤ Cδ β + 2ε. This proves that is tight is the space D([0, T ]; R) of discontinuous functions with the Skorohod topology (see, e.g., [5]) and then on D([0, T ]; G). Hence, we easily deduce the convergence of Z ε to the SNOB in D([0, T ]; G).  6. Simulation of the SNOB. It is easy to simulate a discretized process X in the same way it is easy to simulate the Brownian motion. Following Proposition 2, we draw a random variate with density p(δt, x, ·) when x is close enough to 0. For this, we use a Brownian bridge technique to check if the process reaches 0± before δt (see, e.g., [2] and [25], Section B.2, for an example of application and further references). Thisinvolve the inverse Gaussian distribution 2 ). Random variates with IG (λ, μ) whose density is rμ,λ (x) = 2πλx 3 exp( −λ(x−μ) 2μ2 x IG distribution could be simulated by the methods proposed in [7], page 148 and [30]. We simulate the local time using the following representation under P0 [26, 27]:  dist  0 Lt (B), |Bt | = (l, l − H ) where l := 1 2 H  + V + H2  with H ∼ N (0, t) and V ∼ exp(1/2t) independent from H . The generic algorithm to simulate the process at time δt when at point x at time 0 is the following: √ 1. Set y := x√+ δtG with G a random variate whose distribution is N (0, 1). 2. If |x| ≥ 4 δt , then return y (here, we neglect the exponentially small probability that the process crosses 0 between the times 0 and δt). 3. If xy > 0, then decide with probability exp(−2|xy|/δt) if the path X has crossed 0. • If no crossing occurs, then return y. • If a crossing occurs, draw g ∼ IG (|x|/|y|, x 2 /2δt), so that z := δtg/(1 + g) is a realization of the first hitting time of 0 for a Brownian bridge with B0 = x and Bδt = y. Then go the step 5. 4. If xy < 0, then draw g ∼ IG (−|x|/|y|, x 2 /2δt) and set z := δtg/(1 + g), the first time the Brownian bridge reaches 0. Go to step 5. 5. Set r := δt − z. For two √ independent random variates H ∼ N (0, r) and V ∼ exp(1/2r), set l := (H + V + H 2 )/2. 1740 A. LEJAY 6. For U ∼ U (0, 1) independent from V and H , set s := sgn(x) if exp(−κl) ≥ 2U − 1. Otherwise, set s := − sgn(x). 7. Return s(l − H ). An application to the estimation of a macroscopic estimation parameter in the context of a simplified problem related to brain imaging may be found in [23]. The results are satisfactory, unless κ is too small due to a problem of rare event simulation. Acknowledgements. The author is indebted to Jing-Rebecca Li and Denis Grebenkov for having proposed this research and interesting discussions about it. The author also wishes to thank warmly his wife, Claire Nivlet, for having suggested the name of the process. REFERENCES [1] A NDREWS , S. S. (2009). Accurate particle-based simulation of adsorption, desorption and partial transmission. Phys. Biol. 6 046015. [2] BALDI , P. (1995). Exact asymptotics for the probability of exit from a domain and applications to simulation. Ann. Probab. 23 1644–1670. MR1379162 [3] BASS , R. F., B URDZY, K. and C HEN , Z.-Q. (2008). On the Robin problem in fractal domains. Proc. Lond. Math. Soc. (3) 96 273–311. MR2396121 [4] BASS , R. F. and C HEN , Z.-Q. (2003). Brownian motion with singular drift. Ann. Probab. 31 791–817. MR1964949 [5] B ILLINGSLEY, P. (1999). Convergence of Probability Measures, 2nd ed. Wiley, New York. MR1700749 [6] B OBROWSKI , A. (2012). From diffusions on graphs to Markov chains via asymptotic state lumping. Ann. Henri Poincaré 13 1501–1510. MR2966471 [7] D EVROYE , L. (1986). Nonuniform Random Variate Generation. Springer, New York. MR0836973 [8] E RBAN , R. and C HAPMAN , S. J. (2007). Reactive boundary conditons for stochastic simulation of reaction–diffusion processes. Phys. Biol. 4 16–28. [9] É TORÉ , P. (2006). On random walk simulation of one-dimensional diffusion processes with discontinuous coefficients. Electron. J. Probab. 11 249–275 (electronic). MR2217816 [10] F ELLER , W. (1954). Diffusion processes in one dimension. Trans. Amer. Math. Soc. 77 1–31. MR0063607 [11] F IEREMANS , E., N OVIKOV, D. S., J ENSEN , J. H. and H ELPERN , J. A. (2010). Monte Carlo study of a two-compartment exchange model of diffusion. NMR Biomed. 23 711–724. [12] F REIDLIN , M. I. and W ENTZELL , A. D. (1993). Diffusion processes on graphs and the averaging principle. Ann. Probab. 21 2215–2245. MR1245308 [13] F REIDLIN , M. I. and W ENTZELL , A. D. (1994). Random perturbations of Hamiltonian systems. Mem. Amer. Math. Soc. 109 viii+82. MR1201269 [14] G ALLAVOTTI , G. and M C K EAN , H. P. (1972). Boundary conditions for the heat equation in a several-dimensional region. Nagoya Math. J. 47 1–14. MR0317658 THE SNAPPING OUT BROWNIAN MOTION 1741 [15] G REBENKOV, D. S. (2006). Partially reflected Brownian motion: A stochastic approach to transport phenomena. In Focus on Probability Theory (L. R. Velle, ed.) 135–169. Nova Sci. Publ., New York. MR2553673 [16] G RIEGO , R. J. and M ONCAYO , A. (1970). Random evolutions and piecing out of Markov processes. Bol. Soc. Mat. Mexicana (2) 15 22–29. MR0365723 [17] I KEDA , N., NAGASAWA , M. and WATANABE , S. (1966). A construction of Markov processes by piecing out. Proc. Japan Acad. 42 370–375. MR0202197 [18] I TÔ , K. and M C K EAN , H. P. J R . (1996). Diffusion Processes and Their Sample Paths. Springer, Berlin. [19] K ARLIN , S. and TAYLOR , H. M. (1981). A Second Course in Stochastic Processes. Academic Press, New York. MR0611513 [20] K ATO , T. (1995). Perturbation Theory for Linear Operators. Springer, Berlin. MR1335452 [21] KOPYTKO , B. I. and P ORTENKO , M. I. (2009). The problem of pasting together two diffusion processes and classical potentials. Theory Stoch. Process. 15 126–139. MR2598532 [22] L EJAY, A. (2006). On the constructions of the skew Brownian motion. Probab. Surv. 3 413– 466. MR2280299 [23] L EJAY, A. (2015). Estimation of the mean residence time in cells surrounded by semipermeable membranes by a Monte Carlo method. Research report No. Inria, RR-8709. [24] L EJAY, A. and M ARTINEZ , M. (2006). A scheme for simulating one-dimensional diffusion processes with discontinuous coefficients. Ann. Appl. Probab. 16 107–139. MR2209338 [25] L EJAY, A. and P ICHOT, G. (2012). Simulating diffusion processes in discontinuous media: A numerical scheme with constant time steps. J. Comput. Phys. 231 7299–7314. MR2969713 [26] L ÉPINGLE , D. (1993). Un schéma d’Euler pour équations différentielles stochastiques réfléchies. C. R. Acad. Sci. Paris Sér. I Math. 316 601–605. MR1212213 [27] L ÉPINGLE , D. (1995). Euler scheme for reflected stochastic differential equations. Math. Comput. Simulation 38 119–126. MR1341164 [28] L E G ALL , J.-F. (1984). One-dimensional stochastic differential equations involving the local times of the unknown process. In Stochastic Analysis and Applications (Swansea, 1983). Lecture Notes in Math. 1095 51–82. Springer, Berlin. MR0777514 [29] M EYER , P. A. (1975). Renaissance, recollements, mélanges, ralentissement de processus de Markov. Ann. Inst. Fourier (Grenoble) 25 465–497. MR0415784 [30] M ICHAEL , J. R., S HUCANY, W. R. and H AAS , R. W. (1976). Generating random variates using transformations with multiple roots. Amer. Statist. 30 88–90. [31] PAPANICOLAOU , V. G. (1990). The probabilistic solution of the third boundary value problem for second order elliptic equations. Probab. Theory Related Fields 87 27–77. MR1076956 [32] P ORTENKO , N. I. (2000). A probabilistic representation for the solution to one problem of mathematical physics. Ukrainian Math. J. 52 1457–1469. [33] S ÁNCHEZ -PALENCIA , E. (1980). Non-Homogeneous Media and Vibration Theory. Lecture Notes in Physics 127. Springer, Berlin. [34] S IEGRIST, K. (1981). Random evolution processes with feedback. Trans. Amer. Math. Soc. 265 375–392. MR0610955 [35] S INGER , A., S CHUSS , Z., O SIPOV, A. and H OLCMAN , D. (2007/08). Partially reflected diffusion. SIAM J. Appl. Math. 68 844–868. MR2375298 1742 A. LEJAY [36] S TROOCK , D. W. (1988). Diffusion semigroups corresponding to uniformly elliptic divergence form operators. In Séminaire de Probabilités, XXII. Lecture Notes in Math. 1321 316–347. Springer, Berlin. MR0960535 I NSTITUT É LIE C ARTAN DE L ORRAINE UMR 7502 U NIVERSITÉ DE L ORRAINE VANDŒUVRE - LÈS -NANCY, F-54500 F RANCE AND I NSTITUT É LIE C ARTAN DE L ORRAINE UMR 7502 CNRS VANDŒUVRE - LÈS -NANCY, F-54500 F RANCE AND TOSCA I NRIA V ILLERS - LÈS -NANCY, F-54600 F RANCE E- MAIL : [email protected]