Extreme Value Statistics and Arcsine Laws of Brownian Motion in the Presence of a Permeable Barrier

The Arcsine laws of Brownian motion are a collection of results describing three different statistical quantities of one-dimensional Brownian motion: the time at which the process reaches its maximum position, the total time the process spends in the positive half-space and the time at which the process crosses the origin for the last time. Remarkably the cumulative probabilities of these three observables all follows the same distribution, the Arcsine distribution. But in real systems, space is often heterogeneous, and these laws are likely to hold no longer. In this paper we explore such a scenario and study how the presence of a spatial heterogeneity alters these Arcsine laws. Specifically we consider the case of a thin permeable barrier, which is often employed to represent diffusion impeding heterogeneities in physical and biological systems such as multilayer electrodes, electrical gap junctions, cell membranes and fragmentation in the landscape for dispersing animals. Using the Feynman-Kac formalism and path decomposition techniques we are able to find the exact time-dependence of the probability distribution of the three statistical quantities of interest. We show that a permeable barrier has a large impact on these distributions at short times, but this impact is less influential as time becomes long. In particular, the presence of a barrier means that the three distributions are no longer identical with symmetry about their means being broken. We also study a closely related statistical quantity, namely, the distribution of the maximum displacement of a Brownian particle and show that it deviates significantly from the usual half-Gaussian form.


Introduction
Diffusion is ubiquitous, appearing as a transport mechanism in many physical, chemical and biological systems. Often these systems are littered with spatial heterogeneities which impede or hamper the diffusive motion. One of the most common forms of these heterogeneities are permeable barriers, such that diffusive particles either pass through or are reflected upon interaction with the barrier. These barriers appear at various scales inhibiting diffusive movement in many physical and biological contexts, from multilayer electrodes [1,2,3] and water transport in rock pores [4] to drug delivery systems [5,6]. Permeable material can be found in many examples in cell biology whose function is to regulate the flux of biochemicals between spatial regions [7], such as the bilayer plasma membrane of eukaryotes [8,9,10] and the electrical gap junctions in neurons [11,12]. As well as being prominent in magnetic imaging techniques of water molecules diffusing through cellular compartments [13,14]. Permeable substances also present themselves at larger scales such as heterogeneous landscapes, e.g. habitat type [15,16] or the presence of roads and fences [17], affecting the dispersal of animal movement at ecological scales. All these examples make it apparent that it is necessary to build a mathematical understanding of how the diffusive movement statistics is affected by the presence of a permeable barrier. Here we do so by investigating how the extreme value statistics (EVS) and the closely related Arcsine laws of Brownian motion (BM) change with permeable barriers.
The celebrated Arcsine laws of BM are a landmark result from Lévy [18] describing the probability distribution of three observables of a BM trajectory, x(τ ), starting at the origin, over a time interval τ ∈ [0, t]: (1.) t m (t), the time at which the trajectory reaches its maximum value, x(t m ) = M , (2.) t r (t), the total time the trajectory spends in the positive region, x(τ ) > 0, (3.) t ℓ (t), the time at which the trajectory crosses the origin for the last time. These quantities are displayed for a sample BM trajectory in figure 1. The remarkable feature of these quantities, t m (t), t r (t) and t ℓ (t) is that they all have the same cumulative distribution function [18,19], for i ∈ {m, r, ℓ}. Then the probability densities of these quantities is given by, which displays the counterintuitive property that the density diverges at t i = 0 and t i = t, which means the value of t i is more likely to be either extreme, with the mean, t/2, being the minimum. The first Arcsine law in particular has gained a lot of interest due to its prominent role in the field of extreme value statistics [20,21], where one is often interested in the maximum position, M (t), as well as the time of this event, t m (t). In this case the joint density is sought, which in the BM case is given by [22], where D is the diffusion constant. Marginalization over M of equation (2) recovers equation (1) whereas over t m , the following half-Gaussian distribution is found [21] P(M, t|0) = e − M 2 4Dt √ πDt .
(3) In recent literature there has been a large effort to extend these EVS and Arcsine laws to when the underlying stochastic process is not the simple BM. EVS studies include, various generalizations of BM [22,23,24,25,26,27,28,29,30,31,32,33], as well as other stochastic processes such as Bessel [22], Lévy flights [34,35], random acceleration [36,37] and run-and-tumble [38,39,40,41]. More recently, the time between the maximum and minimum of a stochastic process has been studied as well [42,43]. In addition, the other two Arcsine laws have been studied, together or separately, in numerous contexts, such as constrained and resetting BM [24,44], BM in random mediums [45], bounded regions [46,47,48], and with spatial and temporal heterogeneity [29,49], as well as other stochastic processes such as continuous time random walks [50], random acceleration [51], run-and-tumble [38], fractional BM [33], subdiffusion [52,53,54] and random systems with quenched disorder [55]. Despite this vast literature no such study on EVS and Arcsine laws of Brownian motion in the presence of permeable barriers has appeared. Here we provide the first such study.
The paper is structured as follows. In section 2 we introduce the key equations in modelling diffusion in the presence of permeable barriers. In section 3 we derive the joint density of M (t) and t m (t) and study the marginal densities. Section 4 and 5 are devoted to the density of t r (t) and t ℓ (t), respectively, while we summarize our work in section 6.

Diffusion through Permeable Barriers
We consider a Brownian particle, x(t), undergoing diffusion in one dimension in the presence of a permeable barrier at the origin. Classically, this has been modelled by the diffusion equation (DE) along with the following permeable barrier condition (PBC) [56,57], where P (x, t) represents the probability density of the position of the Brownian particle at time t, with J(x, t) = −D∂ x P (x, t) being the probability current and the parameter κ representing the permeability of the barrier, respectively, with κ → 0 representing an impenetrable barrier (reflecting boundary) and no barrier when κ → ∞. The notation 0 ± indicates the right or left-hand side of the permeable barrier, respectively. As the DE with condition (4) does not lend itself to study quantities other than P (x, t) the present authors have introduced a new fundamental equation [58] with which one reformulates the problem in terms of a modified DE that accounts for the PBC in an inhomogeneous term (taken here at the origin for simplicity), where δ ′ (x) represents the derivative of the Dirac-δ function. The convenience of equation (5) is that for any initial condition localized at x 0 , P (x, 0) = δ(x − x 0 ) can be solved exactly in the Laplace domain as [58] P in terms of the Green's function of Brownian motion in the absence of a permeable barrier, Generalizations to the case when an external potential is present are also possible with G 0 (x, t|x 0 ) becoming the Green's function of the associated Smoluchowski equation. In equation (6) we have used the Laplace variable ϵ to indicate that for an arbitrary time-dependent function, f (t), has its Laplace transformed expression given by f (ϵ) = ∞ 0 e −ϵt f (t)dt. The notation P (x, t|x 0 ) indicates the localized initial condition at x 0 and J 0 (x, t|x 0 ) = −D∂ x G 0 (x, t|x 0 ) is defined as the free probability current where ∂ x 0 G 0 (x, ϵ|0) implies ∂ ∂x 0 G 0 (x, ϵ|x 0 )| x 0 =0 . Further convenience in employing the formalism associated with equation (5) is due to the ability to write an equivalent Fokker-Planck (FP) representation, namely the homogeneous (Itô) FP equation [58], is the Fokker-Planck operator, with the spatially dependent drift and diffusion coefficients given by A(x) = −(D 2 /κ)δ ′ (x) and B(x) = D − (D 2 /κ)δ(x). As will become apparent later, we are mainly interested in the backward FP (Kolmogorov) equation in terms of the initial variables, t 0 and x 0 , and exploiting the time homogeneity of the process, we have [59] As was shown in Ref. [58], L is in fact self-adjoint, such that the backward FP equation, is given by which implies that we can solve equation (8) through the same procedure for which we obtain equation (6). We study the time-dependent joint distribution, P(M, t m , t|x 0 ), of the maximum position M = x(t m ) and the time t m at which this occurs for a Brownian particle in the presence of a permeable barrier at the origin. Since we consider the particle starting from x 0 = 0 the presence of the permeable barrier leads to two different joint densities P(M, t m , t|0 + ) and P(M, t m , t|0 − ).

Extreme value
To find these joint densities we proceed by using a path decomposition approach [24,36,29,27], where we exploit the Markovian nature of the process to split the trajectories of the Brownian particle into two parts. The first part is {x(τ ) : τ ∈ [0, t m ]} and the second part is {x(τ ) : τ ∈ [t m , t]}, see figure 1. Clearly the first part of the trajectory can be expressed as the probability of reaching M for the first time at t m , which is the first-passage time distribution F(M, t m |0 ± ). The second part is the probability that the particle starting at M does not reach this position again in the remaining time, this is the survival probability S(M, t − t m |M ). As S(M, t − t m |M ) = 0, we remedy this by using the quantity S(M + ε, t − t m |M ) and taking ε → 0 + [24,29].
We now proceed to find the two quantities, F(x, t|x 0 ) and S(x, t|x 0 ). This calculation was performed in detail in Ref. [58], where using equation (8) a backward equation was found for S(x, t|x 0 ), which was used to find F(x, t|x 0 ) in the Laplace domain through F(x, ϵ|x 0 ) = 1 − ϵ S(x, ϵ|x 0 ): We now proceed to calculate the joint distributions for x 0 = 0 ± .
Interestingly, if we take the small t limit of P(M, t m , t|0 + ), which corresponds to ϵ, p → ∞ in the Laplace domain for equation (12), we obtain to leading order equation (13). This implies that for t → 0 the permeable barrier acts as if it is fully reflecting. This can be understood by considering for very small times the particle does not have enough time to interact with the barrier and pass through, meaning that essentially no trajectories reach the other side of the barrier.
Although the double inverse Laplace transform of equation (12) for arbitrary κ looks highly non-trivial, significant ground can be made (see Appendix B). We find P(M, t m , t|0 + ) in terms of two different scaling functions such that, where, and with From equation (15) one can see we no longer have the transformation symmetry of t m → t − t m , meaning the symmetry about t/2 that one observes in the barrier free case, equation (2), is broken.
For x 0 = 0 − we have the added complexity due to the chance that the particle will not cross the barrier throughout the whole time period, leading to the maximum occurring at the initial position x 0 . The probability of this occurring is given by (from equation (9)) where This means that the maximum position is the origin and it occurs at time t m = 0. Then from equation (10) we have, (20) Similar to equation (11) we perform a double Laplace transform and obtain Expanding to first order in ε and ensuring normalization we once again find N (ε) = ε, leading to As in equation (12) one can see in the limit κ → ∞ equation (22) reduces to the barrier free case, equation (2). In the fully reflecting limit κ → 0 one recovers P(M, t m , t|0 − ) = δ(M )δ(t m ), since no particles pass through, meaning the maximum position will be at the origin being reached instantly. As in the x 0 = 0 + case, we see from equation (22) that in short time limit one recovers the fully reflecting limit up to leading order in t m . However, using Q m (M, p, ϵ|0 + ) for p, ϵ → ∞, equation (13), we find a better approximation which after using L −1 , equation (23) becomes where (17) and where h(y, z) is defined in equation (18). Again we see the t m → t − t m symmetry breaking as in the x 0 = 0 + case.

Marginal Density
To find the marginal density P(M, t|0 ± ) one integrates over t m , i.e. p → 0 in equations (12) and (22) and after taking the double Laplace inverse (see Appendix C) we get where and with j ± (y, z) defined in equations (C.3) and (C.5) and h(y, z) defined in equation (18). We have verified that in the barrier free case (κ → ∞) one recovers equation (3). From equations (27), (28) and (29), we find that for t >> D/κ 2 with M ∼ D/κ, the integrands in (28) and (29) are dominated by small z, so we expand j ± (y, z)h 2 (y, z) to first order in z and then compute the integral to obtain, and Equations (30) and (31) should be compared with the barrier free case, We plot equation (27) in figure 2 for the two initial conditions, x 0 = 0 ± whilst varying the dimensionless parameter, κ 2 t/D, and compare to stochastic simulations. One can see for certain parameter values, namely for small permeabilities, the distributions become non-monotonic, and we have a bi-modal shape. In the x 0 = 0 + case this feature can be explained by considering for small κ, the particle is unlikely to pass through the barrier and thus being more likely to move away from the barrier leading to a maximum occurring for M > 0. Whereas the peak at the origin is caused by the particle managing to pass through the barrier but never returning. The presence of a local minimum near M = 0 is thus caused by the less likely scenario in which the particle stays near the barrier, constantly interacting with it, but without venturing too far from the origin. The case x 0 = 0 − can be explained similarly except there is a non-zero probability that the particle never crosses the barrier leading to the global maximum always occurring at the origin (the Dirac-δ function).

Marginal Density
Here we consider the marginal density P(t m , t|0 ± ) = ∞ 0 P(M, t m , t|0 ± )dM , where from equations (15) and (25), we obtain the scaling relation, where The integral in equation (33) is hard to compute, but we can study the asymptotic forms of F ± m (τ 1 , τ 2 ). Firstly, we study the short time asymptotics, t << D/κ 2 , which can be approximated by the marginal over M of equations (14) and (24). By using the series form of the Jacobi theta functions and integrating over we have Now we study the long time asymptotics of P(t m , t|0 ± ), which corresponds to t >> D/κ 2 . We find it more convenient to do this in the Laplace domain, where we use equations (12) and (22). There are two regimes which give the full picture of the asymptotics for t → ∞: keeping t m finite in the first regime and keeping t − t m finite in the second. This corresponds to ϵ → 0 with p >> ϵ and p ∼ ϵ in the Laplace domain, respectively. Expanding around ϵ, p → 0 to leading order in equations (12) and (22) and integrating over M one can see that we obtain, Q m (p, ϵ|0 ± ) = 1/ ϵ(p + ϵ), which after double inverse Laplace transforming recovers the Arcsine distribution, (1). For the case p >> ϵ, let us expand Q m (p, ϵ|0 + ) around ϵ → 0 keeping p finite, then we find which gives after performing the integral and inverse Laplace transforming with respect to ϵ, Since equation (37) is very difficult invert, we can investigate the asymptotic dependence for t m → 0. This corresponds to p → ∞, thus by expanding Eq. (37) to leading order for p → ∞ and using arctan(1) = π/4, we find Q m (p, t|0 + ) ∼ √ π 2 √ pt , and performing the Laplace inversion with respect to p we get, For the case x 0 = 0 − , using equation (22), we find, It is clear that equations (38) and (39) deviate from the Arcsine law, ∼ 1/π(t m t) −1/2 , showing the permeable barrier has an influence for t → ∞ when t m is finite. And for the x 0 = 0 − case the singularity at t m → 0 is provided by the Dirac-δ function. We plot equation (32) in figure 3 to show the excellent match to simulations. One can see that for the different initial conditions x 0 = 0 ± , P(t m , t|0 ± ) has two different shapes. When x 0 = 0 + and for small permeabilities one can see the minimum of the distribution is just after t m = 0, which is due to the small likelihood of the particle reaching its maximum and then passing through the barrier and never returning. The peak at t m = 0 indicates the particle instantly crosses the barrier and never returns. For x 0 = 0 − and for small permeabilities the distribution is strongly weighted by the Dirac-δ function indicating the particle never crosses the barrier. However, if the barrier is crossed, it is very unlikely for the particle to cross again meaning it is more likely to reach the maximum at t m = t.

Residence Time in Positive Half-Space t r
Here we study the second Arcsine law, namely the residence time the particle, starting at the origin, spends in the region x > 0 when a permeable barrier is placed at the origin. Once again, due to the presence of this permeable barrier, we have two distinct initial positions, x 0 = 0 + and x 0 = 0 − . To proceed we utilize the fact that the residence time can be written as the functional, t r (t) = t 0 Θ(x(t ′ ))dt ′ , where Θ(z) is the Heaviside step function. From Feynman-Kac theory [62,63], we know that Q r (p, t|x 0 ) = ∞ 0 e −ptr P(t r , t|x 0 )dt r = e −ptr(t) x 0 = x(0) satisfies the following backward Feynman-Kac equation, where we use the backward Fokker-Planck operator in equation (7). Using the selfadjoint nature of this backward Fokker-Planck operator, we may write equation (40) as with the initial condition Q r (p, 0|x 0 ) = 1. In the Laplace domain, t → ϵ, the solution to (41) is given by (see Appendix A), where G(x, p, t|x 0 ) is the Green's function of the barrier free Feynman-Kac equation, i.e. (43) can be found by solving in the Laplace domain,

Time of the Last Crossing of the Origin t ℓ
Finally, we consider the third Arcsine law, the probability density, P(t ℓ , t|0 ± ), of the last time the Brownian particle crosses the origin, which in our case is equivalent to passing through the permeable barrier. Unlike the previous two Arcsine laws, the probability distribution is the same for either x 0 = 0 + or x 0 = 0 − , as the crossing process is symmetric about the origin. Thus, for simplicity we take x 0 = 0 + . To find P(t ℓ , t|0 + ) we exploit the Markovian nature of the process and use a path decomposition approach similar to the method used in section 3. Again we split the trajectory into two parts, {x(τ ) : τ ∈ [0, t ℓ ]} and {x(τ ) : τ ∈ [t ℓ , t]}, where the first part is the trajectory that crosses the origin at t ℓ , and the second part is that the trajectory does not reach the origin for the remaining time t − t ℓ . However, the presence of a permeable barrier at the origin adds further complexity, due to the probability of reaching 0 + and 0 − being different. Therefore, we write P(t ℓ , t|0 + ) as, where the sum in equation (56) comes from the last passage being an up crossing or a down crossing. We take the limit ε → 0 + to find the double Laplace transform of P(t ℓ , t|0 + ), i.e. Q ℓ (p, ϵ|0 + ) = ∞ 0 dte −ϵt ∞ 0 dt ℓ e −pt ℓ P(t ℓ , t|0 + ), such that where P (x, ϵ|x 0 ) is given by equation (6) and S(±ε, ϵ|0 ± ) is found to first order in ε as previously. To find the normalization factor N (ϵ), we use the normalization condition, which indicates that due to the presence of the permeable barrier, there is a probability that the particle will not cross the origin in a time t. After substituting these expressions into (57) and using (58) we find N (ε) = ε/D, giving using standard inverse Laplace transform relations [60], we obtain where for f + (τ 1 , τ 2 ) defined in equation (50). The striking feature here is that the scaling function that fully describes P(t ℓ , t|0 + ) is also part of the scaling function describing the residence time density, P(t r , t|0 + ). Comparing the short time limit, t << D/κ 2 , and the long time limit, t >> D/κ 2 for t ℓ << t , to that of P(t r , t|0 + ), i.e. equations (53) and (55), we see that the peaks at the origin of both distributions have the same time dependence. This feature can be understood by considering that for instantaneous last crossing times, t ℓ << t, this crossing event is almost certainly going to be the first and last crossing, which corresponds to the residence time being equivalent to the last crossing time t r = t ℓ . This breaks down as t ℓ gets larger, causing the different dependencies of the distributions. In this case equation (63) is not valid for t − t ℓ → 0 as f + (τ 1 , 0) = 1/ √ πτ 1 .
We plot equation (60) in figure 5 to compare with stochastic simulations, where we see an excellent match. In comparing different κ values the main feature is the presence of a sharper peak at t ℓ = 0 the smaller the permeability, indicating that once the particle crosses the barrier it is unlikely to do so again. Differently from the Arcsine distribution (1), there is no divergence at t ℓ = t, which is due to there being a non-zero probability of not crossing the barrier for every interaction. We note that instead if one is interested in the distribution of the particle returning to the same side of the barrier for the last time, one would recover the Arcsine distribution.

Conclusion
In summary, we have investigated the extreme-value statistics and Arcsine laws of Brownian motion in the presence of a permeable barrier at the origin, using an inhomogeneous diffusion equation which accounts for the presence of a permeable barrier. The presence of this barrier requires considering two initial positions, the right and left-hand side of the barrier, i.e. x 0 = 0 + and x 0 = 0 − , respectively.
Firstly, using a path-decomposition technique we have obtained the joint density of the maximum displacement, M (t), and the time to reach it, t m (t), for the two initial conditions, P(M, t m , t|0 ± ). We have found that this quantity can be represented by the multiplication of two different scaling functions, indicating the distribution is no longer symmetric under a t m → t − t m transformation. For the x 0 = 0 − case a Dirac-δ centered at t m = 0 is also present due to the probability of not passing through the barrier in finite time. At short times, t << D/κ 2 , P(M, t m , t|0 + ) can be asymptotically approximated by the distribution for when the barrier is fully reflecting, κ = 0, which is given in terms of Jacobi Theta functions. This approximation is valid due to a very small number of particles mananging to go through the barrier for short times. For P(M, t m , t|0 − ) we find a stronger approximation than just the reflecting barrier distribution, which is δ(M )δ(t m ), which accounts for the very few particles that do cross the barrier.
We have also investigated the respective marginal distributions, P(M, t|0 ± ) and P(t m , t|0 ± ). The presence of the barrier has a large impact on the monotonicity of the distribution, P(M, t|0 ± ), where for certain permeabilities a maximum appears away from the origin. At long times, t >> D/κ 2 , the distribution is still dependent on κ. For P(t m , t|0 ± ) the distribution remains asymmetric for long times, such that for large t m the usual Arcsine distribution is recovered, whereas we get a different dependence for t m → 0.
For the rest of the paper we have investigated the other two Arcsine laws, namely the distributions of the residence time, t r (t), and the last crossing of the origin, t ℓ (t). Using Feynman-Kac theory we have calculated P(t r , t|0 ± ) analytically and have found the dependence in terms of a scaling function. For x 0 = 0 + and x 0 = 0 − we have a Dirac-δ located at t r = t and t r = 0, respectively, due to the particle never crossing the barrier. In the short time limit we have found that the scaling functions for x 0 = 0 ± are equivalent under the transformation t r → t − t r , and in the long time limit we recover the Arcsine distribution.
Finally, we have studied P(t ℓ , t|0 ± ), which is equivalent for either initial condition x 0 = 0 + or x 0 = 0 − , since t ℓ (t) is a crossing event. Taking into account up and down crossings we have found P(t ℓ , t|0 ± ) in terms of known functions. Interestingly the scaling function that describes this distribution happens to be part of the scaling function which describes P(t r , t|0 + ), where the peak at the origin of both distributions have the same dependence, because for t ℓ << t the first crossing very likely corresponds to the last crossing, which implies t ℓ = t r . As t ℓ becomes larger this is no longer the case, and we do not observe a divergence as t ℓ → t, due to the non-zero probability of not crossing the barrier for any interaction.
Possible extensions of this study would include the analysis of how a permeable barrier affects the time, T , between the maximum and minimum of the process [42,43]. It would be interesting to see if the presence of a barrier breaks the symmetry around T = 0 and whether the initial position i.e. x 0 = 0 ± leads to differences in the respective distributions. Another interesting avenue to explore is changing the underlying Brownian motion to a different stochastic process, such as anomalous subdiffusion (where an equation akin to (5) has been found for this case [65]), and what is the impact of a permeable barrier to the unperturbed statistics. Futhermore, in a recent study the joint distribution of the first passage time to a target and the number of distinct sites visited when the target is reached of a random walk has been obtained [66]. It would be of interest to quantify how the presence of a permeable barrier affects this distribution, where specifically one would be able to glean the impact of a spatial heterogeneity on the kinetics and geometry of space exploration. tional facilities of the Advanced Computing Research Centre, University of Bristolhttp://www.bristol.ac.uk/acrc/.

Appendix A. Solution of the Feynman-Kac Equation
Here we show how the solution of equation (41), Q r (p, t|x 0 ), can be represented in terms the Green's function, G(x, p, t|x 0 ), of the barrier free Feynman-Kac equation, (43). If we take the last term on the right-hand side (RHS) of equation (41) as an inhomogeneous term, we may construct the solution as follows, Using Q r (p, 0|y) = 1 with Q(p, t|x 0 ) = ∞ −∞ G(y, p, t|x 0 ) and Laplace transforming, t → ϵ, we have Then by taking the derivative of both sides of equation (A.2) with respect to x 0 and setting x 0 = 0, we find To perform the Laplace inversion of Q m (M, p, ϵ|0 + ), we write equation (12) as then P(M, t m , t|0 + ) is given by the two Bromwich integrals, where γ 1 and γ 2 are greater than the real part of all singularities of G + κ D M, s 1 and H κ D M, s 2 , respectively and G + (y, s 1 ) and H(y, s 2 ) are given by, and H(y, To find P(M, t m , t|0 + ) we require the following Laplace inversions L −1 s 1 →τ 1 G + (y, s 1 ) and L −1 s 2 →τ 2 H(y, s 2 ) , then P(M, t m , t|0 − ) only requires the Laplace inversion, From the definition of G + (y, s 1 ) in equation (B.3) we see that G + (y, s 1 ) has no poles but has a branch point at s 1 = 0. Thus, by taking the branch cut to be the negative real axis, then G + (y, s 1 ) is analytic inside the contour, C, in figure B1, meaning that C e s 1 τ 1 G + (y, s 1 )ds 1 = 0. Then by taking R → ∞ and α → 0 for the contour C the Laplace inversion is given by, because in the limit R → ∞ the contributions from C 2 and C 6 vanish, and for α → 0 the contribution from C 4 is zero. Therefore, all we require is finding the integrals C 3 and C 5 . For C 3 we let s 1 = ze iπ , giving and for C 5 we let s 1 = ze −iπ , leading to Taking R → ∞ and α → 0, substituting back into equation (B.5) and converting the complex exponentials to trigonometric functions, we obtain G + (y, τ 1 ) as defined in equation (16). Similarly, since G − (y, s 1 ) = ( √ s 1 + 1) −1 G + (y, s 1 ), one can see that by altering the above calculations it leads to the definition of G − (y, τ 1 ) in equation (26). because in the limit R → ∞ and α → 0 the contributions from C 2 , C 6 and C 4 vanish, and we have After taking R → ∞, α → 0 and substituting these expressions into equation (B.8) one can see that we recover H(y, τ 2 ) in equation (17).