First Exit Times of Harmonically Trapped Particles: A Didactic Review
aa r X i v : . [ m a t h - ph ] N ov First Exit Times of Harmonically Trapped Particles:A Didactic Review
Denis S. Grebenkov
E-mail: [email protected] Laboratoire de Physique de la Mati`ere Condens´ee (UMR 7643),CNRS – Ecole Polytechnique, 91128 Palaiseau, France St. Petersburg National Research University of Information Technologies,Mechanics and Optics, 197101 St. Petersburg, Russia
Abstract.
We revise the classical problem of characterizing first exit times of aharmonically trapped particle whose motion is described by one- or multi-dimensionalOrnstein-Uhlenbeck process. We start by recalling the main derivation steps of apropagator using Langevin and Fokker-Planck equations. The mean exit time, themoment-generating function, and the survival probability are then expressed throughconfluent hypergeometric functions and thoroughly analyzed. We also present arapidly converging series representation of confluent hypergeometric functions thatis particularly well suited for numerical computation of eigenvalues and eigenfunctionsof the governing Fokker-Planck operator. We discuss several applications of first exittimes such as detection of time intervals during which motor proteins exert a constantforce onto a tracer in optical tweezers single-particle tracking experiments; adhesionbond dissociation under mechanical stress; characterization of active periods of trendfollowing and mean-reverting strategies in algorithmic trading on stock markets;relation to the distribution of first crossing times of a moving boundary by Brownianmotion. Some extensions are described, including diffusion under quadratic double-wellpotential and anomalous diffusion.PACS numbers: 02.50.Ey, 05.10.Gg, 05.40.-a, 02.30.Gp
Keywords : first exit time; first passage time; survival probability; harmonic potential;double-well potential; Ornstein-Uhlenbeck process; confluent hypergeometric function;Kummer and Tricomi functions; parabolic cylinder function; optical tweezers; quantumharmonic oscillator
Submitted to:
J. Phys. A: Math. Gen.
ONTENTS Contents1 Introduction 22 First exit time distribution 3
Conclusion 32Appendices 33References 431. Introduction
First passage time (FPT) distributions have found numerous applications in appliedmathematics, physics, biology and finance [1, 2, 3]. The FPT can characterizethe time needed for an animal to find food; the time for an enzyme to localizespecific DNA sequence and to initiate biochemical reaction; the time to exit froma confining domain (e.g., a maze); or the time to buy or sell an asset when itsprice deviation from the mean exceeds a prescribed threshold. The FPT distributionhas been studied for a variety of diffusive processes, ranging from ordinary diffusion(Brownian motion) to continuous-time random walks (CTRWs) [4, 5, 6, 7, 8, 9, 10, 11],fractional Brownian motion [12, 13, 14, 15], L´evy flights [16, 17, 18], surface-mediateddiffusion [19, 20, 21, 22] and other intermittent processes [23, 24], diffusion in scale
ONTENTS exit time(FET) distribution of a multi-dimensional Ornstein-Uhlenbeck process from a ball[49, 50, 51]. The probability distribution can be found through the inverse Fourier (resp.Laplace) transform of the characteristic (resp. moment-generating) function for whichexplicit representations in terms of special functions are well known [30, 39]. Althoughthe problem is formally solved, the solution involves confluent hypergeometric functionsand thus requires subtle asymptotic methods and computational hints. The aim of thereview is to provide a didactic self-consistent description of theoretical, numerical andpractical aspects of this problem.First, we recall the main derivation steps of the FET distribution, from the Langevinequation (Sec. 2.1), through forward and backward Fokker-Planck (FP) equations (Sec.2.2, 2.3), to spectral decompositions based on the eigenvalues and eigenfunctions of theFP operator (Sec. 2.4). This general formalism is then applied to describe the first exittimes of harmonically trapped particles in one dimension: the mean exit time (Sec. 2.5),the survival probability (Sec. 2.6), and the moment-generating function (Sec. 2.7). Inparticular, we analyze the asymptotic behavior of the mean exit time and eigenvalues indifferent limits (e.g., strong trapping potential, large constant force, etc.). Extensions tothe radial Ornstein-Uhlenbeck process in higher-dimensional cases for both interior andexterior problems are presented in Sec. 2.8 and Sec. 2.9, respectively. Although mostof these results are classical, their systematic self-contained presentation and numericalillustrations are missing.Section 3 starts from the summary of computational hints for computing confluenthypergeometric functions while technical details are reported in Appendix B. We discussthen three applications: (i) calibration of optical tweezers’ stiffness in single-particletracking experiments and detection of eventual constant forces exerted on a tracerby motor proteins (Sec. 3.2), (ii) adhesion bond dissociation under mechanical stress(Sec. 3.3), and (iii) distribution of triggering times of trend following strategies inalgorithmic trading on stock markets (Sec. 3.4). We also illustrate a direct relation tothe distribution of first crossing times of a moving boundary by Brownian motion (Sec.3.5). Finally, we present several extensions of the spectral approach, including diffusionunder quadratic double-well potential (Sec. 3.6) and anomalous diffusion (Sec. 3.7).Many technical details are summarized in Appendices.
2. First exit time distribution
We first recall the standard theoretical description of harmonically trapped particles byLangevin and Fokker-Planck equations [46, 52, 53]. We start with one-dimensionalOrnstein-Uhlenbeck process and then discuss straightforward extensions to higherdimensions.
ONTENTS We consider a diffusing particle of mass m trapped by a harmonic potential of strength k and pulled by a constant force F . The thermal bath surrounding the particle resultsin its stochastic trajectory which can be described by Langevin equation [52] m ¨ X ( t ) = − γ ˙ X ( t ) + F ( X ( t )) + ξ ( t ) , (1)where − γ ˙ X ( t ) is the viscous Stokes force ( γ being the drag constant), F ( X ( t )) = − kX ( t ) + F includes the externally applied Hookean and constant forces, and ξ ( t )is the thermal driving force with Gaussian distribution such that h ξ ( t ) i = 0 and h ξ ( t ) ξ ( t ′ ) i = 2 k B T γδ ( t − t ′ ), with k B ≃ . · − J/K being the Boltzmann constant, T the absolute temperature (in degrees Kelvin), δ ( t ) the Dirac distribution, and h . . . i denoting the ensemble average or expectation. In the overdamped limit ( m = 0), onegets the first-order stochastic differential equation˙ X ( t ) = 1 γ [ F ( X ( t )) + ξ ( t )] = kγ (ˆ x − X ( t )) + ξ ( t ) γ , X (0) = x , (2)where ˆ x = F /k is the stationary position (mean value), and x is the starting position.The Langevin equation can also be written in a conventional (dimensionless) stochasticform [39, 40] dX t = µ ( X t , t ) dt + σ ( X t , t ) dW t , X = x , (3)where W t is the standard Wiener process (Brownian motion), µ ( x, t ) and σ ( x, t ) are thedrift and volatility which in general can depend on x and t . In our case, the volatilityis constant, while the drift is a linear function of x , µ ( x, t ) = (ˆ x − x ) θ , i.e. dX t = θ (ˆ x − X t ) dt + σdW t , X = x , (4)where θ = kδ/γ , and σ = √ Dδ , with δ being a time scale, and D = k B T /γ thediffusion coefficient. This stochastic differential equation defines an Ornstein-Uhlenbeck(OU) process, with mean ˆ x , variance σ , and rate θ . An integral representation of Eq.(4) reads X t = x e − θt + ˆ x (1 − e − θt ) + σ t Z e θ ( t ′ − t ) dW t ′ . (5)One can see that X t is a Gaussian process with mean h X t i = x e − θt + ˆ x (1 − e − θt ) andcovariance h X t X t ′ i = σ θ ( e − θ | t − t ′ | − e − θ ( t + t ′ ) ).The discrete version of Eq. (4) with a fixed time step δ is known as auto-regressivemodel AR(1): X n = (1 − kδ/γ ) X n − + F δ/γ + √ Dδ ξ n , (6)where ξ n are standard iid Gaussian variables with mean zero and unit variance.This discrete scheme can be used for numerical generation of stochastic trajectories.An extension of the above stochastic description to multi-dimensional processes isstraightforward. ONTENTS The Langevin equation (2) expresses the displacement ˙ X ( t ) δ over a short time step δ in terms of the current position X ( t ). In other words, the distribution of the nextposition is fully determined by the current position, the so-called Markov property. Sucha Markov process can be characterized by a propagator or a transition density, i.e., theconditional probability density p ( x, t | x , t ) of finding the particle at x at time t , giventhat it was at x at earlier time t . The propagator can be seen as a “fraction” of pathsfrom x to x among all paths started at x (of duration t − t ) which formally writesas the average of the Dirac distribution δ ( X ( t ) − x ) over all random paths started from x : p ( x, t | x , t ) = h δ ( X ( t ) − x ) i X ( t )= x . The Markov property implies the Chapman-Kolmogorov (or Smoluchowski) equation p ( x, t | x , t ) = ∞ Z −∞ dx ′ p ( x, t | x ′ , t ′ ) p ( x ′ , t ′ | x , t ) ( t < t ′ < t ) , (7)which expresses a simple fact that any continuous path from X ( t ) = x to X ( t ) = x can be split at any intermediate time t ′ into two independent paths, from x to x ′ , andfrom x ′ to x .As a function of the arrival state ( x and t ), the propagator satisfies the forwardFokker-Planck (FP) equation [52, 53]. We reproduce the derivation of this equationfrom [54] which relies on the evaluation of the integral I = ∞ Z −∞ dx h ( x ) h p ( x, t + δ | x , t ) − p ( x, t | x , t ) i for any smooth function h ( x ) with compact support. One has I = ∞ Z −∞ dx h ( x ) ∞ Z −∞ dx ′ p ( x, t + δ | x ′ , t ) p ( x ′ , t | x , t ) − ∞ Z −∞ dx ′ h ( x ′ ) p ( x ′ , t | x , t ) ∞ Z −∞ dx p ( x, t + δ | x ′ , t )= ∞ Z −∞ dx ∞ Z −∞ dx ′ p ( x, t + δ | x ′ , t ) p ( x ′ , t | x , t ) h h ( x ) − h ( x ′ ) i , where the first term was represented using Eq. (7), while the normalization of theprobability density p ( x, t + δ | x ′ , t ) allowed one to add the integral over x in the secondterm. Expanding h ( x ) into a Taylor series around x ′ and then exchanging the integrationvariables x and x ′ , one gets I = ∞ Z −∞ dx p ( x, t | x , t ) ∞ X n =1 (cid:18) d n dx n h ( x ) (cid:19) n ! ∞ Z −∞ dx ′ p ( x ′ , t + δ | x, t ) ( x ′ − x ) n , ONTENTS n times, dividing by δ and taking the limit δ → ∞ Z −∞ dx h ( x ) ∂p ( x, t | x , t ) ∂t = ∞ Z −∞ dx h ( x ) ∞ X n =1 ( − n d n dx n (cid:16) D ( n ) ( x ) p ( x, t | x , t ) (cid:17) , where the left hand side is the limit of I/δ as δ →
0, and D ( n ) ( x ) = 1 n ! lim δ → δ ∞ Z −∞ dx ′ p ( x ′ , t + δ | x, t ) ( x ′ − x ) n . (8)Since the above integral relation is satisfied for arbitrary function h ( x ), one deduces theso-called Kramers-Moyal expansion: ∂p ( x, t | x , t ) ∂t = ∞ X n =1 ( − n d n dx n (cid:16) D ( n ) ( x ) p ( x, t | x , t ) (cid:17) . (9)Here we assumed that the process is time homogeneous, i.e., p ( x, t | x , t ) is invariantunder time shift: p ( x, t | x , t ) = p ( x, t + t ′ | x , t + t ′ ) that implies the time-independenceof D ( n ) ( x ).The density p ( x ′ , t + δ | x, t ) in Eq. (8) characterizes the displacement between X ( t ) = x and X ( t + δ ) = x ′ which can be written as X ( t + δ ) − X ( t ) ≃ δγ [ F ( x ) + ξ ( t )]for small δ according to the Langevin equation (2). After discretization in units of δ , the thermal force ξ ( t ) becomes a Gaussian variable with mean zero and variance2 k B T γ/δ . As a consequence, the displacement x ′ − x is also a Gaussian variable withmean ( δ/γ ) F ( x ) and variance ( δ/γ ) k B T γ/δ , i.e., p ( x ′ , t + δ | x, t ) = 1 √ πDδ exp (cid:18) − ( x ′ − x − F ( x ) δ/γ ) Dδ (cid:19) for small δ . Substituting this density into Eq. (8) and evaluating Gaussian integrals,one gets D (1) ( x ) = F ( x ) /γ , D (2) = D , and D ( n ) = 0 for n > ∂∂t p ( x, t | x , t ) = L x p ( x, t | x , t ) , L x = − ∂ x F ( x ) γ + D∂ x , (10)where L x is the Fokker-Planck operator acting on the arrival point x . This equationis completed by the initial condition p ( x, t | x , t ) = δ ( x − x ) at t = t , with a fixedstarting point x . Note that the forward FP equation can be seen as the probabilityconservation law, ∂∂t p ( x, t | x , t ) = − ∂ x J ( x, t | x , t ) , where J ( x, t | x , t ) = F ( x ) γ p ( x, t | x , t ) − D∂ x p ( x, t | x , t ) is the probability flux. Setting J = 0, one can solves the first-order differential equation to retrieve the equilibriumsolution p eq ( x ) = Zw ( x ), where Z is the normalization factor, and w ( x ) = exp x Z dx ′ F ( x ′ ) k B T = exp (cid:18) − V ( x ) k B T (cid:19) = exp (cid:18) − kx k B T + F xk B T (cid:19) , (11) ONTENTS V ( x ) = − R x dx ′ F ( x ′ ) is the potential associated to the force F ( x ). This is thestandard Boltzmann-Gibbs equilibrium distribution.When the FP operator L x has a discrete spectrum, the probability density admitsthe spectral decomposition p ( x, t | x , t ) = ∞ X n =0 v n ( x ) v n ( x ) ˜ w ( x ) e − λ n ( t − t ) (12)over the eigenvalues λ n and eigenfunctions v n ( x ) of L x : L x v n ( x ) + λ n v n ( x ) = 0 ( n = 0 , , , . . . ) (13)(eventually with appropriate boundary conditions, see below). The weight ˜ w ( x ) =1 /w ( x ) ensures the orthogonality of eigenfunctions: Z dx ˜ w ( x ) v m ( x ) v n ( x ) = δ m,n , (14)while the closure (or completeness) relation reads ∞ X n =0 v n ( x ) v n ( x ) ˜ w ( x ) = δ ( x − x ) . (15)This relation implies the initial condition p ( x, t | x , t ) = δ ( x − x ). As for the Langevinequation, an extension to the multi-dimensional case is straightforward. In particular,the derivative ∂ x is replaced by the gradient operator, while ∂ x becomes the Laplaceoperator [53, 55]. The forward FP equation describes the evolution of the probability density p ( x, t | x , t )from a given initial state (here, the starting point x at time t ). Alternatively, if theparticle is found at the arrival point x at time t (or, more generally, in a prescribedsubset of states), one can interpret p ( x, t | x , t ) as the conditional probability densitythat the particle is started from x at time t knowing that it arrived at x at later time t .As a function of x and t , this probability density satisfies the backward Fokker-Planck(or Kolmogorov) equation [53]: − ∂∂t p ( x, t | x , t ) = L ∗ x p ( x, t | x , t ) , (16)where the backward FP operator L ∗ is adjoint to the forward FP operator L (i.e.,( L f, g ) = ( f, L ∗ g ) for any two functions f and g from an appropriate functional space).Eq. (10) implies L ∗ x = F ( x ) γ ∂ x + D∂ x = kγ (ˆ x − x ) ∂ x + D∂ x . (17)Note that this operator acts on the starting point x while the sign minus in front oftime derivative reflects the backward time direction. Eq. (16) is easily obtained bydifferentiating the Champan-Kolmogorov equation (7) with respect to the intermediatetime t ′ . ONTENTS u n ( x ) of the backward FP operator L ∗ are simply u n ( x ) = v n ( x ) /w ( x ).As a consequence, one can rewrite the spectral decomposition (12) as p ( x, t | x , t ) = ∞ X n =0 u n ( x ) u n ( x ) w ( x ) e − λ n ( t − t ) , (18)with the weight w ( x ) from Eq. (11). The eigenfunctions u n ( x ) are as well orthogonal: Z dx w ( x ) u m ( x ) u n ( x ) = δ m,n , (19)while the closure (or completeness) relation reads ∞ X n =0 u n ( x ) u n ( x ) w ( x ) = δ ( x − x ) . (20)This relation implies the terminal condition p ( x, t | x , t ) = δ ( x − x ) at t = t . In contrastto Eq. (12), the weight w ( x ) in the spectral representation (18) depends on the fixed arrival point x , while the backward FP operator L ∗ x acts on eigenfunctions u n ( x ).When there is no force term, the operator L is self-adjoint, L = L ∗ , and theprobability density is invariant under time reversal: p ( x, t | x , t ) = p ( x , t | x, t ). Thisproperty does not hold in the presence of force.Finally, the backward FP equation is closely related to the Feynman-Kac formula fordetermining distributions of various Wiener functionals [56, 57, 58, 59, 60]. For instance,we already mentioned that the probability density p ( x, t | x , t ) can be understood asthe conditional expectation: p ( x, t | x , t ) = h δ ( X ( t ) − x ) i X ( t )= x . More generally, forgiven functions ψ ( x ), f ( x , t ) and U ( x , t ), the conditional expectation u ( x , t ) = (cid:28) exp (cid:18) − t Z t dt ′ U ( X ( t ′ ) , t ′ ) (cid:19) ψ ( X ( t )) (21)+ t Z t dt ′ f ( X ( t ′ ) , t ′ ) exp (cid:18) − t ′ Z t dt ′′ U ( X ( t ′′ ) , t ′′ ) (cid:19)(cid:29) X ( t )= x satisfies the backward FP equation − ∂∂t u ( x , t ) = L ∗ x u ( x , t ) − U ( x , t ) u ( x , t ) + f ( x , t ) , (22)subject to the terminal condition u ( x , t ) = ψ ( x ) at a later time t > t . In this review, we study the random variable τ = inf { t > | X ( t ) | > L } , i.e., the firstexit time of the process X ( t ) from an interval [ − L, L ] when started from x at t = 0.The cumulative distribution function of τ is also known as the survival probability S ( x , t ) = P { τ > t } up to time t of a particle which started from x . The notionof survival is associated to disappearing of the particle that hit either endpoint, due ONTENTS S ( x , t ) can be expressedthrough the probability density p ( x, t | x ,
0) of moving from x to x in time t withoutvisiting the endpoints ± L during this motion . Alternatively, p ( x, t | x ,
0) can be seenas the conditional probability density of starting from point x at time t = 0 undercondition to be at x at time t . This condition includes the survival up to time t , i.e., notvisiting the endpoints ± L . The probability density p ( x, t | x ,
0) satisfies the backwardFP equation with Dirichlet boundary condition at x = ± L : p ( x, t | ± L,
0) = 0. Thiscondition simply states that a particle started from either endpoint has immediatelyhit this endpoint, i.e. not survived. Note that this condition is preserved during allintermediate times t ′ due to the Chapman-Kolmogorov equation (7).Since the survival probability S ( x , t ) ignores the actual position x at time t , onejust needs to average the density p ( x, t | x ,
0) over x : S ( x , t ) = L Z − L dx p ( x, t | x ,
0) = ∞ X n =0 u n ( x ) e − λ n t L Z − L dx u n ( x ) w ( x ) , (23)where the spectral decomposition (18) was used. The eigenfunctions u n ( x ) of thebackward FP operator should satisfy Dirichlet boundary condition at x = ± L : u n ( ± L ) = 0. Eq. (23) also implies the backward FP equation ∂S ( x , t ) ∂t = L ∗ x S ( x , t ) , (24)which is completed by the initial condition S ( x ,
0) = 1 (the particle exists at thebeginning) and Dirichlet boundary condition S ( ± L, t ) = 0 (the process is stopped uponthe first arrival at either endpoint of the confining interval [ − L, L ]). Since the processis homogeneous in time, p ( x, t | x , t ) depends on t − t and thus ∂∂t p ( x, t | x , t ) = − ∂∂t p ( x, t | x , t ) that allows one to write the left-hand side of the backward FP equation(24) with the plus sign. Note that the characterization of first passage times throughthe backward FP equation goes back to the seminal work in 1933 by Pontryagin et al. [61]. Similar equations emerge in quantum mechanics when one searches for eigenstatesof a particle trapped by a short-range harmonic potential [62] (see also Appendix D forquantum harmonic oscillator).The FET probability density is q ( x , t ) = − ∂S ( x ,t ) ∂t , while the moment-generatingfunction is given by its Laplace transform: h e − sτ i = ∞ Z dt e − st q ( x , t ) ≡ ˜ q ( x , s ) , (25)with tilde denoting Laplace-transformed quantities. The Laplace transform of Eq. (24)yields the equation ( L ∗ x − s ) ˜ S ( x , s ) = − q ( x , s ) = 1 − s ˜ S ( x , s ), one gets( L ∗ x − s )˜ q ( x , s ) = 0 , (26)with Dirichlet boundary condition ˜ q ( ± L, s ) = 1.
ONTENTS h τ m i x can be found in one of standard ways:(i) from the moment-generating function, h τ m i x = ( − m lim s → ∂ m ∂s m ˜ q ( x , s ); (27)(ii) from the spectral representation of the survival probability: h τ m i x = m ! ∞ X n =0 u n ( x ) λ − mn L Z − L dx u n ( x ) w ( x ); (28)(iii) from recurrence partial differential equations (PDEs) L ∗ x h τ m i x = − m h τ m − i x , (29)with Dirichlet boundary conditions [30].In what follows, we focus on the mean exit time h τ i x , the moment-generating function˜ q ( x , s ), and the survival probability S ( x , t ) for harmonically trapped particles. The mean exit times of diffusive processes were studied particularly well because oftheir practical importance and simpler mathematical analysis (see [1, 2, 63, 64, 65] andreferences therein). In fact, the mean exit time, h τ i x = ∞ Z dt t q ( x , t ) = ∞ Z dt S ( x , t ) , (30)satisfies the simpler equation than the time-dependent PDE (16): L ∗ x h τ i x = − , (31)with Dirichlet boundary conditions at x = ± L . The double integration and imposedboundary conditions yield [66] ‡h τ i x = 1 D (cid:26)(cid:20)(cid:16) L Z − L dxw ( x ) (cid:17) − L Z − L dxw ( x ) x Z dx ′ w ( x ′ ) (cid:21) x Z − L dxw ( x ) − x Z − L dxw ( x ) x Z dx ′ w ( x ′ ) (cid:27) . (32)Substituting w ( x ) from Eq. (11), one gets h τ i x = L D √ π κ (cid:26) erf( i √ κ ( x /L − ϕ )) + erf( i √ κ (1 + ϕ ))erf( i √ κ (1 − ϕ )) + erf( i √ κ (1 + ϕ )) (33) × √ κ (1 − ϕ ) Z √ κ ( − − ϕ ) dz e z erf( z ) − √ κ ( x /L − ϕ ) Z √ κ ( − − ϕ ) dz e z erf( z ) (cid:27) , ‡ In [66], the sign minus in front of U ( z ) in the second integral in the numerator of the first term inEq. (7.7) is missing. ONTENTS −1 −0.5 0 0.5 100.20.40.60.81 z 〈 τ 〉 z / 〈 τ 〉 (a) κ = 0.1 κ = 1 κ = 10 −1 −0.5 0 0.5 100.20.40.60.8 z 〈 τ 〉 z (b) ϕ = 0 ϕ = 1 ϕ = 10 Figure 1.
Mean exit time h τ i z as a function of z = x /L : for different κ at fixed ϕ = 0 (a) and for different ϕ at fixed κ = 1 (b) . The timescale L /D is set to 1. Forplot (a) , the mean exit time is divided by its maximal value (at z = 0) in order torescale the curves. Circles indicate the mean exit time L D (1 − z ) without trapping( k = 0). where erf( z ) is the error function, and κ and ϕ are two dimensionless parameterscharacterizing the trapping harmonic potential and the pulling constant force,respectively κ ≡ kL k B T , ϕ ≡ ˆ xL = F kL . (34)Throughout the paper, we consider ϕ ≥ ϕ < ϕ → − ϕ and x → − x . For large κ or ϕ , one can use an equivalentrepresentation (A.1) provided in Appendix A.1.Several limiting cases are of interest: • When ϕ = 0 (i.e., F = 0), Eq. (33) is reduced to h τ i x = L D √ π κ √ κ Z √ κ x /L dz e z erf( z ) . (35) • In the limit k →
0, one gets a simpler expression h τ i x = L Dη (cid:16) − x /L − e − ηx /L − e − η e η − e − η (cid:17) , (36)where η = F L/ ( k B T ) is another dimensionless parameter. If F = 0, one retrieves theclassical result for Brownian motion: h τ i x = L D (cid:16) − ( x /L ) (cid:17) . (37) • For small κ , the Taylor expansion of Eq. (33) yields h τ i x ≃ L − x D (cid:16) κ − ϕ ( x /L ) + ( x /L ) O ( κ ) (cid:17) . (38)We emphasize that the limits κ → k → ϕ → ∞ according to Eq. (34). ONTENTS • In the opposite limit of large κ , four cases can be distinguished (seeAppendix A.1): h τ i x ≃ L D √ π e κ κ / ( ϕ = 0), √ π e κ (1 − ϕ ) κ / (1 − ϕ ) (0 < ϕ < κ ln √ κ (1 − x /L )0 . . . . ( ϕ = 1),12 κ ln ϕ − x /Lϕ − ϕ > x not too close to ± L . Note that the limit of the second asymptotic relation (for 0 < ϕ <
1) as ϕ → ϕ = 0 by factor 2. In fact, when ϕ >
0, it is much moreprobable to reach the right endpoint than the left one, and h τ i x characterizes mainlythe exit through the right endpoint at large κ . In turn, when ϕ = 0, both endpoints areequivalent that doubles the chances to exit and thus reduces by factor 2 the mean exittime. Note that the first two relations (up to a numerical prefactor) can be obtained bythe Kramers theory of escape from a potential well [66, 67]. The last relation in Eqs.(39) can be retrieved from the last line of Eq. (7.9) of Ref. [66].The behavior of the mean exit time h τ i x as a function of the starting point x isillustrated on Fig. 1. The increase of κ at fixed ϕ = 0 transforms the spatial profile ofthe mean exit time from the parabolic shape (37) at κ = 0 to a Π-shape at large κ (Fig.1a). In other words, the dependence on the starting point becomes weak at large κ . Atthe same time, the height of the profile rapidly grows with κ according to Eq. (39). Onthe opposite, the spatial profile becomes more skewed and sensitive to the starting pointas ϕ increases at fixed κ = 1, while the height is decreasing (Fig. 1b). As expected, theconstant force breaks the initial symmetry of the harmonic potential and facilitates theescape from the trap.Figure 2 shows how the mean exit time h τ i from the center varies with κ and ϕ .When there is no constant force ( ϕ = 0), one observes a rapid (exponential) growth atlarge κ , in agreement with Eq. (39) (shown by circles). The presence of a moderateconstant force (with 0 < ϕ <
1) slows down the increase of the mean exit time. Forinstance, at ϕ = 0 . h τ i exhibits a broad minimum at intermediate values of κ , but itresumes growing at larger values of κ . In turn, for ϕ ≥
1, there is no exponential growthwith κ , and the mean exit time slowly decreases, as expected from Eq. (39). Since theconstant force shifts the minimum of the harmonic potential from 0 to ˆ x = F /k , theborder value ϕ = 1 corresponds to the minimum ˆ x at the exit position (ˆ x = L ). For ϕ <
1, the harmonic potential keeps the particle away from the exit and thus greatlyincreases the mean exit time. In turn, for ϕ >
1, the harmonic potential attracts theparticle to ˆ x which is outside the interval [ − L, L ] and thus speeds up the escape.Although we considered the FET from a symmetric interval [ − L, L ] for convenience,shifting the coordinate by ˆ x allows one to map the original problem to the FET from a ONTENTS −1 −1 κ 〈 τ 〉 (a) ϕ = 0 ϕ = 0.5 ϕ = 1 −1 −2 −1 ϕ 〈 τ 〉 (b) κ = 0.1 κ = 1 κ = 10 Figure 2.
Mean exit time h τ i as a function of κ for fixed ϕ (a) and as a functionof ϕ for fixed κ (b) . The timescale L /D is set to 1. On both plots, circles show theexponential asymptotic relation in Eq. (39) for large κ and 0 ≤ ϕ <
1. On plot (a) ,crosses present the asymptotic Eq. (38) for small κ , to which the next-order term,2 κ /
45, is added. On plot (b) , crosses present the logarithmic asymptotic relation inEq. (39) for ϕ > nonsymmetric interval [ − a, b ] with a = L (1 + ϕ ) and b = L (1 − ϕ ), with the startingpoint x being shifted by Lϕ to vary from − a to b . As a consequence, the choice ofthe symmetric interval [ − L, L ] is not restrictive, and all the results can be recast for ageneral interval [ − a, b ] by shifts. The survival probability is fully determined by the eigenvalues and eigenfunctions of thebackward FP operator. The eigenvalue equation (13) reads Du ′′ − ( k/γ )( x − ˆ x ) u ′ + λu = 0 . (40)A general solution of this equation is well known [68] u ( z ) = c M (cid:18) − α κ , , κ ( z − ϕ ) (cid:19) + c ( z − ϕ ) M (cid:18) − α κ + 12 , , κ ( z − ϕ ) (cid:19) , (41)where z = x/L is the dimensionless coordinate, λ = Dα /L , c and c are arbitraryconstants, and M ( a, b, z ) = F ( a, b, z ) = ∞ X n =0 a ( n ) z n b ( n ) n ! (42)is the confluent hypergeometric function of the first kind (also known as Kummerfunction), with a (0) = 1 and a ( n ) = a ( a + 1) . . . ( a + n −
1) = Γ( a + n )Γ( a ) , where Γ( z ) isthe gamma function. The first and second terms in Eq. (41) are respectively symmetricand antisymmetric functions with respect to ϕ .To shorten notations, we set m (1) α,κ ( z ) ≡ M (cid:18) − α κ , , κz (cid:19) , (43) m (2) α,κ ( z ) ≡ zM (cid:18) − α κ + 12 , , κz (cid:19) , (44) ONTENTS u ( z ) = c m (1) α,κ ( z − ϕ ) + c m (2) α,κ ( z − ϕ ) . (45)The Dirichlet boundary conditions read c m (1) α,κ ( − − ϕ ) + c m (2) α,κ ( − − ϕ ) = 0 (at x = − L ) ,c m (1) α,κ (1 − ϕ ) + c m (2) α,κ (1 − ϕ ) = 0 (at x = L ) . In the special case ϕ = 1, one gets c = 0, and the eigenvalues are determined fromthe equation m (2) α,κ (2) = 0. In general, for ϕ = 1, one considers the determinant of theunderlying 2 × D α,κ,ϕ = m (1) α,κ ( − − ϕ ) m (2) α,κ (1 − ϕ ) − m (2) α,κ ( − − ϕ ) m (1) α,κ (1 − ϕ ) . (46)Setting this determinant to 0 yields the equation on α : D α n ,κ,ϕ = 0 , (47)where α n ( n = 0 , , , . . . ) denote all positive solutions of this equation (for fixed κ and ϕ ). The eigenfunctions read then u n ( z ) = β n √ L h c (1) n m (1) α n ,κ ( z − ϕ ) − c (2) n m (2) α n ,κ ( z − ϕ ) i , (48)where c (1) n = m (2) α n ,κ (1 − ϕ ) , c (2) n = m (1) α n ,κ (1 − ϕ ) , (49)and the normalization constant is β − n = − ϕ Z − − ϕ dz e − κz h c (1) n m (1) α n ,κ ( z ) − c (2) n m (2) α n ,κ ( z ) i . (50)Multiplying Eq. (40) by w ( x ) and integrating from a to b , one gets b Z a dx u ( x ) w ( x ) = Dλ [ u ′ ( a ) w ( a ) − u ′ ( b ) w ( b )] . (51)The derivative of the Kummer function can be expressed through Kummer functions,in particular, ∂ z m (1) α,κ ( z ) = α κz (cid:16) m (1) α,κ ( z ) − m (1) √ α − κ,κ ( z ) (cid:17) , (52) ∂ z m (2) α,κ ( z ) = (cid:18) κz − − α κ (cid:19) m (2) α,κ ( z ) + (cid:18) α κ (cid:19) m (2) √ α +4 κ,κ ( z ) , (53)from which one gets explicit formulas for u ′ n ( a ) and u ′ n ( b ) and thus for the integral inEq. (51). We get therefore S ( x , t ) = ∞ X n =0 w n e − Dtα n /L h c (1) n m (1) α n ,κ ( x /L − ϕ ) − c (2) n m (2) α n ,κ ( x /L − ϕ ) i , (54) ONTENTS t q ( ,t ) (a) κ = 0 κ = 0.1 κ = 1 κ = 2 t q ( ,t ) (b) ϕ = 0 ϕ = 0.1 ϕ = 0.5 ϕ = 0.9 Figure 3.
FET probability density q (0 , t ) for several κ at fixed ϕ = 0 (a) andfor several ϕ at fixed κ = 1 (b) . The timescale L /D is set to 1. The spectraldecomposition (57) is truncated after 30 terms. where w n = β n e − κ α n [ v ( − − v (1)] , (55)with v ( z ) = e κϕz (cid:16) c (1) n ∂ z m (1) α n ,κ ( z − ϕ ) − c (2) n ∂ z m (2) α n ,κ ( z − ϕ ) (cid:17) . (56)Taking the derivative with respect to time, one obtains the FET probability density q ( x , t ) = DL ∞ X n =0 w n α n e − Dtα n /L h c (1) n m (1) α n ,κ ( x /L − ϕ ) − c (2) n m (2) α n ,κ ( x /L − ϕ ) i . (57)In the limit κ →
0, functions m (1) α,κ ( z ) and m (2) α,κ ( z ) approach cos( αz ) and sin( αz ),respectively, so that eigenfunctions from Eq. (48) become u n ( z ) = β n √ L sin( α (1 − z )),while the determinant in Eq. (46) is reduced to sin(2 α ), from which α n = π ( n + 1) / ϕ vanishes, and one retrieves the classical result forBrownian motion S ( x , t ) = 2 ∞ X n =0 ( − n e − Dtπ ( n +1 / /L π ( n + 1 /
2) cos( π ( n + 1 / x /L ) . (58)Only symmetric eigenfunctions with α n = π ( n + 1 /
2) contribute to this expression.For centered harmonic potential ( ϕ = 0), Eq. (47) is reduced to m (1) α n ,κ (1) m (2) α n ,κ (1) = 0 , (59)which determines two sequences of zeros: α n, from m (1) α n, ,κ (1) = 0, and α n, from m (2) α n, ,κ (1) = 0. As a consequence, one can consider separately two sequences ofsymmetric and antisymmetric eigenfunctions: m (1) α n, ,κ ( z ) and m (2) α n, ,κ ( z ). According toEq. (23), integration over arrival points removes all the terms containing antisymmetriceigenfunctions. This simpler situation is considered as a particular case in Sec. 2.8.Figure 3 illustrates the behavior of the probability density q ( x , t ). For fixed ϕ = 0,an increase of κ increases the mean exit time and makes the distribution wider. Notethat the most probable FET remains almost constant. The opposite trend appears for ONTENTS −1 −0.5 0 0.5 100.20.40.60.81 z S ( z ,t ) (a) t = 0.002t = 0.01t = 0.1t = 0.5 −1 −0.5 0 0.5 100.20.40.60.81 z S ( z ,t ) (b) t = 0.002t = 0.01t = 0.1t = 0.5 Figure 4.
Survival probability S ( x , t ) as a function of the starting point z = x /L ,with κ = 1, and ϕ = 0 (a) and ϕ = 0 . (b) . The timescale L /D is set to 1. Thespectral decomposition (54) is truncated after 30 terms. variable ϕ at fixed κ = 1: an increase of ϕ diminishes the mean exit time and makesthe distribution narrower. This is expected because a strong constant force would drivethe particle to one exit and dominate over stochastic part.Figure 4 shows the dependence of the survival probability S ( x , t ) on the startingpoint x . At short times, the survival probability is close to 1 independently of x ,except for close vicinity of the endpoints. As time increases, S ( x , t ) is progressivelyattenuated. The spatial profile is symmetric for centered harmonic potential ( ϕ = 0),and skewed to the left in the presence of a positive constant force ( ϕ = 0 . Since any linear combination of functions in Eq. (45) satisfies Eq. (26), with s = − Dα /L , one can easily find the moment-generating function ˜ q ( x , s ) by imposingthe boundary condition ˜ q ( ± L, s ) = 1:˜ q ( x , s ) = A (1) α,κ,ϕ D α,κ,ϕ m (1) α,κ ( x /L − ϕ ) + A (2) α,κ,ϕ D α,κ,ϕ m (2) α,κ ( x /L − ϕ ) , (60)where A (1) α,κ,ϕ = m (2) α,κ (1 − ϕ ) − m (2) α,κ ( − − ϕ ) ,A (2) α,κ,ϕ = m (1) α,κ ( − − ϕ ) − m (1) α,κ (1 − ϕ ) , and D α,κ,ϕ is defined by Eq. (46). Setting a = L (1 + ϕ ) and b = L (1 − ϕ ), one retrievesthe moment-generating function of the FET of an Ornstein-Uhlenbeck process froman interval [ − a, b ] reported in [39] (p. 548, 3.0.1), in which Eq. (60) is written morecompactly in terms of two-parametric family S ( ν, a, b ) of parabolic cylinder functions(see Appendix B.1). For symmetric interval [ − a, a ], a similar expression for the moment-generating function was provided in [30].It is worth noting that the probability density q ( x , t ) could be alternatively foundby the inverse Laplace transform of Eq. (60). For this purpose, one determines the ONTENTS s n of ˜ q ( x , s ) in the complex plane which are given by zeros α n of D α,κ,ϕ accordingto Eq. (47). In other words, one has s n = − Dα n /L , and the residue theorem yields q ( x , t ) = 4 κDL ∞ X n =0 e − Dtα n /L h A (1) α n ,κ,ϕ D ′ α n ,κ,ϕ m (1) α n ,κ ( x /L − ϕ ) + A (2) α n ,κ,ϕ D ′ α n ,κ,ϕ m (2) α n ,κ ( x /L − ϕ ) i , (61)where D ′ α,κ,ϕ denotes the derivative of D α,κ,ϕ with respect to s = − α / (4 κ ). Comparingthe above formula to Eq. (57), one gets another representation for coefficients w n w n = 4 κα n A (1) α n ,κ,ϕ D ′ α n ,κ,ϕ , (62)where we used the identity c (1) n A (2) α n ,κ,ϕ = − c (2) n A (1) α n ,κ,ϕ , with c (1 , n from Eq. (49).Two alternative representations (55) and (62) allow one to compute the normalizationcoefficients β n without numerical integration in Eq. (50). In higher dimensions, we consider the FET of a multi-dimensional Ornstein-Uhlenbeckprocess from a ball of radius L . For centered harmonic potential (i.e., F = 0), thederivation follows the same steps as earlier. In fact, the integration of the probabilitydensity p ( x, t | x ,
0) over the arrival point x in the multi-dimensional version of Eq.(23) removes the angular dependence of the survival probability so that the eigenvalueequation is reduced to the radial part (cid:18) D (cid:20) ∂ r + d − r ∂ r (cid:21) − krγ ∂ r (cid:19) u n ( r ) + λ n u n ( r ) = 0 , (63)where d is the space dimension. In other words, we consider the FPT of the radialOrnstein-Uhlenbeck process to the level L . In turn, the analysis for non-centeredharmonic potential with F = 0 is much more involved in higher dimensions due toangular dependence, and is beyond the scope of this review. Survival probability.
A solution of Eq. (63) is given by the Kummer function which isregular at r = 0 u n ( r ) = β n L d/ M (cid:16) − α n κ , d , κ ( r/L ) (cid:17) ( n = 0 , , , . . . ) , (64)where β n is the normalization factor: β − n = Z dz z d − e − κz (cid:20) M (cid:16) − α n κ , d , κz (cid:17)(cid:21) . (65)The eigenvalues λ n = Dα n /L are determined by the positive zeros α n of the equation M (cid:16) − α n κ , d , κ (cid:17) = 0 . (66) ONTENTS S ( r , t ) = ∞ X n =0 w n e − Dtα n /L M (cid:18) − α n κ , d , κ ( r /L ) (cid:19) , (67)where w n = β n e − κ κ M (cid:18) − α n κ + 1 , d , κ (cid:19) , (68)and we used the identity Z dz z d − e − κz M (cid:18) − α n κ , d , κz (cid:19) = e − κ κ M (cid:18) − α n κ + 1 , d , κ (cid:19) . (69)The FET probability density is then q ( r , t ) = DL ∞ X n =0 w n α n e − Dtα n /L M (cid:18) − α n κ , d , κ ( r /L ) (cid:19) . (70)In the limit κ →
0, one can use the identity (see Appendix B.1)lim κ → M (cid:18) − α κ , d , κz (cid:19) = Γ( d/ J d/ − ( αz )( αz/ d/ − = cos( αz ) ( d = 1) J ( αz ) ( d = 2)sin( αz ) αz ( d = 3) (71)to retrieve the classical results for Brownian motion (here J n ( z ) is the Bessel functionof the first kind). In particular, one retrieves α n = π ( n + 1 /
2) in one dimension and α n = π ( n + 1) in three dimensions (with n = 0 , , , . . . ). For the one-dimensional case,we retrieved only the zeros of symmetric eigenfunctions that contribute to the survivalprobability (cf. discussion in Sec. 2.6). Moment-generating function.
The moment-generating function, obeying Eq. (63) with − s instead of λ n , is˜ q ( r , s ) = M ( sL κD , d , κ r L ) M ( sL κD , d , κ ) , (72)in agreement with [39] (p. 581, 2.0.1). This function satisfies the boundary condition˜ q ( L, s ) = 1 and is regular at r = 0. The Laplace inversion of this expression yieldsanother representation of the probability density q ( r , t ) = 4 κDL ∞ X n =0 e − Dtα n /L M ( − α n κ , d , κ r L ) M ′ ( − α n κ , d , κ ) , (73)where M ′ ( a, b, z ) denotes the derivative of M ( a, b, z ) with respect to a . Comparing thisrelation to Eq. (70), the coefficients w n from Eq. (68) can also be identified as w n = 4 κα n M ′ ( − α n κ ; d ; κ ) . (74)As mentioned above, two expressions (68, 74) for w n can be used to compute thenormalization constants β n without numerical integration in Eq. (65). ONTENTS Mean exit time.
Following the same steps as in Sec. 2.5, one gets the mean exit timefor the higher-dimensional case h τ i r = L D κ √ κ Z √ κr /L dr r − d e r r Z dr r d − e − r , (75)where we imposed Dirichlet boundary condition at r = L and the regularity conditionat r = 0. In the limit κ →
0, one retrieves the classical result h τ i r = ( L − r ) / (2 dD ).In the opposite limit κ ≫
1, one gets h τ i r ≃ L D Γ( d/ e κ κ d/ ( κ ≫ , (76)which is applicable for any r not too close to L . The behavior of the mean exit timefor general spherically symmetric potentials is discussed in [66].In Appendix A.2, the asymptotic behavior of the smallest eigenvalue λ = Dα /L is obtained: λ ≃ DL κ d/ Γ( d/ e − κ ( κ ≫ , (77)which is just the inverse of the above asymptotic relation for the mean exit time. Whilethe first eigenvalue exponentially decays with κ , the other eigenvalues linearly grow with κ (see Appendix A.2): λ n ≃ DL κn ( κ ≫ . (78)As a consequence, the gap between the lowest eigenvalue λ and the next eigenvalue λ grows linearly with κ . For t ≫ /λ , the contribution of all excited eigenmodes becomesnegligible as compared to the lowest mode, and the first exit time follows approximatelyan exponential law, P { τ > t } ≃ exp( − t/ h τ i ), with the mean h τ i from Eq. (76).We illustrate this behavior for three-dimensional case on Fig. 5a which presentsthe first three eigenvalues λ n as functions of κ . Note that the correction term to theasymptotic line for λ is significant even for κ = 10 (see Appendix A.2 for details). Forcomparison, Fig. 5b shows the first three eigenvalues for the exterior problem discussedin the next subsection. For the exterior problem, when the process is started outside the interval [ − L, L ] (oroutside the ball of radius L in higher dimensions), the FET is also referred to as thefirst passage time to the boundary of this domain: τ = inf { t > | X ( t ) | < L } . Whilethe mean exit time and the probability distribution can be found in a very similarway (see below), their properties are very different from the earlier considered interiorproblem. For the sake of simplicity, we only consider the centered harmonic potential(i.e., F = 0), although the noncentered case in one dimension can be treated similarly.In one dimension, the domain ( −∞ , − L ) ∪ ( L, ∞ ) is split into two disjointsubdomains so that τ is in fact the first passage time to a single barrier, either at ONTENTS κ λ n (a) n = 0n = 1n = 2 κ λ n (b) n = 0n = 1n = 2 Figure 5.
First three eigenvalues λ n of the Fokker-Planck operator as functions of κ , for the interior problem (a) and the exterior problem (b) in three dimensions. Thetimescale L /D is set to 1. On plot (a) , crosses present the asymptotic relation (77)for the first eigenvalue λ while thin solid lines indicate the asymptotic relation 4 κn for higher eigenvalues ( n = 1 , , . . . ). At κ = 0, one retrieves the eigenvalues π ( n + 1) for Brownian motion. On plot (b) , thin lines indicate the asymptotic behavior (A.18)at small κ . x = L (if started from x > L ), or at x = − L (if started from x < − L ). This situationis described in Appendix C. Mean exit time.
Following the steps of Sec. 2.5, one obtains the mean exit time h τ i r = L D κ √ κ r /L Z √ κ dr r − d e r ∞ Z r dr r d − e − r , (79)where we imposed Dirichlet boundary condition at r = L and the regularity conditionat infinity.For even dimensions d , the change of integration variables yields the explicit formula h τ i r = L Dκ r /L ) + d − X j =1 Γ( d ) κ − j j Γ( d − j ) (1 − ( r /L ) − j ) (80)(we use the convention that P nj =1 a n is zero if n < h τ i r = L Dκ ln( r /L ) ( d = 2) . (81)For odd d , repeated integration by parts yields h τ i r = L D κ (cid:26) √ π √ κ z Z √ κ dz e z erfc( z ) + d − X j =1 Γ( d )Γ( d − j ) − Γ( j + )Γ( ) ! − z − j j κ j (82)+ e κ erfc( √ κ ) d − X j =1 Γ (cid:16) d − j (cid:17) κ j − d/ − e κz erfc( √ κ z ) d − X j =1 Γ (cid:16) d − j (cid:17) ( κz ) j − d/ (cid:27) , ONTENTS z = r /L , and erfc( z ) = 1 − erf( z ). Note that for d = 1, all terms vanish exceptthe integral.For large r or large κ , the leading asymptotic term is L Dκ ln( r /L ) for all dimensions(for odd dimensions, this term comes from the integral). When r /L approaches 1, themean exit time vanishes as c ( r /L − c depends on κ and d .When κ →
0, the mean exit time diverges: h τ i r ≃ L D Γ( d )2( d −
2) (1 − ( r /L ) − d ) κ − d/ ( d = 2) . (83)(for d = 2, see Eq. (81)). This divergence is expected because, for the exterior problem,the mean exit time for Brownian motion is infinite in all dimensions, irrespectively ofits recurrent or transient character.Finally, the mean exit time for non-centered harmonic potential (i.e., F = 0) inone dimension reads for x > L as h τ i x = L D √ π κ √ κ ( x /L − ϕ ) Z √ κ (1 − ϕ ) dz e z erfc( z ) . (84)In the limit of large κ , two asymptotic regimes are distinguished:(i) when ϕ <
1, the upper and lower limits of integration go to infinity so that themean exit time behaves as h τ i x ≃ L D κ ln x /L − ϕ − ϕ ; (85)(ii) when ϕ >
1, the lower limit of integration goes to −∞ , and the mean exit timeexponentially diverges as h τ i x ≃ L D √ π e κ ( ϕ − κ / ( ϕ − . (86)Both regimes are similar to that of the interior problem considered in Sec. 2.5. Probability distribution.
The moment-generating function ˜ q ( r , s ) for the exteriorproblem satisfies the same equation (63), with − s instead of λ n , as ˜ q ( r , s ) from Eq.(72) for the interior problem. In order to ensure the regularity condition at infinity (as r → ∞ ), one replaces M ( a, b, z ) by the confluent hypergeometric function of the secondkind (also known as Tricomi function): U ( a, b, z ) = Γ(1 − b )Γ( a − b + 1) M ( a, b, z ) + Γ( b − a ) z − b M ( a − b + 1 , − b, z ) (87)(for integer b , this relation is undefined but can be extended by continuity, seeAppendix B.2). For a >
0, the function U ( a, b, z ) vanishes as z → ∞ , in contrastto an exponential growth of M ( a, b, z ) according to Eqs. (B.7, B.8). In turn, U ( a, b, z )exhibits non-analytic behavior near z = 0, U ( a, b, z ) ≃ Γ(1 − b )Γ( a − b +1) + Γ( b − a ) z − b + . . . , thatlimited its use for the interior problem. ONTENTS q ( r , s ) = U ( sL κD , d , κ r L ) U ( sL κD , d , κ ) ( r ≥ L ) , (88)in agreement with [39] (p. 581, 2.0.1).Denoting by α n the positive zeros of the equation U (cid:18) − α n κ , d , κ (cid:19) = 0 , (89)the inverse Laplace transform yields the FET probability density: q ( r , t ) = 4 κDL ∞ X n =0 e − Dtα n /L U ( − α n κ , d , κ r L ) U ′ ( − α n κ , d , κ ) , (90)where U ′ ( a, b, z ) is the derivative with respect to a . Its integral over time is the survivalprobability: S ( r , t ) = ∞ X n =0 w n e − Dtα n /L U (cid:16) − α n κ , d , κ r L (cid:17) , (91)with w n = 4 κα n U ′ ( − α n κ , d , κ ) . (92)Alternatively, one can use the eigenvalues λ n = Dα n /L and the correspondingeigenfunctions u n ( r ) = β n L d/ U (cid:18) − α n κ , d , κ ( r/L ) (cid:19) , (93)where β n is the normalization factor: β − n = ∞ Z dz z d − e − κz (cid:20) U (cid:16) − α n κ , d , κz (cid:17)(cid:21) . (94)The normalization factors β n diverge as κ → w n = β n e − κ κ U (cid:18) − α n κ + 1 , d , κ (cid:19) . (95)The asymptotic behavior of eigenvalues as κ → In spite of apparent similarities between the interior and the exterior problems, there isa significant difference in spectral properties of two problems. This difference becomesparticularly clear in the limit κ → ONTENTS α n → π ( n + 1) / κ → κ → λ n vanish as κ → κ > κ = 0 are drastically different. One cansee that the asymptotic behavior of the eigenvalues λ n is quite different for the interiorand the exterior problems.
3. Discussion
In this section, we discuss computational hints for confluent hypergeometric functions(Sec. 3.1), three applications in biophysics and finance (Sec. 3.2, 3.3, 3.4), relation tothe distribution of first crossing times of a moving boundary by Brownian motion (Sec.3.5), diffusion under quadratic double-well potential (Sec. 3.6), and further extensions(Sec. 3.7).
The probability distribution of first exit times involves confluent hypergeometricfunctions M ( a, b, z ) (for interior problem) and U ( a, b, z ) (for exterior problem). Forinstance, the eigenvalues of the Fokker-Planck operator are obtained through zeros α n ofthe equation M ( − α n κ , d , κ ) = 0 or similar. As a consequence, one needs to compute thesefunctions for large | a | = α n / (4 κ ). Although the series in the definition (42) of M ( a, b, z )converges for all z , numerical summation becomes inaccurate for large | a | , and otherrepresentations of confluent hypergeometric functions are needed. In Appendix B.2, wediscuss an efficient numerical scheme for rapid and accurate computation of M ( a, b, z )for large | a | and moderate z which relies on the expansion (B.14). Moreover, we showthat this scheme is as well applicable for computing the derivative of M ( a, b, z ) withrespect to a which appears in Eq. (73) or similar after the inverse Laplace transform.For non-integer b , the Tricomi function U ( a, b, z ) is expressed through M ( a, b, z ) byEq. (87) that allows one to apply the same numerical scheme for the exterior problemin odd dimensions d . Although the Tricomi function for integer b can be obtained bycontinuation, the derivation of its rapidly converging representation is more subtle. Inpractice, one can compute U ( a, b, z ) for an integer b by extrapolation of a sequence U ( a, b ε , z ) computed for non-integer b ε approaching b as ε → z and moderate | a | , one can use integral representations of M ( a, b, z ) and U ( a, b, z ) (see Appendix B.2). The same algorithms can also be applied tocompute parabolic cylinder function D ν ( z ) and Whittaker functions (see Appendix B.1).The Matlab code for computing both Kummer and Tricomi functions is available § § See http://pmc.polytechnique.fr/pagesperso/dg/confluent/confluent.html . ONTENTS The Langevin equation (2) can describe the thermal motion of a small tracer in aviscous medium. The Hookean force − kX ( t ) incorporates the harmonic potential ofan optical tweezer which are used to trap the tracer in a specific region of the medium[69]. Optical trapping strongly diminishes the region accessible to the tracer and thusenables to reduce the field of view and to increase the acquisition rate up to few MHz[70, 71, 72, 73, 74]. At the same time, trapping affects the intrinsic dynamics of thetracer and may screen or fully remove its features at long times. The choice of thestiffness k is therefore a compromise between risk of loosing the tracer from the field ofview (too small k ) and risk of suppressing important dynamical features (too large k ).The FET statistics can then be used for estimating the appropriate stiffness due to aquantitative characterization of escape events. For instance, one can choose the stiffnessto ensure that the mean exit time strongly exceeds the duration of experiment, or thatthe escape probability is below a prescribed threshold.Another interesting option consists in detecting events in which a constant force isapplied to the tracer. In living cells, such events can mimic the action of motor proteinsthat attach to the tracer and pull it in one direction [75, 76, 77, 78, 79]. The presenceof a constant force facilitates the escape from the optical trap while higher fraction ofescape events (as compared to the case without constant force) can be an indicator ofsuch active transport mechanisms.Originally, the idea of fast escape in the case of comparable Hookean and externalforces was used to estimate the force generated by a single protein motor [75]. A“trap and escape” experiment consisted in trapping a single organelle moving alongmicrotubules at strong stiffness and then gradually reducing it until the organelle escapesthe trap. Repeating such measurement, one can estimate the “escape power” kL as ameasure of the driving force F when ϕ = F / ( kL ) ∼
1, where L is the size of the trap.In this way, the driving force generated by a single (presumably dynein-like) motor wasestimated to be 2.6 pN [75]. This approximate but direct way of force measurementrelies on the drastic change in the mean exit time behavior at ϕ = 1 according to Eq.(39).Interestingly, these mechanisms can even be detected from a single trajectory.When there is no constant force, the mean-square displacement (MSD) of a trappedtracer, h (∆ X ) i , is known to approach the constant level 2 k B T /k [52, 80, 81]. In otherwords, the square root of the long-time asymptotic MSD determines the typical size ℓ k = p k B T /k = √ Dτ k of the confining region due to optical trapping, with τ k = γ/k .Setting L = qℓ k , one gets κ = q , i.e., the dimensionless parameter κ can be interpretedas the squared ratio between the exit distance L and the characteristic size of the trap ℓ k . The above analysis showed that a tracer can rapidly reach levels which are below orslightly above ℓ k . However, significantly longer explorations are extremely improbable.In fact, according to Eq. (76), the mean exit time for κ ≫ h τ i ≃ τ k Γ( d/ e κ κ d/ ( κ ≫ , (96) ONTENTS L = qℓ k = √ κ √ Dτ k . For large enough t (i.e., t ≫ τ k ), the contributions ofall excited eigenstates vanish, and the survival probability exhibits a mono-exponentialdecay: S (0 , t ) ≃ exp( − t/ h τ i ), where we replaced the smallest eigenvalue λ by 1 / h τ i for κ ≫ w ≈ τ k ≪ t ≪ h τ i , the survival probability remains therefore close to1. A constant force F pulling the tracer from the optical trap strongly affects themean exit time and the survival probability. The dimensionless parameter ϕ from Eq.(34) is the ratio between the new stationary position ˆ x of the trajectory and the exitlevel L = √ κ ℓ k : ϕ = F kL = ˆ xℓ k √ κ = F √ Dτ k k B T √ κ . (97)For large ϕ , the mean exit time can be approximated according to Eq. (39) as h τ i ≃ τ k ln ϕϕ − ≃ τ k /ϕ , i.e., it becomes smaller than τ k , and much smaller than themean exit time from Eq. (96) without force. As expected, exit events would be observedmuch more often in the presence of strong constant force.For a long acquired trajectory, one can characterize how often different levelsare reached. Strong deviations from the expected statistics (given by the survivalprobability) would suggest the presence of a constant force. To illustrate this idea,we simulate the thermal motion of a spherical tracer of radius a = 1 µ m submerged inwater and trapped by an optical tweezer with a typical stiffness constant k = 10 − N/m[73, 74]. The Stokes relation implies γ = 6 πaη ≈ . · − kg/s, from whichthe diffusion coefficient is D = k B T /γ ≈ . · − m /s at T = 300 K (with η ≈ − kg/m/s being the water viscosity). The characteristic trapping time is τ k = γ/k ≈ . ℓ k = p k B T /k ≈
91 nm. Figure6a shows one simulated trajectory of the tracer. According to Eq. (33), the mean exittimes from the intervals ( − ℓ k , ℓ k ) and ( − ℓ k , ℓ k ) are 27 . ± ℓ k and only few crossings of levels ± ℓ k . For comparison, we generated another trajectoryfor which a constant force F = 0 . ϕ = 2 . / √ κ ) is applied between 0 . . . . x = 200 nm so that the trajectory crosses the level ℓ k and remains above this level for whole duration of the forced period. Once the force isswitched off, the trajectory returns to its initial regime with zero mean. One can seethat the use of FET statistics presents a promising perspective for design and analysis ofsingle-particle tracking experiments, while Bayesian techniques can be further applied toget more reliable results [82, 83]. Note that the FPT statistics have also been suggestedas robust estimators of diffusion characteristics [84] (see also [26]). Another method ONTENTS t (s) X ( t ) ( n m ) (a) t (s) X ( t ) ( n m ) (b) Figure 6.
Two simulated trajectories of a spherical tracer submerged in water underthe optical trapping: (a) no constant force; (b) constant force F = 0 . . . ± ℓ k = 91 nm, while dashed red line shows the stationary levelˆ x (equal to 0 for zero force and 200 nm for F = 0 . relying on the time evolution of the tracer probability distribution was proposed forsimultaneously extracting the restoring-force constant and diffusion coefficient [85].At the same time, we emphasize that this perspective needs further analysis.First, we focused on normal diffusion in a harmonic potential while numerous single-particle tracking experiments evidenced anomalous diffusion in living cells and polymersolutions [73, 86, 87, 88, 89]. Several theoretical models have been developed todescribe anomalous processes such as continuous-time random walks (CTRW), fractionalBrownian motion (fBm), and generalized Langevin equation [5, 6, 78, 79]. While anextension of the presented results is rather straightforward for CTRW (Sec. 3.7), theFPT problems for non-Markovian fBm or generalized Langevin equation are challengingdue to lack of equivalent Fokker-Planck formulation. Second, the quadratic profile isan accurate approximation for optical trapping potential only for moderate deviationsfrom the center of the laser beam [69], while the spatial profile can be more complicatedfor strong deviations. In other words, an accurate description of the tracer escape mayrequire more sophisticated analysis. Finally, the inference of constant forces from asingle trajectory may present some statistical challenges because different escape eventscan be correlated. We briefly mention another biophysical example of bond dissociation. Adhesion betweencells or of cells to surfaces is mediated by weak noncovalent interactions. Whilea reversible bond between two molecules can break spontaneously (due to thermalfluctuations), an external force is needed to rupture multiple bonds that link two cellstogether [90]. The dynamics of bond rupture can be seen as the first exit time problemin which exit or escape occurs when the intermolecular distance exceeds an effectiveinteraction radius. Bell suggested to apply the kinetic theory of the strength of solids
ONTENTS t b = t exp[( E b − r b F ) / ( k B T )] , (98)where E b is the bond energy, r b is the range of the minimum of the binding free energy, F is the applied external force per bond, and t is the lifetime at the critical force E b /r b at which the minimum of the free energy vanishes [90]. This relation became a canonicaldescription of adhesion bond dissociation under force.If the binding potential can be approximated as quadratic, then the lifetime of abond is precisely the mean exit time h τ i of a harmonically trapped particle. In thatcase, the second asymptotic relation in Eqs. (39) implies the quadratic dependence onthe force, h τ i ∼ e κ (1 − ϕ ) , where ϕ = F / ( kr b ), and κ = E b / ( k B T ) = kr b / (2 k B T ). Inother words, Eq. (98) is retrieved only for weak forces when the quadratic term ϕ canbe neglected. However, in the regime where the bond is most likely to break, the appliedforce is large, and the mean exit time may have completely different asymptotics (see,e.g., the last line of Eqs. (39) for ϕ > Algorithmic trading is another field for applications of FETs. In algorithmic trading,a set of trading rules is developed in order to anticipate the next price variation of anasset from its earlier (historical) prices [97]. Although the next price is random (andthus unpredictable), one aims to catch some global or local trends which can be inducedby collective behavior of multiple traders or macroeconomic tendences [98, 99, 100].Many trading strategies rely on the exponential moving average ¯ p n of the earlier prices p k [101, 102, 103, 104]¯ p n = λ ∞ X k =0 (1 − λ ) k p n − k , (99)where 0 < λ ≤ p n and the “anticipated” average price¯ p n , δ n ≡ p n − ¯ p n = (1 − λ ) ∞ X k =0 (1 − λ ) k r n − k , ( r n = p n − p n − ) , (100)can be seen as an indicator of a new trend. For independent Gaussian price variations r n , writing δ n +1 = (1 − λ ) δ n + (1 − λ ) r n +1 , one retrieves Eq. (6) for a discrete versionof an Ornstein-Uhlenbeck process, where ˆ x = µ (1 − λ ) /λ is related to the mean price ONTENTS µ , θ = − ln(1 − λ ), and σ = σ − λ ) √ θ √ − (1 − λ ) is proportional to the standarddeviation (volatility) σ of price variations.The indicator δ n can be used in both mean-reverting and trend following strategies.In the mean-reverting frame, if δ n exceeds a prescribed threshold L , this is a trigger to sellthe asset at its actual (high) price, in anticipation of its return to the expected (lower)level ¯ p n in near future. Similarly, the event δ n < − L triggers buying the asset. In theopposite trend following frame, the condition δ n > L is interpreted as the beginning ofa strong trend and thus the signal to buy the asset at its actual price, in anticipation ofits further growth (similarly for δ n < − L ). In other words, the same condition δ n > L (or δ n < − L ) can be interpreted differently depending on the empirical knowledgeon the asset behavior. Whatever the strategy is used, the statistics of crossing of theprescribed levels ± L is precisely the FET problem. Theoretical results in Sec. 2 can helpto characterize durations between buying and selling moments. In particular, the choiceof the threshold L is a compromise between execution of too frequent buying/sellingtransactions (i.e., higher transaction costs) at small L and missing intermediate trends(i.e. smaller profits) at large L . We also note that Ornstein-Uhlenbeck processes oftenappear in finance to model, e.g., interest rates (Vasicek model) and currency exchangerates [105, 106]. A general frame of using eigenfunctions for pricing options is discussedin [107]. The first exit time problem can be extended to time-evolving domains [108, 109, 110,111, 112]. For instance, one can investigate the first passage time of Brownian motionto a time-dependent barrier L ( t ), τ = inf { t > X ( t ) = L ( t ) } , or the first exit timefrom a symmetric “envelope” [ − L ( t ) , L ( t )], τ = inf { t > | X ( t ) | = L ( t ) } . Althoughthe survival probability S ( x , t ) satisfies the standard diffusion equation with Dirichletboundary condition, the boundary L ( t ) evolves with time. For a smooth L ( t ), setting S ( x , t ) = v ( z, t ) with a new space variable z = x /L ( t ) yields ∂v ( z, t ) ∂t = DL ( t ) ∂ z v ( z, t ) − L ′ ( t ) L ( t ) z ∂ z v ( z, t ) , (101)with Dirichlet boundary condition v ( ± , t ) = 0 at two fixed endpoints (here we focus onthe exit time). Setting a new time variable T = ln( L ( t ) /L (0)), the above equation canalso be written as the backward Fokker-Planck equation with time-dependent diffusioncoefficient D ( T ) = D L ′ ( t ) L ( t ) = D e − T L (0) L ′ ( L − ( L (0) e T )) and a centered harmonic potential: ∂v ( z, T ) ∂T = D ( T ) ∂ z v ( z, T ) − z ∂ z v ( z, T ) . (102)In higher dimensions, the second derivative ∂ z is simply replaced by the radial Laplaceoperator ∂ r + d − r ∂ r .In general, the above equation does not admit explicit solutions. A notableexception is the case of square-root boundaries which has been thoroughly investigated ONTENTS L ( t ) = p b ( t + t ) (with b > t > L ′ ( t ) L ( t ) = b so that D ( T ) is independent of T (or t ). In other words, oneretrieves the backward Fokker-Planck problem (17, 24) with ˆ x = 0, k/γ = 1, L = 1,and D replaced by D/b . Its exact solution is given by Eq. (73) for d -dimensional case: q ( z , T ) = 2 ∞ X n =0 e − T ν n M ( − ν n , d , bz D ) M ′ ( − ν n , d , b D ) , (103)where z = r /L (0) = r / √ bt denotes the rescaled starting point r , and ν n = α n / (4 κ )are zeros of M ( − ν, d , b D ) = 0. Changing back T to t , one gets p ( z , t ) = 1 t ∞ X n =0 (1 + t/t ) − ν n − M ( − ν n , d , bz D ) M ′ ( − ν n , d , b D ) , (104)This expression in a slightly different form was provided for d = 1 in [118]. Note alsothat the probability P { sup 0, (107)where two minima are located at − x and x (with x > x > k and k aretwo spring constants, and v = ( k x − k x ) is a constant ensuring the continuity ofthe potential at x = 0. The resulting Langevin equation remains linear, in contrastto other bistable potentials such as a quartic potential (e.g., V ( x ) = ax + bx + cx ).The diffusive dynamics under double-well potentials was thoroughly investigated by ONTENTS k i . However,neither Kummer, nor Tricomi function is appropriate to represent the solution in thiscase. In fact, the Kummer function M ( a, / , z ) rapidly grows at infinity, while theTricomi function U ( a, / , z ) behaves as √ π Γ( a +1 / − √ π Γ( a ) | z | + . . . for small z , i.e., itsderivative is discontinuous at 0. A convenient representation can still be obtained asa linear combination of two Kummer functions, in which the rapid growth of thesefunctions is compensated. This is precisely the case of parabolic cylinder functions D ν ( z ) and D ν ( − z ) which vanish as z → ∞ (resp., z → −∞ ) but rapidly grow as z → −∞ (resp. z → ∞ ) unless ν is a nonnegative integer [see Eqs. (B.3), (D.2), (D.3)].An eigenfunction can therefore be written as u ( x ) = c e κ ( x/x +1) / D ν (cid:16) −√ κ ( x/x + 1) (cid:17) , x ≤ c e κ ( x/x − / D ν (cid:16) √ κ ( x/x − (cid:17) , x ≥ 0, (108)where κ i = k i x i / (2 k B T ), ν i = λx i / (2 κ i D ), and λ , c , c are determined by normalizationand two interface conditions at x = 0. The continuity of the eigenfunction at x = 0 canbe satisfied by choosing c = β e κ / D ν ( −√ κ ) , c = β e κ / D ν ( −√ κ ) , where β is a normalization constant.The second interface condition is deduced from the orthogonality of eigenfunctionswith two weights w , from Eq. (11) for positive and negative semi-axes, w i ( x ) = exp( κ i [1 − ( x/x i ± ]) , where plus (resp., minus) corresponds to i = 1 (resp., i = 2). The orthogonality imposesthe interface condition Dw (0) u ′ (0 − ) − Dw (0) u ′ (0 + ) = 0 , (109)where the same diffusion coefficient D is assumed on both sides. Since w i (0) = 1, oneretrieves the standard flux continuity equation, u ′ (0 − ) = u ′ (0 + ), yielding an equationdetermining the eigenvalues λ : x D ν ( −√ κ ) h κ D ν ( −√ κ ) + √ κ D ν +1 ( −√ κ ) i + x D ν ( −√ κ ) h κ D ν ( −√ κ ) + √ κ D ν +1 ( −√ κ ) i = 0 , where we used the identity ∂∂z D ν ( z ) = z D ν ( z ) − D ν +1 ( z ), and λ appears in ν i = λx i / (2 κ i D ). The smallest eigenvalue λ = 0 corresponds to the steady state.The normalization constant β is found according to β − = e κ + κ x D ν ( −√ κ ) √ κ ∞ Z −√ κ dz D ν ( z ) + x D ν ( −√ κ ) √ κ ∞ Z −√ κ dz D ν ( z ) , (110) ONTENTS −4 −2 0 2 400.20.40.60.811.2 x p ( x ,t | x , ) (a) t = 0.1t = 1t = 5t = ∞ −4 −2 0 2 400.20.40.60.811.2 x p ( x ,t | x , ) (b) t = 0.1t = 1t = 5t = ∞ Figure 7. Evolution of the probability density p ( x, t | x , t ) for diffusion underquadratic double-well potential (sketched by black dotted line) with two minima at ± x = x = 1), and κ = 2, κ = 1. Dashed vertical line indicates the starting pointat t = 0: x = 2 (a) and x = − (b) . Symbols represent normalized histograms ofarrival positions obtained by Monte Carlo simulations of an adapted version of Eq. (6)(with time step δ = 10 − and 10 sample trajectories), while lines show the spectraldecomposition (18) with 50 terms. We set D = 1. in which both integrals can be partly computed by using the identity [130] ∞ Z dz D ν ( z ) = √ π / ψ ( − ν ) − ψ ( − ν )Γ( − ν ) , (111)where ψ ( z ) = Γ ′ ( z ) / Γ( z ) is the digamma function. The lowest eigenfunctioncorresponding to λ = 0, is constant, u ( x ) = β , with β − = √ π (cid:20) x e κ (1 + erf( √ κ )) √ κ + x e κ (1 + erf( √ κ )) √ κ (cid:21) . As a consequence, one retrieves the equilibrium Boltzmann-Gibbs distribution, p eq ( x ) = p ( x, ∞| x , 0) = u ( x ) u ( x ) w ( x ) = β w ( x ).Figure 7 illustrates the evolution of the probability density p ( x, t | x , t ) for diffusionunder quadratic double-well potential with two minima at ± x = x = 1) andtwo different dimensionless strengths: κ = 2 and κ = 1. One can see how the initialDirac distribution, concentrated at x = 2 (Fig. 7a) or x = − p eq ( x ) (shown by black solid line). Other diffusion characteristics can be deduced fromthe probability density. The spectral approach is a general tool for computing FETs and other first-passagequantities. We briefly mention four straightforward extensions.(i) In one dimension, one can easily derive the splitting probability H ( x ), i.e., theprobability to exit from one endpoint (e.g., x = L ) before the other ( x = − L ). Thesplitting probability is governed by the stationary equation L ∗ x H ( x ) = 0 so that H ( x ) ONTENTS α = 0. Two constants c and c are setby boundary conditions H ( L ) = 1 and H ( − L ) = 0, from which H ( x ) = erf( i √ κ ( x /L − ϕ )) + erf( i √ κ (1 + ϕ ))erf( i √ κ (1 − ϕ )) + erf( i √ κ (1 + ϕ )) , (112)where we used M (0 , b, z ) = 1 and M (1 / , / , z ) = √ π erf( i √ z )2 i √ z . Note that this expressioncan be recognized in Eq. (33) for the mean exit time.(ii) Dirichlet boundary conditions, q ( ± L, t ) = 0, were imposed on the FETprobability density at both endpoints in order to stop the process whenever it exits fromthe interval. One can consider other boundary value problems, e.g., with one reflectingendpoint or one/two semi-reflecting points. In this case, Dirichlet boundary condition atone or both endpoints is replaced by Neumann or Robin boundary conditions [46]. Forinstance, the condition ∂∂x q ( x , t ) = 0 at x = − L describes the reflecting barrier at − L .The Robin boundary condition, ∂∂x q ( x , t ) + hq ( x , t ) = 0, allows one to consider partialabsorptions/reflections for modeling various transport mechanisms on the boundaryand to switch continuously between Neumann (pure reflections) and Dirichlet (pureabsorptions) cases by varying h from 0 to infinity [131, 132, 133, 134, 135, 136, 137].The solution can be obtained in the same way.(iii) The first passage time to a single barrier can be deduced from the first exittime from an interval by sending one endpoint to infinity (see Appendix C).(iv) A straightforward extension of the spectral approach allows one to deduce FETsof continuous-time random walks (CTRW) [5, 6]. In this model, long stalling periodsbetween moves result in anomalous subdiffusion, when the mean-square displacementevolves sublinearly with time: h ( X ( t ) − X (0)) i ≃ D α t α , with the exponent 0 < α < D α . The same derivations can be formallyrepeated for the fractional Fokker-Planck equation that governs the survival probabilityof CTRWs. In practice, it is sufficient to replace s/D by s α /D α in the Laplace domainthat in time domain yields the replacement of exponential functions exp( − λ n t ) byMittag-Leffler functions E α ( − λ n D α t α /D ) in spectral decompositions such as Eq. (23)or similar. As expected for CTRWs, the mean exit time diverges due to long stallingperiods while the survival probability exhibits a power law decay t − α at long timesinstead of the exponential decay for normal diffusion. Conclusion We revised the classical problem of finding first exit times for harmonically trappedparticles. Although the explicit formulas for the moment-generating function h e − sτ i canbe found in standard textbooks (e.g., [39]), the computation of the probability densityand the survival probability through the inverse Laplace transform requires substantialanalysis of confluent hypergeometric functions. For didactic purpose, we reproduced themain derivation steps and resulting spectral decompositions that involve the eigenvaluesand eigenfunctions of the governing Fokker-Planck operator. We also provided explicit ONTENTS Acknowledgments The author acknowledges partial support by the grant ANR-13-JSV5-0006-01 of theFrench National Research Agency. Appendix A. Limit of large κ Appendix A.1. Mean exit time For large κ or ϕ , Eq. (33) is not appropriate for numerical computation of the meanexit time because integrals and error functions are exponentially large. Re-arrangingthese terms, one can rewrite Eq. (33) as h τ i z = L D √ π κ (cid:26) e − κϕ (1+ z ) − κ (1 − z ) D ( √ κ ( z − ϕ )) D ( √ κ (1+ ϕ )) e − κϕ D ( √ κ (1 − ϕ )) D ( √ κ (1+ ϕ )) (A.1) × √ κ (1+ ϕ ) Z √ κ ( ϕ − dz e z erfc( z ) − √ κ (1+ ϕ ) Z √ κ ( ϕ − z ) dz e z erfc( z ) (cid:27) , ONTENTS z = x /L , and D ( x ) is the Dawson function: D ( x ) = e − x x Z dt e t , (A.2)which is related to the error function of imaginary argument aserf( ix ) = 2 i √ π e x D ( x ) . (A.3)For large x , the Dawson function decays as D ( x ) ≃ x + 14 x + 38 x + . . . . (A.4)The relation (A.1) allows one to compute the mean exit time in the limit of large κ and/or ϕ . In fact, since the Dawson function vanishes for large argument, the ratioin front of the first integral in Eq. (A.1) becomes exponentially close to 1 so that h τ i z ≃ L D √ π κ √ κ ( ϕ − z ) Z √ κ ( ϕ − dz e z erfc( z ) . (A.5)Three situations can be considered separately.(i) If ϕ > 1, the upper and lower limits of the above integral are positive and large sothat h τ i z ≃ L D κ ln ϕ − z ϕ − κ ≫ , (A.6)where we used the asymptotic relation b Z a dz e z erfc( z ) ≃ ln( b/a ) √ π ( a, b ≫ . (A.7)Note that Eq. (A.6) is accurate already for ϕ & κ ≥ < ϕ < κ → ∞ , the lower limit goes to −∞ , and the integral exponentiallydiverges: h τ i ≃ L D √ π e κ (1 − ϕ ) κ / (1 − ϕ ) ( κ ≫ 1) (A.8)(here the starting point is set to 0, but the result holds for all z not too close to1). This relation is valid for any 0 < ϕ < 1. Setting formally ϕ = 0, one getsthe relation which is twice larger than the asymptotic Eq. (39) derived for ϕ = 0.The missing factor 2 can be retrieved from the ratio in front of the first integral inEq. (A.1). The difference between the cases ϕ = 0 and ϕ > ϕ > x = L is closer to the minimum position ˆ x than theleft endpoint x = − L . When κ is large, the probability of large deviations fromˆ x rapidly decays with the distance so that the probability of exiting through the ONTENTS ϕ = 0 (and thus ˆ x = 0), both endpoints are equivalentthat doubles the chances to exit and thus twice reduces the mean exit time.(iii) In the marginal case ϕ = 1, the integral in Eq. (A.5) grows logarithmically with κ .One can split the integral by an intermediate point ¯ z ≫ ¯ z Z dz e z erfc( z ) + √ κ (1 − z ) Z ¯ z dz e z erfc( z ) ≃ √ π ln √ κ (1 − z ) c (¯ z ) , where c (¯ z ) ≡ ¯ z exp −√ π ¯ z Z dz e z erfc( z ) −→ . . . . (¯ z → ∞ ) . We get therefore h τ i z ≃ L D κ ln √ κ (1 − z )0 . . . . ( κ ≫ . (A.9)This asymptotic relation is accurate starting from √ κ (1 − z ) & Appendix A.2. Eigenvalues (interior problem) For large κ , we search for positive solutions α n of the equation M ( − α n κ , b, κ ) = 0 in theform: α n / (4 κ ) = n − ε , where ε is a small parameter, and n = 0 , , , . . . One gets then0 = M ( − n + ε ; b ; κ ) ≃ n X j =0 ( − n )( − n + 1) . . . ( − n + j − κ j b ( j ) j !+ ε ( − n n ! ∞ X j = n +1 ( − n + j − κ j b ( j ) j ! + O ( ε ) , from which the small parameter ε can be determined as ε ≃ − S ( − n n ! S , (A.10)where S and S denote two above sums. The second sum can be written as S = ∞ X j = n +1 ( − n + j − κ j b ( j ) j ! = Γ( b ) ∞ X j =0 κ j + n +1 Γ( b + j + 1 + n )( j + 1) . . . ( j + 1 + n ) . This expression can be obtained by integrating n + 1 times the Mittag-Leffler function E ,b + n +1 ( κ ) which asymptotically behaves as E ,b + n +1 ( κ ) ≃ κ − b − n e κ (1 + O (1 /κ )) as κ ≫ 1. Since the integration does not change the leading term, one concludes that S ≃ Γ( b ) κ − b − n e κ (1 + O (1 /κ )) ( κ ≫ . Keeping the highest-order term in the first sum, S ≃ ( − n κ n /b ( n ) , one gets ε ≃ − κ b +2 n n !Γ( b + n ) e − κ , ONTENTS α n as κ ≫ α n ≃ κ h n + κ b +2 n e − κ n !Γ( b + n ) i ( n = 0 , , , . . . ) . (A.11)In particular, the smallest solution α exponentially decays with κ , α ≃ κ b Γ( b ) e − κ ( κ ≫ , (A.12)while the other eigenvalues grow linearly with κ : α n ≃ κn ( κ ≫ , n = 1 , , . . . ) , (A.13)and the first-order correction ε decays exponentially fast. This asymptotic behaviorcan be related to equidistant energy levels of a quantum harmonic oscillator (seeAppendix D).Since α rapidly vanishes, the first eigenfunction approaches the unity: M ( − α κ , d , κz ) → 1. As a consequence, the normalization constant is simply β ≈ κ d/ / Γ( d/ 2) so that w ≃ 1, because M (1 , b, z ) = Γ( b ) E ,b ( z ). Appendix A.3. Eigenvalues (exterior problem) For the exterior problem, we consider the asymptotic behavior of solutions of U ( − α κ , b, κ ) = 0 as κ → 0. For non-integer b , one can use Eq. (87) to write inthe lowest order in κ U (cid:16) − α κ , b, κ (cid:17) ≃ Γ(1 − b )Γ(1 − b − α κ ) + Γ( b − − α κ ) κ − b . (A.14)For b < κ − b is a small parameter so that the first term has to be small. Setting1 − b − α κ = − n + ε (with n = 0 , , , . . . ) one gets ε = ( − n − Γ( b − n ! Γ( b − − n )Γ(1 − b ) κ − b , (A.15)from which α n ≃ κ (cid:16) − b + n + ( − n Γ( b − n ! Γ( b − − n )Γ(1 − b ) κ − b + . . . (cid:17) . (A.16)In turn, if b > κ − b is a large parameter so that the second term has to be small.Setting − α κ = − n + ε , one gets ε = ( − n − Γ(1 − b ) n ! Γ(1 − b − n )Γ( b − κ b − , (A.17)from which α n ≃ κ (cid:16) n + ( − n Γ(1 − b ) n ! Γ(1 − b − n )Γ( b − κ b − + . . . (cid:17) . (A.18)For integer b , the analysis is more subtle and is beyond the scope of this paper. Wejust checked numerically that α ∝ κ b as κ → b = 1 and b = 2 that corresponds todimensions d = 2 and d = 4. ONTENTS Appendix B. Confluent hypergeometric functions For the sake of completeness, we summarize selected relations between special functionsthat are often used to describe first passage times of Ornstein-Uhlenbeck processes (see[68] for details). After that, we describe a rapidly converging representation of confluenthypergeometric functions. Appendix B.1. Relations The Kummer confluent hypergeometric function M ( a, b, z ) = F ( a ; b ; z ) defined in Eq.(42), satisfies the Kummer’s equation: zy ′′ + ( b − z ) y ′ − ay = 0 . (B.1)For b = 1 / 2, this equation is also related to the Weber’s equation y ′′ − ( z / c ) y = 0 , (B.2)which has two independent solutions: e − z / M ( c/ / , / , z / 2) (even) and ze − z / M ( c/ / , / , z / 2) (odd). These solutions are often expressed through theparabolic cylinder function D ν ( z ), which satisfies Eq. (B.2) with ν = − c − / D ν ( z ) = cos( πν )Γ( ν ) √ π − ν/ e − z / M (cid:18) − ν , , z (cid:19) (B.3)+ sin( πν )Γ( ν ) √ π − ( ν +1) / e − z / zM (cid:18) − ν , , z (cid:19) = 2 ν/ e − z / U (cid:18) − ν , , z (cid:19) (B.4)(the last relation is valid only for ℜ{ z } ≥ M ( a, b, z ) and U ( a, b, z ) are also related tothe Whittaker functions M a,b ( z ) and W a,b ( z ) [68] M a,b ( z ) = e − z/ z b +1 / M (1 / b − a, b, z ) ,W a,b ( z ) = e − z/ z b +1 / U (1 / b − a, b, z ) . The following relations help to analyze the Brownian motion limit [39]lim κ → M (cid:16) a κ , b + 1 , κx (cid:17) = 2 b Γ( b + 1) I b ( √ xa )( xa ) b/ , (B.5)lim κ → κ b Γ (cid:16) a κ (cid:17) U (cid:16) a κ , b + 1 , κx (cid:17) = 2 − b K b ( √ xa )( x/a ) b/ , (B.6)where I ν ( z ) and K ν ( z ) are the modified Bessel functions of the first and second kind,respectively.The asymptotic expansions for large | z | (and fixed a and b ) are [68] (Sec. 13.5): M ( a, b, z ) ≃ e z z a − b Γ( b )Γ( a ) N − X n =0 ( b − a ) ( n ) (1 − a ) ( n ) n ! z n + O ( | z | − N ) ! (B.7) ONTENTS e ± πia z − a Γ( b )Γ( b − a ) N − X n =0 ( a ) ( n ) (1 + a − b ) ( n ) n ! ( − z ) n + O ( | z | − N ) ! ,U ( a, b, z ) ≃ z − a N − X n =0 ( a ) ( n ) (1 + a − b ) ( n ) n ! ( − z ) n + O ( | z | − N ) ! ( | arg( z ) | < π/ , (B.8)where the upper [resp., lower] sign in the second line is taken if − π/ < arg( z ) < π/ − π/ < arg( z ) ≤ − π/ N , N , and N are truncation orders. Appendix B.2. ComputationSeries representations. The computation of the Kummer function M ( a, b, z ) by directseries summation in Eq. (42) is not convenient for large | a | . For this case, two equivalentrepresentations were proposed:(i) M ( a, b, z ) = Γ( b ) e z/ b − ∞ X n =0 A n z n J b − n ( p z (2 b − a ))( p z (2 b − a )) b − n , (B.9)where the coefficients A n are defined by A = 1 , A = 0 , A = b/ , nA n = ( n − b ) A n − + (2 a − b ) A n − (see [68], Sec. 13.3.7, and [138], Sec. 4.8). Note that the coefficients A n depend on a and grow with | a | .(ii) M ( a, b, z ) = Γ( b ) e z/ b − ∞ X n =0 p n ( b, z ) J b − n ( p z (2 b − a ))( p z (2 b − a )) b − n , (B.10)where p n ( b, z ) are the Buchholz polynomials in b and z (see [139], Sec. 7.4). Thesepolynomials are less explicit than the coefficients A n , but they are independent of a . As a consequence, this representation is particularly convenient for large | a | .The recurrence relations for the Buchholz polynomials were derived in [140]: p n ( b, z ) = ( iz ) n n ! [ n/ X k =0 (cid:18) n k (cid:19) f k ( b ) g n − k ( z ) , (B.11)where the polynomials f k ( b ) and g k ( z ) are defined recursively by f k ( b ) = − (cid:18) b − (cid:19) k − X j =0 (cid:18) k − j (cid:19) k − j | B k − j ) | k − j f j ( b ) , f ( b ) = 1 , (B.12) g k ( z ) = − iz [( k − / X j =0 (cid:18) k − j (cid:19) j +1 | B j +1) | j + 1 g k − j − ( z ) , g ( z ) = 1 , (B.13)and B j are the Bernoulli numbers. Using the recurrence relations between Besselfunctions, νx J ν ( x ) = J ν − ( x ) + J ν +1 ( x ), one can express J b − n ( x ) = P n (1 /x ) J b − ( x ) + Q n (1 /x ) J b ( x ) , where the polynomials P n ( z ) and Q n ( z ) are defined recursively P ( z ) = 1 , P ( z ) = 0 , P n +1 ( z ) = 2( b − n ) zP n ( z ) − P n − ( z ) ,Q ( z ) = 0 , Q ( z ) = 1 , Q n +1 ( z ) = 2( b − n ) zQ n ( z ) − Q n − ( z ) . ONTENTS x andmoderate z : M ( a, b, z ) = e z/ ∞ X n =0 p n ( b, z ) (cid:20) F b ( x ) P n (1 /x ) x n + G b ( x ) Q n (1 /x ) x n − (cid:21) , (B.14)where x = p z (2 b − a ), and F b ( x ) = Γ( b )2 b − x − b J b − ( x ) , G b ( x ) = Γ( b )2 b − x − b J b ( x ) . (B.15)In particular, for b = d/ 2, one has d F b ( x ) G b ( x )1 cos( x ) sin( x ) /x J ( x ) J ( x ) /x x ) /x (sin( x ) − x cos( x )) /x (B.16)The above recursive relations allow one to compute rapidly the polynomials p n ( b, z ), P n (1 /x ) and Q n (1 /x ). The series (B.14) can be truncated after 5-10 terms when | az | islarge enough, and z is not too large (see [140] for several examples).According to Eq. (87), one can apply this method to compute the Tricomiconfluent hypergeometric function U ( a, b, z ) for non-integer b . Other series expansionsfor U ( a, b, z ) are discussed in [141, 142]. For integer b , one can substitute b ε = b + ε intoEq. (87) and then take the limit ε → 0. This extension by continuity yields [68] U ( a, b, z ) = ( − b ( b − a − b + 1) (cid:26) M ( a, b, z ) ln z + ( b − a ) b − X k =0 ( a − b + 1) ( k ) z k − b +1 k ! (2 − b ) ( k ) + ∞ X k =0 a ( k ) z k k ! b ( k ) (cid:16) ψ ( a + k ) − ψ (1 + k ) − ψ ( b + k ) (cid:17)(cid:27) , b = 1 , , . . . , where ψ ( z ) = Γ ′ ( z ) / Γ( z ) is the digamma function, and the intermediate sum is omittedfor b = 1. In practice, one can apply the above numerical scheme to rapidly compute U ( a, b ε , z ) through M ( a, b ε , z ) with several non-integer b ε approaching the integer b , andthen to extrapolate them in the limit b ε → b .Taking the derivative of Eq. (B.10) with respect to a and using the relation J ′ ν ( x ) = νx J ν ( x ) − J ν +1 ( x ), one obtains ∂∂a M ( a, b, z ) = Γ( b ) e z/ b z ∞ X n =0 p n ( b, z ) J b + n ( p z (2 b − a ))( p z (2 b − a )) b + n (B.17)or, equivalently, ∂∂a M ( a, b, z ) = 2 ze z/ ∞ X n =0 p n ( b, z ) (cid:20) F b ( x ) P n +1 (1 /x ) x n +1 + G b ( x ) Q n +1 (1 /x ) x n (cid:21) , (B.18)with x = p z (2 b − a ). This expression allows one to rapidly compute the coefficients w n in the spectral representation of the survival probability. Similar relation can bederived for ∂∂a U ( a, b, z ) using Eq. (87) for non-integer b . Finally, one can also apply theseformulas for computing the parabolic cylinder function D ν ( z ) and its derivative ∂∂ν D ν ( z )which are used to characterize the first passage time to a single barrier (Appendix C). ONTENTS Integral representations. The above scheme is convenient for large | a | and moderate | z | . However, if | a | is moderate while | z | is large, the numerical convergence of theabove series is slowed down due to a rapid growth of Buchholz polynomials with z . Inaddition, the computation of the Tricomi function U ( a, b, z ) as a linear combination (87)of two large Kummer functions can result in significant round-off errors at large z . Inthis case, one can apply a different technique which relies on integral representations ofconfluent hypergeometric functions.For the Kummer function M ( a, b, z ), one can use the following integralrepresentation for ℜ{ b − a } > k M ( a, b, z ) = e z z − b Γ( b )Γ( b − a ) ∞ Z dt e − t t b − − a J b − (2 √ zt ) . (B.19)This representation is convenient for computing eigenvalues and eigenfunctions because a = − α / (4 κ ) < b = d/ > U ( a, b, z ) has an integral representation for positive a [68] U ( a, b, z ) = 1Γ( a ) ∞ Z dt e − zt t a − (1+ t ) b − a − ( ℜ{ a } > , ℜ{ z } > . (B.20)When a is negative, one can use the recurrence relation to increase a : U ( a − , b, z ) + ( b − a − z ) U ( a, b, z ) + a ( a + 1 − b ) U ( a + 1 , b, z ) = 0 . (B.21)Applying this relation repeatedly, one gets U ( a, b, z ) = p n ( a, b, z ) U ( a + n, b, z ) + q n ( a, b, z ) U ( a + n + 1 , b, z ) , (B.22)where the polynomials p n ( a, b, z ) and q n ( a, b, z ) can be rapidly computed throughrecurrence relations: p n ( a, b, z ) = q n − ( a, b, z ) − ( b − a + n ) − z ) p n − ( a, b, z ) , p = 1 ,q n ( a, b, z ) = − ( a + n )( a + n + 1 − b ) p n − ( a, b, z ) , q = 0 . Choosing n such that a + n > 0, one can express U ( a, b, z ) in terms of U ( a + n, b, z )and U ( a + n + 1 , b, z ) which are found by numerical integration of Eq. (B.20). If z is too large, it is convenient to divide each recurrence relation by z and to considerthem as polynomials of 1 /z . The resulting value can be compared with the asymptoticexpansion (B.8). Appendix C. First passage time to a single barrier The first passage times (one-barrier problem) for harmonically trapped particles haveattracted more attention than the first exit times (two-barrier problem) [29, 32, 33, 34,35, 36]. In general, the first passage time τ ℓ to a single barrier at ℓ > k See http://dlmf.nist.gov/13.16.E3 ONTENTS x lies on the right to ℓ (i.e., x > ℓ ), thisproblem is equivalent to the exterior problem to reach the interval [ − ℓ, ℓ ] from outside(see Sec. 2.9). In turn, if 0 < x < ℓ , the FPT to a single barrier ℓ can be deduced fromthe FET from the interval [ − a, ℓ ] in the limit a → ∞ .In order to illustrate this point, we focus on the moment-generating function ˜ q ( x , s )given by Eq. (60), with ℓ = L (1 − ϕ ) and a = L (1 + ϕ ). Setting 1 − ϕ = ε , we considerthe limit ε → 0, for which L = ℓ/ε → ∞ and 1 − ϕ = ε → a → ∞ while ℓ iskept fixed. The asymptotic behavior of functions m (1 , α,κ from Eq. (43) as ε → m (1) α,κ (1 − ϕ ) ≃ M (cid:16) a, , y (cid:17) ,m (2) α,κ (1 − ϕ ) ≃ εM (cid:16) a + 12 , , y (cid:17) ,m (1) α,κ ( − − ϕ ) ≃ Γ(1 / a ) (4 y/ε ) a − / e y/ε ,m (2) α,κ ( − − ϕ ) ≃ − / a + 1 / 2) (4 y/ε ) a − e y/ε , where we replaced α and κ by − Ds/L and kL / (2 Dγ ), introduced short notations a = sγ/ (2 k ) and y = kℓ / (2 Dγ ), and used the asymptotic relation (B.7) for the lasttwo functions. Substituting the above expressions into Eq. (60), one deduces in thelimit ε → q ( x , s ) = M ( a, , y ) + 2 Γ( a +1 / a ) √ y M ( a + 1 / , , y ) M ( a, , y ) + 2 Γ( a +1 / a ) √ yM ( a + 1 / , , y ) , (C.1)where y = kx / (2 Dγ ). Using Eq. (B.3), one can alternatively write the moment-generating function as˜ q ( x , s ) = exp (cid:18) k ( x − ℓ )4 Dγ (cid:19) D − sγ/k (cid:16) − x q kDγ (cid:17) D − sγ/k (cid:16) − ℓ q kDγ (cid:17) (0 ≤ x ≤ ℓ ) . (C.2)Note also that Eq. (88) for the exterior case x > ℓ can also be written in terms of theparabolic cylinder function D ν ( z ) according to Eq. (B.4):˜ q ( x , s ) = exp (cid:18) k ( x − ℓ )4 Dγ (cid:19) D − sγ/k (cid:16) x q kDγ (cid:17) D − sγ/k (cid:16) ℓ q kDγ (cid:17) ( x > ℓ ) , (C.3)in agreement with [39] (see also [42, 32]).The inverse Laplace transform yields the probability density p x,a ( t ) [42, 35] q ( x , t ) = − kγ exp (cid:18) k ( x − ℓ )4 Dγ (cid:19) ∞ X n =1 D ν n (cid:16) ± x q kDγ (cid:17) D ′ ν n (cid:16) ± ℓ q kDγ (cid:17) e − ν n kt/γ , (C.4)where 0 < ν < ... < ν n < ... are the zeros of the function D ν ( ± ℓ p k/ ( Dγ )), and D ′ ν n ( z )is the derivative of D ν ( z ) with respect to ν , evaluated at point ν = ν n [42] (p. 154). The ONTENTS x > ℓ and x < ℓ , respectively. Both D ν ( z ) and D ′ ν ( z ) can be rapidly evaluated by the numerical scheme presented in Appendix B.2.In the special case ℓ = 0, the FET probability density gets a simple explicit form: q ( x , t ) = x √ πD (cid:18) k/γ sinh( kt/γ ) (cid:19) / exp (cid:18) − kx Dγ e − kt/γ sinh( kt/γ ) + kt γ (cid:19) . (C.5)In the limit k → 0, one retrieves the classical formula for the FPT of Brownian motionat the origin q ( x , t ) = x √ πDt exp (cid:18) − x Dt (cid:19) ( k = 0) . (C.6) Appendix D. Quantum harmonic oscillator The eigenvalue problem (40) with b = 1 / m and frequency ω [120]. In fact, the eigenstates ψ n andenergies E n of the Hamiltonian H = ˆ p m + mω x m satisfy the time-independent Schr¨odingerequation h − ~ m ∂ x + mω x i ψ ( x ) = Eψ ( x ) , (D.1)where ˆ p = − i ~ ∂ x is the momentum operator, and ~ is the (reduced) Planck constant. Interms of the dimensionless coordinate z = x p mω/ ~ , the above Schr¨odinger equationis reduced to the Weber’s equation (B.2), with c = − E ~ ω . Setting ψ ( z ) = e − z / ˜ u ( z )yields ˜ u ′′ − z ˜ u ′ − ( c + 1 / u = 0, from which the rescaling u ( x ) = ˜ u ( p k/ ( Dγ )( x − ˆ x ))implies Eq. (40), with λ = − kγ ( c + 1 / λ = kγ ( E ~ ω − ).If no boundary condition is imposed, the non-normalized eigenstate is simply ψ ( x ) = D ν ( x p mω/ ~ ), where D ν ( z ) is the parabolic cylinder function (seeAppendix B.1), and ν = − c − / 2. One can check that D ν ( z ) ≃ e − z / z ν h − ν ( ν − z + O ( z − ) i ( z ≫ , (D.2) D ν ( z ) ≃ e − z / z ν h − ν ( ν − z + O ( z − ) i (D.3) − √ π Γ( − ν ) e πiν e z / z − ν − h ν + 1)( ν + 2)2 z + O ( z − ) i ( z ≪ − . In order to eliminate the unphysical rapid growth of the eigenstate as z → −∞ , oneneeds to impose ν = n with n = 0 , , , . . . to remove the last term, from which oneretrieves the quantized energies of the quantum harmonic oscillator: E n = ~ ω ( n + 1 / H n ( z ) ψ n ( x ) = 1 √ n n ! (cid:16) mωπ ~ (cid:17) / exp (cid:18) − mωx ~ (cid:19) H n (cid:18)r mω ~ x (cid:19) , where the usual normalization prefactor is included, and we used D n ( z ) =2 − n/ e − z / H n ( z/ √ ONTENTS L → ∞ , we retrieve the asymptotic behavior λ n ≃ kγ n or, equivalently, α n = L D λ n ≃ κn as κ → ∞ . Note that the prefactor 2 κ is twice smaller than that ofEq. (A.13) because the latter relation accounts only for symmetric eigenfunctions thatcontribute to the survival probability.Imposing Dirichlet boundary condition at x = ± L corresponds to setting infinitepotential outside the interval [ − L, L ] (and keeping the harmonic potential inside). Theeigenvalue problem for a quantum oscillator in such potential is equivalent to the analysisof the first exit time distribution in Sec. 2.6. [1] Redner S 2001 A Guide to First Passage Processes (Cambridge: Cambridge University press)[2] Metzler R, Oshanin G, and Redner S (Eds.) 2014 First-Passage Phenomena and TheirApplications (Singapore: World Scientific).[3] B´enichou O and Voituriez R 2014, From first-passage times of random walks in confinement togeometry-controlled kinetics, Phys. Rep. Phys. Rep. Phys. Rep. J. Phys. A: Math. Gen. R161-R208[7] Condamin S, B´enichou O, and Klafter J 2007, First-Passage Time Distributions for Subdiffusionin Confined Geometry, Phys. Rev. Lett. Phys. Rev.E Phys. Rev. E J. Chem. Phys. Phys. Rev. E Commun. Math. Phys. Europhys. Lett. EuroPhys. Lett. , 20008[15] Sanders LP and Ambj¨ornsson T 2012, First passage times for a tracer particle in single filediffusion and fractional Brownian motion, J. Chem. Phys. Phys. Rev. Lett. Phys. Rev. E J.Phys. A Phys. Rev. Lett. ONTENTS of surface-mediated diffusion in spherical domains, J. Stat. Phys. J. Stat. Phys. Phys. Rev. E J. Phys. A Rev.Mod. Phys. Nature , 77-80[26] Condamin S, Tejedor V, Voituriez R, B´enichou O, and Klafter J 2008, Probing microscopic originsof confined subdiffusion by first-passage observables, Proc. Nat. Acad. Sci. USA Phys. Rev. E J. Sound. Vibr. Phys. Rev. Ann. Math. Statist. J. Stat. Phys. J. Appl. Prob. Finance Stochast. Finance Stochast. Stoch. Models Quant.Finance Phys. Rev. E Phys. Rev. Lett. Handbook of Brownian Motion: Facts and Formulae (Birkhauser Verlag, Basel-Boston-Berlin)[40] Revuz D and Yor M 1999 Continuous martingales and Brownian motion , Third Ed. (Springer)[41] Itˆo K and McKean HP 1996 Diffusion Processes and Their Sample Paths Mathematical Methods For Financial Markets (Springer-Verlag)[43] Crank J 1975 The Mathematics of Diffusion , 2nd Ed. (Clarendon, Oxford)[44] Carslaw HS and Jaeger JC 1959 Conduction of Heat in Solids , 2nd Ed. (Clarendon, Oxford)[45] Hughes BD 1995 Random Walks and Random Environments (Clarendon Press, Oxford)[46] Gardiner C 2004 Handbook of Stochastic Methods for Physics, Chemistry and Natural Sciences Classical Potential Theory and Its Probabilistic Counterpart (Springer, Verlag,Berlin, Heidelberg, and New York)[48] van Kampen NG 2007 Stochastic Processes in Physics and Chemistry , 3ed. (Amsterdam: Elsevier)[49] Uhlenbeck GE and Ornstein LS 1930, On the theory of the Brownian motion, Phys. Rev. ONTENTS [50] Wang MC and Uhlenbeck GE 1945, On the Theory of the Brownian Motion II, Rev. Mod. Phys. Lectures inApplied Mathematics and Informatics ed. L. M. Ricciardi (Manchester University Press)[52] Coffey WT, Kalmykov YP, and Waldron JT 2004 The Langevin Equation: With Applications toStochastic Problems in Physics, Chemistry and Electrical Engineering The Fokker-Planck equation: methods of solution and applications , 3rd Ed. (Berlin:Springer)[54] Kolpas A, Moehlis J and Kevrekidis IG 2007, Coarse-grained analysis of stochasticity-inducedswitching between collective motion states, Proc. Nat. Acad. Sci. SIAMRev. Trans. Am. Math. Soc. Proc. 2nd Berkeley Symp. Math. Stat. Prob. edited by J. Neyman (Universityof California Press, Berkeley) 189-215[58] Freidlin M 1985 Functional Integration and Partial Differential Equations , Annals of MathematicsStudies (Princeton University Press, Princeton, New Jersey)[59] Simon B 1979 Functional integration and quantum physics (Academic press, London, New York,San Francisco)[60] Bass RF 1998 Diffusions and Elliptic Operators (Springer, New York)[61] Pontryagin L, Andronov A, and Witt A 1933, On the statistical treatment of dynamical systems, Zh. Eksp. Teor. Fiz. 172 [Reprinted in Noise in Nonlinear Dynamical Systems , 1989, ed. byF. Moss and P. V. E. McClintock (Cambridge University Press, Cambridge), Vol. 1, p. 329][62] Castro LB and de Castro AS 2013, Trapping of a particle in a short-range harmonic potentialwell, J. Math. Chem. J. Stat. Phys. J. Chem. Phys. J. Stat. Phys. Rev.Mod. Phys. Physica Handbook of Mathematical Functions (Dover Publisher, NewYork)[69] Rohrbach A, Tischer C, Neumayer D, Florin E-L, and Stelzer EHK 2004, Trapping and trackinga local probe with a photonic force microscope, Rev. Sci. Instrum. Phys. Rev. E Nat.Phys. Eur. Phys. J. E Phys. Rev. E ONTENTS [75] Ashkin A, Sch¨utze K, Dziedzic JM, Euteneuer U, and Schliwa M 1990, Force generation oforganelle transport measured in vivo by an infrared laser trap, Nature Science Trends Cell. Biol. Soft Matter Rev. Mod. Phys. Phys. Rev. E Phys. Rev. E Inv. Problems , 1053-1065 (2004)[83] Masson J-B, Casanova D, T¨urkcan S, Voisinne G, Popoff MR, Vergassola M, and AlexandrouA 2009, Inferring maps of forces inside cell membrane microdomains, Phys. Rev. Lett. ,048103[84] Kenwright DA, Harrison AW, Waigh TA, Woodman PG, and Allan VJ 2012, First-passage-probability analysis of active transport in live cells, Phys. Rev. E Phys. Rev. E Phys. Rev. Lett. Phys. Rev. Lett. Acta Phys. Pol. B Phys. Rev. Lett. Science Biophys. J. Biophys. J. Nature Phys. Rev. Lett. Ann. Rev. Biophys. Biomol. Struct. Surf. Sci. Rep. Time Series Analysis: Forecasting and Control Thirded. (Prentice-Hall)[98] Covel MW 2009 Trend Following (Updated Edition): Learn to Make Millions in Up or DownMarkets (Pearson Education, New Jersey)[99] Clenow AF 2013 Following the Trend: Diversified Managed Futures Trading (Wiley & Sons,Chichester UK)[100] Bouchaud J-P and Potters M 2003 Theory of Financial Risk and Derivative Pricing: From ONTENTS Statistical Physics to Risk Management (Cambridge University Press)[101] Chan LKC, Jegadeesh N, and Lakonishok J 1996, Momentum Strategies, J. Finance J. Finan. Econ. J. Finance Options, Futures and Other Derivatives (Upper Saddle River, NJ: Prentice Hall).[106] Vasicek O 1977, An Equilibrium Characterisation of the Term Structure, J. Financ. Econom. Operations Research J.Appl. Prob. J. Appl. Prob. Boundary crossing of Brownian motion (Vol. 40 of Lecture Notes in Statistics,Springer-Verlag, Berlin)[111] Wang L and P¨otzelberger K 1997, Boundary crossing probability for Brownian motion and generalboundaries, J. Appl. Prob. Ann.Appl. Probab. Proc 5th Berkeley Symp MathStatist Prob Ann. Math. Statist. Theory Prob. Appl. J. Appl. Prob. Adv. Appl. Prob. J. Appl. Prob. Commun.Stat. Theory Meth. Quantum Mechanics Introduction to Quantum Mechanics J. Stat. Phys. J. Stat. Phys. Z. Phys. B J.Stat. Phys. J. Stat. Phys. Z. Phys. B: Conden. Mat. Sov. Phys. JETP ONTENTS particle in a double well potential, Phys. Rev. E Table of Integrals, Series, and Products , (New York: AcademicPress)[131] Collins FC and Kimball GE 1949, Diffusion-controlled reaction rates, J. Colloid Sci. J. Chem. Phys. , 4009[133] Sano H and Tachiya M 1979, Partially diffusion-controlled recombination, J. Chem. Phys. Radiat. Phys. Chem. Phys. Rev.Lett. Focus on Probability Theory Ed. L. R. Velle, pp. 135-169 (Hauppauge NY: NovaScience Publishers)[137] Singer A, Schuss Z, Osipov A, and Holcman D 2008, Partially reflected diffusion, SIAM J. Appl.Math. The Special Functions and their Approximations (Academic Press, New York)[139] Buchholz H 1969 The Confluent Hypergeometric Function (Springer-Verlag, Heidelberg)[140] Abad J and Sesma J 1995, Computation of the Regular Confluent Hypergeometric Function, Mathematica J. J. Comput. Appl. Math. Numer. Math.41