Extreme value statistics and Arcsine laws of Brownian motion in the presence of a permeable barrier

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–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–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(\tau)$, starting at the origin, over a time interval $\tau\in[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(\tau)\gt0$, (3.) $t_\ell(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_\ell(t)$ is that they all have the same cumulative distribution function [18, 19],

$$\begin{equation*} \mathbb{P}[t_i(t)\leqslant \tau|x(0) = 0] = \frac{2}{\pi} \arcsin\left(\sqrt{\frac{\tau}{t}} \right), \end{equation*}$$

for $i\in\{m,r,\ell\}$. Then the probability densities of these quantities is given by,

$$\begin{equation} \mathcal{P}(t_i,t|0) = \frac{1}{\pi \sqrt{t_i(t-t_i)}}, \end{equation} \tag{ 1 }$$

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.

Zoom In Zoom Out Reset image size

The first Arcsine law in particular has gained a lot of interest due to its prominent role in the field of EVS [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],

$$\begin{equation} \mathcal{P}(M,t_m,t|0) = \frac{M}{2 \pi D t_m^{3/2}\sqrt{t-t_m}}e^{\frac{-M^2}{4 D t_m}}, \end{equation} \tag{ 2 }$$

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]

$$\begin{equation} \mathcal{P}(M,t|0) = \frac{e^{-\frac{M^2}{4D t}}}{\sqrt{\pi D t}}. \end{equation} \tag{ 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–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–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–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–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. Sections 4 and 5 are devoted to the density of $t_r(t)$ and $t_\ell(t)$, respectively, while we summarize our work in section 6.

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],

$$\begin{equation} J(0^\pm,t) = \kappa[P(0^-,t)-P(0^+,t)], \end{equation} \tag{ 4 }$$

where $P(x,t)$ represents the probability density of the position of the Brownian particle at time t, with $J(x,t) = -D \partial_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 $\kappa\to \infty$. The notation $0^\pm$ 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),

$$\begin{equation} \frac{\partial P(x,t)}{\partial t} = D \frac{\partial^2 P(x,t)}{\partial x^2} + \frac{D}{\kappa}\delta^{^{\prime}}(x)J(0,t), \end{equation} \tag{ 5 }$$

where $\delta^{^{\prime}}(x)$ represents the derivative of the Dirac-δ function. The convenience of equation (5) is that for any initial condition localized at x₀, $P(x,0) = \delta(x-x_0)$ can be solved exactly in the Laplace domain as [58]

$$\begin{equation} \widetilde{P}(x,\epsilon|x_0) = \widetilde{G}_0(x,\epsilon|x_0) -\partial_{x_0}\widetilde{G}_0(x,\epsilon|0)\frac{\widetilde{J}_0(0,\epsilon|x_0)}{\frac{\kappa}{D}+ \partial_{x_0}\widetilde{J}_0(0,\epsilon|0)}, \end{equation} \tag{ 6 }$$

in terms of the Green's function of Brownian motion in the absence of a permeable barrier, $G_0(x,t|x_0) = \exp\left\{-(x-x_0)^2/4Dt \right\}/\sqrt{4 \pi Dt}$. 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 $\widetilde{f}(\epsilon) = \int_0^\infty \mathrm{e}^{-\epsilon t}f(t)\mathrm{d}t$. The notation $P(x,t|x_0)$ indicates the localized initial condition at x₀ and $J_0(x,t|x_0) = -D\partial_x G_0(x,t|x_0)$ is defined as the free probability current where $\partial_{x_0}\widetilde{G}_0(x,\epsilon|0)$ implies $\frac{\partial}{\partial x_0}\widetilde{G}_0(x,\epsilon|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], $\partial_t P(x,t) = L_x P(x,t)$ where $L_x = -\partial_x A(x) +\partial_x^2 B(x)$ is the FP operator, with the spatially dependent drift and diffusion coefficients given by $A(x) = -(D^2/\kappa)\delta^{^{\prime}}(x)$ and $B(x) = D-(D^2/\kappa)\delta(x)$. As will become apparent later, we are mainly interested in the backward FP (Kolmogorov) equation in terms of the initial variables, t₀ and x₀, $-\partial_{t_0} P(x_0,t_0) = L_{x_0}^\dagger P(x_0,t_0)$, where $L^\dagger_x$ is the formal adjoint of the FP operator i.e. $L^\dagger_x = A(x)\partial_x+B(x) \partial_x^2$. For $P(x_0,t_0 = 0) = \delta(x-x_0)$ and exploiting the time homogeneity of the process, we have [59]

$$\begin{equation} \frac{\partial}{\partial t}P(x,t|x_0) = A(x_0)\frac{\partial}{\partial x_0}P(x,t|x_0)+B(x_0)\frac{\partial^2}{\partial x_0^2} P(x,t|x_0). \end{equation} \tag{ 7 }$$

As was shown in [58], L is in fact self-adjoint, such that the backward FP equation, is given by

$$\begin{equation} \frac{\partial P(x,t|x_0)}{\partial t} = D \frac{\partial^2 P(x,t|x_0)}{\partial x_0^2}-\frac{D^2}{\kappa} \delta^{^{\prime}}(x_0)\partial_{x_0}P(x,t|0), \end{equation} \tag{ 8 }$$

which implies that we can solve equation (8) through the same procedure for which we obtain equation (6).

3.1. Joint probability density $\mathcal{P}(M,t_m,t|0^\pm)$

We study the time-dependent joint distribution, $\mathcal{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 $\mathcal{P}(M,t_m,t|0^+)$ and $\mathcal{P}(M,t_m,t|0^-)$.

To find these joint densities we proceed by using a path decomposition approach [24, 27, 29, 36], 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(\tau): \tau \in [0,t_m]\}$ and the second part is $\{x(\tau): \tau \in [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 $\mathcal{F}(M,t_m|0^\pm)$. 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+\varepsilon,t-t_m|M)$ and taking $\varepsilon\to 0^+$ [24, 29].

We now proceed to find the two quantities, $\mathcal{F}(x,t|x_0)$ and $S(x,t|x_0)$. This calculation was performed in detail in [58], where using equation (8) a backward equation was found for $S(x,t|x_0)$, which was used to find $\mathcal{F}(x,t|x_0)$ in the Laplace domain through $\widetilde{\mathcal{F}}(x,\epsilon|x_0) = 1-\epsilon \widetilde{S}(x,\epsilon|x_0)$:

$$\begin{equation} \widetilde{\mathcal{F}}(x,\epsilon|x_0) = \left\{ \begin{array}{ll} \frac{2\kappa \mathrm{e}^{-|x-x_0|\sqrt{\genfrac{}{}{.3pt}{3}{\epsilon}{D}}}}{\sqrt{D\epsilon}\left[1+\mathrm{e}^{-2|x|\sqrt{\genfrac{}{}{.3pt}{3}{\epsilon}{D}}}\right]+2\kappa}, \ x_0 < 0 < x \ \ \mathrm{or} \ \ x < 0 < x_0,\\ \frac{\left[\sqrt{D \epsilon}+2\kappa\right] \mathrm{e}^{-|x-x_0|\sqrt{\genfrac{}{}{.3pt}{3}{\epsilon}{D}}}+\sqrt{D \epsilon}\mathrm{e}^{-|x+x_0|\sqrt{\genfrac{}{}{.3pt}{3}{\epsilon}{D}}}}{\sqrt{D\epsilon}\left[1+\mathrm{e}^{-2|x|\sqrt{\genfrac{}{}{.3pt}{3}{\epsilon}{D}}}\right]+2\kappa}, \ 0 < x_0 < x \ \ \mathrm{or} \ \ x < x_0 < 0. \end{array} \right . \end{equation} \tag{ 9 }$$

We now proceed to calculate the joint distributions for $x_0 = 0^\pm$.

Case $x_0 = 0^+$.

The joint distribution is given by

$$\begin{equation} \mathcal{P}(M,t_m,t|0^+) = \lim_{\varepsilon\to 0^+} \frac{\mathcal{F}(M,t_m|0^+)S(M+\varepsilon,t-t_m|M)}{N(\varepsilon)}, \end{equation} \tag{ 10 }$$

where $N(\varepsilon)$ ensures the quantity is normalized, alternatively after Laplace transforming with respect to t_m and t ($t_m\to p$ and $t\to \epsilon$) we have,

$$\begin{align} & \widetilde{Q}_m(M,p,\epsilon|0^+) = \lim_{\varepsilon\to 0^+} \frac{1}{N(\varepsilon)} \int_{0}^{\infty} \mathrm{d} t_m \mathrm{e}^{-p t_m} \mathcal{F}(M,t_m|0^+)\\ \nonumber \end{align} \tag{ 11 }$$

$$\begin{align} & \quad\,\,\,\times \int_{0}^{\infty} \mathrm{d} t \mathrm{e}^{-\epsilon t} S(M+\varepsilon,t-t_m|M) = \lim_{\varepsilon\to 0^+} \frac{\widetilde{\mathcal{F}}(M,p+\epsilon|0^+)\widetilde{S}(M+\varepsilon,\epsilon|M)}{N(\varepsilon)}, \end{align}$$

where we have used the notation $Q_m(M,p,t|0^\pm) = \int_0^\infty \mathrm{d}t_m \mathrm{e}^{-pt_m} \mathcal{P}(M,t_m,t|0^\pm)$.

Expanding $\widetilde{S}(M+\varepsilon,\epsilon|M)$ to first order in ε, we have $\widetilde{S}(M+\varepsilon,\epsilon|M)\simeq -\varepsilon\lim_{x\to M}\partial_{x}\widetilde{\mathcal{F}}(x,\epsilon|M)/\epsilon$ and requiring $\mathcal{P}(M,t_m,t|0^\pm)$ to be normalized, namely, $\int_0^\infty \widetilde{Q}_m(M,0,\epsilon|0^\pm)\mathrm{d}M = \epsilon^{-1}$, we find $N(\varepsilon) = \varepsilon$, and then obtain

$$\begin{align} & \widetilde{Q}_m(M,p,\epsilon|0^+)\\ \nonumber & \quad= \frac{2 \mathrm{e}^{M \sqrt{\frac{p+\epsilon }{D}}} \left(\mathrm{e}^{2 M \sqrt{\frac{\epsilon }{D}}} \left(\sqrt{D \epsilon }+2 \kappa \right)-\sqrt{D \epsilon }\right) \left(\sqrt{D (p+\epsilon )}+\kappa \right)}{\sqrt{D \epsilon } \left(\mathrm{e}^{2 M \sqrt{\frac{\epsilon }{D}}} \left(\sqrt{D \epsilon }+2 \kappa \right)+\sqrt{D \epsilon }\right) \left(\mathrm{e}^{2 M \sqrt{\frac{p+\epsilon }{D}}} \left(\sqrt{D (p+\epsilon )}+2 \kappa \right)+\sqrt{D (p+\epsilon )}\right)}. \end{align} \tag{ 12 }$$

One can see from equation (12) that in the no barrier limit, $\kappa\to\infty$, equation (12) correctly reduces to $\widetilde{Q}_m(M,p,\epsilon|0) = \mathrm{e}^{-M \sqrt{(p+\epsilon)/D}}/\sqrt{D \epsilon }$, giving equation (2) after the double Laplace inversion.

Similarly, in the limit of the barrier becoming impermeable, κ → 0, we find,

$$\begin{equation} \widetilde{Q}_m(M,p,\epsilon|0^+) = \frac{\tanh \left(M \sqrt{\frac{\epsilon} {D}}\right) \mathrm{sech}\left(M \sqrt{\frac{p+\epsilon}{D }}\right)}{\sqrt{D \epsilon}} \end{equation} \tag{ 13 }$$

after using the following inverse Laplace transform relations: $\mathcal{L}^{-1}_{\epsilon\to t}\left\{\epsilon^{-1/2} \tanh(a\sqrt{\epsilon}) \right\} = \vartheta_4(0,\mathrm{e}^{-a^2/t})/\sqrt{\pi t}$ and $\mathcal{L}^{-1}_{\epsilon\to t}\left\{\mathrm{sech}(a\sqrt{\epsilon}) \right\} = a\vartheta_1^{^{\prime}}(0,\mathrm{e}^{-a^2/t})/\sqrt{4 \pi t^3}$ [60], equation (13) becomes,

$$\begin{equation} \mathcal{P}(M,t_m,t|0^+) = \frac{M \vartheta_1^{^{\prime}}\left(0,\mathrm{e}^{-\frac{M^2}{D t_m}}\right)\vartheta_4\left(0,\mathrm{e}^{-\frac{M^2}{D (t-t_m)}}\right)}{2 \pi D t_m^{3/2}\sqrt{t-t_m}}, \end{equation} \tag{ 14 }$$

where $\vartheta_1(z,q) = \sum_{n = -\infty}^{\infty}(-1)^{n-1/2}q^{(n+\frac{1}{2})^2}\mathrm{e}^{(2n+1)iz}$ and $\vartheta_4(z,q) = \sum_{n = -\infty}^{\infty}(-1)^n q^{n^2}\mathrm{e}^{2ni z}$ are the Jacobi Theta functions and $\vartheta_n^{^{\prime}}(z,q)$ represents the derivative with respect to z [61].

Interestingly, if we take the small t limit of $\mathcal{P}(M,t_m,t|0^+)$, which corresponds to $\epsilon,p \to \infty$ 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 $\mathcal{P}(M,t_m,t|0^+)$ in terms of two different scaling functions such that,

$$\begin{equation} \mathcal{P}(M,t_m,t|0^+) = \frac{\kappa^3}{D^2}\mathscr{G}^+\left(\frac{\kappa}{D}M,\frac{\kappa^2}{D} t_m\right) \mathscr{H}\left(\frac{\kappa}{D}M,\frac{\kappa^2}{D} (t-t_m)\right), \end{equation} \tag{ 15 }$$

where,

$$\begin{equation} \mathscr{G}^+(y,\tau_1) = \frac{1}{\pi} \int_0^\infty \mathrm{e}^{-\tau_1 z} \sin(y \sqrt{z}) h(y,z) \mathrm{d}z, \end{equation} \tag{ 16 }$$

and

$$\begin{equation} \mathscr{H}(y,\tau_2) = \frac{4}{\pi} \int_0^\infty \frac{\mathrm{e}^{-\tau_2 z}}{\sqrt{z}}h(y,z) \mathrm{d}z, \end{equation} \tag{ 17 }$$

with

$$\begin{equation} h(y,z) = \left[2+z+z\cos(2y\sqrt{z})+2\sqrt{z}\sin(2y\sqrt{z})\right]^{-1}. \end{equation} \tag{ 18 }$$

From equation (15) one can see we no longer have the transformation symmetry of $t_m\to t-t_m$, meaning the symmetry about $t/2$ that one observes in the barrier free case, equation (2), is broken.

Case $x_0 = 0^-$.

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₀. The probability of this occurring is given by (from equation (9))

$$\begin{equation} S(0^\pm,t|0^\mp) = \mathrm{e}^{\frac{\kappa^2t}{D}} \mathrm{erfc}\left(\kappa\sqrt{\frac{t}{D}}\right), \end{equation} \tag{ 19 }$$

where $\mathrm{erfc}(z) = 1-\mathrm{erf}(z)$ where $\mathrm{erf}(z) = (2/\sqrt{\pi})\int_0^z \mathrm{d}u \mathrm{e}^{-u^2}$. This means that the maximum position is the origin and it occurs at time $t_m = 0$. Then from equation (10) we have,

$$\begin{equation} \mathcal{P}(M,t_m,t|0^-) = \lim_{\varepsilon\to 0^+} \frac{1}{N(\varepsilon)} \big[\mathcal{F}(M,t_m|0^-)S(M+\varepsilon,t-t_m|M) +S(0^+,t|0^-)\delta(M)\delta(t_m)\big]. \end{equation} \tag{ 20 }$$

Similar to equation (11) we perform a double Laplace transform and obtain

$$\begin{equation} \widetilde{Q}_m(M,p,\epsilon|0^-) = \lim_{\varepsilon\to 0^+} \frac{1}{N(\varepsilon)}\big[\widetilde{\mathcal{F}}(M,p+\epsilon|0^-)\widetilde{S}(M+\varepsilon,\epsilon|M) +\widetilde{S}(0^+,\epsilon|0^-)\delta(M)\big]. \end{equation} \tag{ 21 }$$

Expanding to first order in ε and ensuring normalization we once again find $N(\varepsilon) = \varepsilon$, leading to

$$\begin{equation} \widetilde{Q}_m(M,p,\epsilon|0^-) = \frac{\kappa}{\sqrt{D(p+\epsilon)}+\kappa}\widetilde{Q}_m(M,p,\epsilon|0^+) +\frac{D\delta(M)}{\kappa \sqrt{D \epsilon }+D \epsilon }. \end{equation} \tag{ 22 }$$

As in equation (12) one can see in the limit $\kappa\to \infty$ equation (22) reduces to the barrier free case, equation (2). In the fully reflecting limit κ → 0 one recovers $\mathcal{P}(M,t_m,t|0^-) = \delta(M)\delta(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 $\widetilde{Q}_m(M,p,\epsilon|0^+)$ for $p,\epsilon \to \infty$, equation (13), we find a better approximation

$$\begin{equation} \widetilde{Q}_m(M,p,\epsilon|0^-)\simeq \frac{\kappa\tanh \left(M \sqrt{\frac{\epsilon} {D}}\right) \mathrm{sech}\left(M \sqrt{\frac{p+\epsilon}{D }}\right)}{D\sqrt{(p+\epsilon)\epsilon}}+\left(\frac{1}{\epsilon}-\frac{\kappa}{\sqrt{D \epsilon^3}} \right)\delta(M), \end{equation} \tag{ 23 }$$

which after using $\mathcal{L}^{-1}_{\epsilon\to t}\left\{\epsilon^{-1/2} \mathrm{sech}(a\sqrt{\epsilon}) \right\} = \theta(\mathrm{e}^{-a^2/t})/\sqrt{\pi t}$ [60], equation (23) becomes

$$\begin{equation} \mathcal{P}(M,t_m,t|0^-)\simeq \frac{\kappa \theta \left(\mathrm{e}^{-\frac{M^2}{D t_m}}\right) \vartheta_4\left(0,\mathrm{e}^{-\frac{M^2}{D (t-t_m)}}\right)}{\pi D \sqrt{t_m(t-t_m)}}+\left(1-\frac{2\kappa\sqrt{t}}{\sqrt{\pi D}}\right)\delta(M)\delta(t_m), \end{equation} \tag{ 24 }$$

where $\theta(q) = \sum_{n = 0}^{\infty}(-1)^nq^{(n+\frac{1}{2})^2}$.

Now we perform the double inverse Laplace transform of equation (22) to obtain (see appendix B)

$$\begin{align} \mathcal{P}(M,t_m,t|0^-)& = \frac{\kappa^3}{D^2}\mathscr{G}^-\left(\frac{\kappa}{D}M,\frac{\kappa^2}{D} t_m\right) \mathscr{H}\left(\frac{\kappa}{D}M,\frac{\kappa^2}{D} (t-t_m)\right)\\\nonumber &\quad +\frac{\kappa^3}{D^2}\mathscr{I}\left(\frac{\kappa}{D}M,\frac{\kappa^2}{D} t_m,\frac{\kappa^2}{D} (t-t_m)\right) \end{align} \tag{ 25 }$$

where $\mathscr{I}(y,\tau_1,\tau_2) = \mathrm{e}^{\tau_2}\mathrm{erfc}(\sqrt{\tau_2})\delta(y)\delta(\tau_1)$, $\mathscr{H}(y,\tau_2)$ is defined in equation (17) and

$$\begin{equation} \mathscr{G}^-(y,\tau_1) = \int_0^\infty \frac{\mathrm{e}^{-\tau_1 z}}{\pi} \left(\sin(y \sqrt{z})+\sqrt{z}\cos(y \sqrt{z})\right) h(y,z) \mathrm{d}z, \end{equation} \tag{ 26 }$$

where $h(y,z)$ is defined in equation (18). Again we see the $t_m\to t-t_m$ symmetry breaking as in the $x_0 = 0^+$ case.

3.2. Marginal density $\mathcal{P}(M,t|0^\pm)$

To find the marginal density $\mathcal{P}(M,t|0^\pm)$ 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

$$\begin{equation} \mathcal{P}(M,t|0^\pm) = \frac{\kappa}{D}\mathcal{I}^\pm\left(\frac{\kappa}{D}M,\frac{\kappa^2}{D}t\right) \end{equation} \tag{ 27 }$$

where

$$\begin{equation} \mathcal{I}^+(y,\tau) = \frac{1}{\pi}\int_0^\infty \frac{\mathrm{e}^{-\tau z}}{\sqrt{z}}j^+(y,z)h^2(y,z)\mathrm{d}z, \end{equation} \tag{ 28 }$$

and

$$\begin{equation} \mathcal{I}^-(y,\tau) = \frac{1}{\pi}\int_0^\infty \frac{\mathrm{e}^{-\tau z}}{\sqrt{z}}j^-(y,z)h^2(y,z)\mathrm{d}z + \mathrm{e}^\tau\mathrm{erfc}(\sqrt{\tau})\delta(y), \end{equation} \tag{ 29 }$$

with $j^\pm(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 ($\kappa \to \infty$) one recovers equation (3).

From equations (27)–(29), we find that for $t\gt\gt D/\kappa^2$ with $M\sim D/\kappa$, the integrands in (28) and (29) are dominated by small z, so we expand $j^\pm(y,z)h^2(y,z)$ to first order in z and then compute the integral to obtain,

$$\begin{equation} \mathcal{P}(M,t|0^+)\simeq \frac{4 \kappa D (\kappa t-2M)-\kappa^2 M^2-2D^2}{4 \sqrt{\pi}\kappa^2 (Dt)^{3/2}} \end{equation} \tag{ 30 }$$

and

$$\begin{equation} \mathcal{P}(M,t|0^-)\simeq \frac{4 \kappa D (\kappa t-\frac{5}{2}M)-\kappa^2 M^2-6D^2}{4 \sqrt{\pi}\kappa^2 (Dt)^{3/2}}+\frac{\delta(M)}{\kappa\sqrt{\pi D t}}. \end{equation} \tag{ 31 }$$

Equations (30) and (31) should be compared with the barrier free case, $\mathcal{P}(M,t|0)\simeq (Dt-M^2/4)/\sqrt{\pi D^3 t^3}$.

We plot equation (27) in figure 2 for the two initial conditions, $x_0 = 0^\pm$ whilst varying the dimensionless parameter, $\kappa^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).

Zoom In Zoom Out Reset image size

3.3. Marginal density $\mathcal{P}(t_m,t|0^\pm)$

Here we consider the marginal density $\mathcal{P}(t_m,t|0^\pm) = \int_0^\infty \mathcal{P}(M,t_m,t|0^\pm)\mathrm{d}M$, where from equations (15) and (25), we obtain the scaling relation,

$$\begin{equation} \mathcal{P}(t_m,t|0^\pm) = \frac{\kappa^2}{D}\mathscr{F}^\pm_m\left(\frac{\kappa^2}{D}t_m,\frac{\kappa^2}{D}(t-t_m)\right), \end{equation} \tag{ 32 }$$

where

$$\begin{equation} \mathscr{F}^\pm_m(\tau_1,\tau_2) = \int_0^\infty \mathscr{G}^\pm(y,\tau_1)\mathscr{H}(y,\tau_2)\mathrm{d}y. \end{equation} \tag{ 33 }$$

The integral in equation (33) is hard to compute, but we can study the asymptotic forms of $\mathscr{F}^\pm_m(\tau_1,\tau_2)$. Firstly, we study the short time asymptotics, $t\lt\lt D/\kappa^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

$$\begin{align} \mathscr{F}_m^+(\tau_1,\tau_2)&\simeq \sum_{n = -\infty}^\infty \sum_{m = -\infty}^\infty \frac{(-1)^{m+n} (2 n+1) \sqrt{\tau_2} }{\pi \sqrt{\tau_1} \left(4 m^2 \tau_1+(2 n+1)^2 \tau_2\right)}\\ \nonumber &= \sum _{n = -\infty }^{\infty } \frac{(-1)^n \mathrm{cosech} \left((n+\frac{1}{2})\pi \sqrt{\frac{\tau_2}{\tau_1}}\right)}{2 \tau_1}, \end{align} \tag{ 34 }$$

where we used $\sum_{m = -\infty}^{\infty} (-1)^m (m^2+z)^{-1} = \pi \mathrm{cosech}(\pi \sqrt{z})z^{-1/2}$. And,

$$\begin{equation} \mathscr{F}_m^-(\tau_1,\tau_2)\simeq \sum_{n = -\infty}^\infty \sum_{m = -\infty}^\infty (-1)^{m+n}(4 \pi m^2 \tau_1+(2 n+1)^2\pi \tau_2)^{-1/2}. \end{equation} \tag{ 35 }$$

Now we study the long time asymptotics of $\mathcal{P}(t_m,t|0^\pm)$, which corresponds to $t\gt\gt D/\kappa^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\to \infty$: keeping t_m finite in the first regime and keeping $t-t_m$ finite in the second. This corresponds to ε → 0 with $p \gt\gt \epsilon$ and $p \sim \epsilon$ in the Laplace domain, respectively. Expanding around $\epsilon,p\to0$ to leading order in equations (12) and (22) and integrating over M one can see that we obtain, $Q_m(p,\epsilon|0^\pm) = 1/\sqrt{\epsilon(p+\epsilon)}$, which after double inverse Laplace transforming recovers the Arcsine distribution, (1). For the case $p\gt\gt\epsilon$, let us expand $\widetilde{Q}_m(p,\epsilon|0^+)$ around ε → 0 keeping p finite, then we find

$$\begin{equation} \widetilde{Q}_m(p,\epsilon|0^+)\simeq \frac{2 (\sqrt{D p}+\kappa)}{\sqrt{D \epsilon}}\int_0^\infty \frac{\mathrm{e}^{-M \sqrt{\frac{p}{D}}}}{\sqrt{D p}\left(1+\mathrm{e}^{-M \sqrt{\frac{p}{D}}}\right)+2 \kappa}\mathrm{d}M, \end{equation} \tag{ 36 }$$

which gives after performing the integral and inverse Laplace transforming with respect to ε,

$$\begin{equation} Q_m(p,t|0^+)\simeq \frac{2 \left(\sqrt{D p}+\kappa \right) \mathrm{arctan}\left(\frac{(D p)^{1/4}}{(\sqrt{D p}+2 \kappa )^{1/2}}\right)}{\sqrt{\pi t}(D p)^{1/4} \left[p \left(\sqrt{D p}+2 \kappa \right)\right]^{1/2}}. \end{equation} \tag{ 37 }$$

Since equation (37) is very difficult invert, we can investigate the asymptotic dependence for $t_m\to 0$. This corresponds to $p\to \infty$, thus by expanding equation (37) to leading order for $p\to \infty$ and using $\mathrm{arctan}(1) = \pi/4$, we find $Q_m(p,t|0^+)\sim \frac{\sqrt{\pi} }{2 \sqrt{p t }}$, and performing the Laplace inversion with respect to p we get,

$$\begin{equation} \mathcal{P}(t_m,t|0^+)\sim \frac{1}{2 \sqrt{t_m t}} \ \mathrm{for} \ t_m\to 0 ,\ t\to \infty. \end{equation} \tag{ 38 }$$

For the case $x_0 = 0^-$, using equation (22), we find,

$$\begin{equation} \mathcal{P}(t_m,t|0^-)\sim \frac{\sqrt{\pi}\kappa}{2 \sqrt{D t}} + \frac{\sqrt{D}}{\kappa \sqrt{\pi t}} \delta(t_m) \ \mathrm{for} \ t_m\to 0 ,\ t\to \infty. \end{equation} \tag{ 39 }$$

It is clear that equations (38) and (39) deviate from the Arcsine law, $\sim\!\!1/\pi(t_m t)^{-1/2}$, showing the permeable barrier has an influence for $t\to \infty$ when t_m is finite. And for the $x_0 = 0^-$ case the singularity at $t_m\to0$ 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^\pm$, $\mathcal{P}(t_m,t|0^\pm)$ 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$.

Zoom In Zoom Out Reset image size

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) = \int_0^t \Theta(x(t^{^{\prime}}))\mathrm{d}t^{^{\prime}}$, where $\Theta(z)$ is the Heaviside step function. From Feynman–Kac theory [62, 63], we know that $Q_r(p,t|x_0) = \int_0^\infty \mathrm{e}^{-p t_r} \mathcal{P}(t_r,t|x_0)\mathrm{d}t_r = \langle \mathrm{e}^{-p t_r(t)}| x_0 = x(0)\rangle$ satisfies the following backward Feynman–Kac equation,

$$\begin{equation} \frac{\partial Q_r(p,t|x_0)}{\partial t} = A(x_0)\frac{\partial Q_r(p,t|x_0)}{\partial x_0}+B(x_0)\frac{\partial^2 Q_r(p,t|x_0)}{\partial x_0^2} -p \Theta(x_0)Q_r(p,t|x_0), \end{equation} \tag{ 40 }$$

where we use the backward FP operator in equation (7). Using the self-adjoint nature of this backward FP operator, we may write equation (40) as

$$\begin{equation} \frac{\partial Q_r(p,t|x_0)}{\partial t} = D\frac{\partial^2 Q_r(p,t|x_0)}{\partial x_0^2}-p\Theta(x_0)Q_r(p,t|x_0)-\frac{D^2}{\kappa} \delta^{^{\prime}}(x_0)\partial_{x_0}Q_r(p,t|0), \end{equation} \tag{ 41 }$$

with the initial condition $Q_r(p,0|x_0) = 1$. In the Laplace domain, $t\to\epsilon$, the solution to (41) is given by (see appendix A),

$$\begin{equation} \widetilde{Q}_r(p,\epsilon|x_0) = \widetilde{\mathcal{Q}}(p,\epsilon|x_0) +\partial_{x_0} \widetilde{\mathcal{Q}}(p,\epsilon|0) \frac{\partial_x \widetilde{\mathcal{G}}(0,p,\epsilon|x_0) }{\frac{\kappa}{D^2}-\partial^2_{x,x_0}\widetilde{\mathcal{G}}(0,p,\epsilon|0)}, \end{equation} \tag{ 42 }$$

where $\mathcal{G}(x,p,t|x_0)$ is the Green's function of the barrier free Feynman–Kac equation, i.e.

$$\begin{equation} \frac{\partial \mathcal{G}(x,p,t|x_0)}{\partial t} = D\frac{\partial^2 \mathcal{G}(x,p,t|x_0)}{\partial x_0^2}-p\Theta(x_0)\mathcal{G}(x,p,t|x_0), \end{equation} \tag{ 43 }$$

with $\mathcal{G}(x,p,0|x_0) = \delta(x_0-x)$ and $\mathcal{Q}(p,t|x_0) = \int_{-\infty}^{\infty}\mathcal{G}(x,p,t|x_0)\mathrm{d}x$.

The solution of equation (43) can be found by solving in the Laplace domain,

$$\begin{equation} \epsilon \widetilde{\mathcal{G}}(x,p,\epsilon|x_0)-\delta(x_0-x) = \left\{ \begin{array}{ll} D \frac{\partial^2}{\partial x_0^2} \widetilde{\mathcal{G}}(x,p,\epsilon|x_0)- p \widetilde{\mathcal{G}}(x,p,\epsilon|x_0) \ \ x_0>0,\\[4pt] D \frac{\partial^2}{\partial x_0^2} \widetilde{\mathcal{G}}(x,p,\epsilon|x_0) \ \ x_0<0. \end{array} \right . \end{equation} \tag{ 44 }$$

Two boundary conditions are required with equation (44) for $x_0\to\pm\infty$. As $\mathcal{P}(t_r,t|x_0\to\infty) = \delta(t_r-t)$ and $\mathcal{P}(t_r,t|x_0\to-\infty) = \delta(t_r)$, this corresponds to $\widetilde{\mathcal{G}}(x,p,\epsilon|x_0\to \infty) = (\epsilon+p)^{-1}$ and $\widetilde{\mathcal{G}}(x,p,\epsilon|x_0\to -\infty) = \epsilon^{-1}$. In addition, we require continuity at the origin for $\widetilde{\mathcal{G}}(x,p,\epsilon|x_0)$ and its derivative [64], i.e. $\widetilde{\mathcal{G}}(x,p,\epsilon|0^-) = \widetilde{\mathcal{G}}(x,p,\epsilon|0^+)$ and $\lim_{x_0\to0^-}\partial_{x_0}\widetilde{\mathcal{G}}(x,p,\epsilon|x_0) = \lim_{x_0\to0^+}\partial_{x_0}\widetilde{\mathcal{G}}(x,p,\epsilon|x_0)$. The presence of the Dirac-δ function on the left-hand side (LHS) means we have two more conditions, such that we have continuity at $x_0 = x$ therefore $\widetilde{\mathcal{G}}(x,p,\epsilon|x^-) = \widetilde{\mathcal{G}}(x,p,\epsilon|x^+)$ and from integrating over the Dirac-δ, we have $\lim_{x_0\to x^-}\partial_{x_0}\widetilde{\mathcal{G}}(x,p,\epsilon|x_0)-\lim_{x_0\to x^+}\partial_{x_0}\widetilde{\mathcal{G}}(x,p,\epsilon|x_0) = 1/D$.

Solving equation (44) with the aforementioned conditions to obtain $\widetilde{\mathcal{G}}(x,p,\epsilon|x_0)$ and inserting into equation (42) gives $\widetilde{Q}_r(p,\epsilon|x_0)$. Setting $x_0 = 0^+$ and $x_0 = 0^-$ leads to,

$$\begin{equation} \widetilde{Q}_r(p,\epsilon|0^+) = \frac{\sqrt{D} \epsilon +\kappa \left(\sqrt{p+\epsilon }+\sqrt{\epsilon }\right)}{\sqrt{D} \epsilon (p+\epsilon )+\kappa \sqrt{\epsilon } \left(\sqrt{\epsilon (p+\epsilon )}+p+\epsilon \right)}, \end{equation} \tag{ 45 }$$

and

$$\begin{equation} \widetilde{Q}_r(p,\epsilon|0^-) = \frac{\sqrt{D} (p+\epsilon )+\kappa \left(\sqrt{p+\epsilon }+\sqrt{\epsilon }\right)}{\sqrt{D} \epsilon (p+\epsilon )+\kappa \sqrt{\epsilon } \left(\sqrt{\epsilon (p+\epsilon )}+p+\epsilon \right)}. \end{equation} \tag{ 46 }$$

Instantly one can see that in the no barrier limit, $\kappa \to \infty$, we recover the Arcsine law, $\widetilde{Q}_r(p,\epsilon|0) = 1/\sqrt{\epsilon(p+\epsilon)}$. For finite κ the double inverse Laplace transform ($p\to t_r$ and $\epsilon \to t$) of equations (45) and (46) can be written in terms of the scaling relation (see appendix D),

$$\begin{equation} \mathcal{P}(t_r,t|0^\pm) = \frac{\kappa^2}{D}\mathscr{F}_r^\pm\left(\frac{\kappa^2}{D}t_r,\frac{\kappa^2}{D}(t-t_r)\right) \end{equation} \tag{ 47 }$$

where,

$$\begin{equation} \mathscr{F}_r^+(\tau_1,\tau_2) = \mathscr{f}^+(\tau_1,\tau_2)+\unicode{x210A}^+(\tau_1,\tau_2)+\mathrm{e}^{\tau_1}\mathrm{erfc}(\sqrt{\tau_1})\delta(\tau_2), \end{equation} \tag{ 48 }$$

and

$$\begin{equation} \mathscr{F}_r^-(\tau_1,\tau_2) = \mathscr{f}^-(\tau_1,\tau_2)+\unicode{x210A}^-(\tau_1,\tau_2)+\mathrm{e}^{\tau_2}\mathrm{erfc}(\sqrt{\tau_2})\delta(\tau_1). \end{equation} \tag{ 49 }$$

The Dirac-δ functions appear due to the particle spending the whole time in the positive region or none of the time in that region, and so are multiplied the probability of never crossing the barrier for the whole time t, i.e. equation (19). The scaling functions in equations (48) and (49) are given by (see appendix D),

$$\begin{equation} \mathscr{f}^+(\tau_1,\tau_2) = \frac{1}{\sqrt{\pi \tau_1}}e^{\tau_2}\mathrm{erfc}(\sqrt{\tau_2}), \end{equation} \tag{ 50 }$$

$$\begin{equation} \mathscr{f}^-(\tau_1,\tau_2) = \frac{2}{\pi}\sqrt{\frac{\tau_2}{\tau_1}}-2\tau_2\mathscr{f}^+(\tau_1,\tau_2), \end{equation} \tag{ 51 }$$

and

$$\begin{equation} \unicode{x210A}^\pm(\tau_1,\tau_2) = \sum_{n = 0}^{\infty}\frac{(-1)^n \tau_1^{n/2}}{\Gamma(\frac{n}{2}+1)} \mathcal{C}_n^\pm(\tau_2), \end{equation} \tag{ 52 }$$

where $\mathcal{C}_n^+(\tau_2)$ and $\mathcal{C}_n^-(\tau_2)$ are detailed in equations (D.22) and (D.29) and $\Gamma(z) = \int_0^\infty t^{z-1}\mathrm{e}^{-t}\mathrm{d}t$ is the Gamma function.

By having the analytical expressions for $\mathcal{P}(t_r,t|0^\pm)$, we can study the asymptotics of this distribution. For the short time asymptotics, $t\lt\lt D/\kappa^2$, we take $\tau_1,\tau_2 \to 0$ in $\mathscr{F}_r^\pm(\tau_1,\tau_2)$ and use $\mathscr{f}^+(\tau_1,\tau_2)\simeq 1/\sqrt{\pi \tau_1}$ and $\unicode{x210A}^+(\tau_1,\tau_2)\simeq 2\sqrt{\tau_1/\tau_2}/\pi$ (see appendix E) to give

$$\begin{equation} \mathcal{P}(t_r,t|0^+)\simeq \frac{\kappa}{\sqrt{\pi D t_r}}+\frac{2\kappa^2 \sqrt{t_r}}{\pi D \sqrt{t-t_r}}+\left(1-2\kappa \sqrt{\frac{t_r}{\pi D}}\right)\delta(t-t_r). \end{equation} \tag{ 53 }$$

Then from equation (51) we have $\mathscr{f}^-(\tau_1,\tau_2)\simeq 2 \sqrt{\tau_2/\tau_1}/\pi$ and $\unicode{x210A}^-(\tau_1,\tau_2)\simeq 1/\sqrt{\pi \tau_2}$ (see appendix E) for $\tau_1,\tau_2\to0$, which leads to the following symmetrical dependence for $t\lt\lt D/\kappa^2$,

$$\begin{equation} \mathcal{P}(t_r,t|0^+) = \mathcal{P}(t-t_r,t|0^-). \end{equation} \tag{ 54 }$$

The long time asymptotics, $t\gt\gt D/\kappa^2$, corresponds to the limit ε → 0 in the Laplace domain, by expanding $Q_r(p,\epsilon|0^\pm)$ in equations (45) and (46) around small ε to leading order whilst keeping $p+\epsilon$ finite, one obtains $Q_r(p,\epsilon|0^\pm)\simeq 1/\sqrt{\epsilon(p+\epsilon)}$, leading to the Arcsine distribution,

$$\begin{equation} \mathcal{P}(t_r,t|0^\pm)\simeq\frac{1}{\pi\sqrt{t_r(t-t_r)}}, \end{equation} \tag{ 55 }$$

showing how the influence of the permeable barrier on $\mathcal{P}(t_r,t|0^\pm)$ wanes over time.

We plot $\mathcal{P}(t_r,t|0^\pm)$ in figure 4 to show the excellent match with stochastic simulations. By varying the strength of the permeability of the barrier the resulting curves move further from the Arcsine distribution, where for smaller κ, the peaks at $t_r = 0$ and $t_r = t$ become sharper. This illustrates the notion that once the particle crosses the barrier it is unlikely to do so again for small permeabilities.

Zoom In Zoom Out Reset image size

Finally, we consider the third Arcsine law, the probability density, $\mathcal{P}(t_\ell,t|0^\pm)$, 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 $\mathcal{P}(t_\ell,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(\tau): \tau \in [0,t_\ell]\}$ and $\{x(\tau): \tau \in [t_\ell,t] \}$, where the first part is the trajectory that crosses the origin at $t_\ell$, and the second part is that the trajectory does not reach the origin for the remaining time $t-t_\ell$. 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 $\mathcal{P}(t_\ell,t|0^+)$ as,

$$\begin{equation} \mathcal{P}(t_\ell,t|0^+) = \lim_{\varepsilon\to 0^+} \frac{1}{N(\varepsilon)}\Big[P(0^-,t_\ell|0^+)S(-\varepsilon,t-t_\ell|0^-) +P(0^+,t_\ell|0^+)S(\varepsilon,t-t_\ell|0^+)\Big], \end{equation} \tag{ 56 }$$

where the sum in equation (56) comes from the last passage being an up crossing or a down crossing. We take the limit $\varepsilon\to 0^+$ to find the double Laplace transform of $\mathcal{P}(t_\ell,t|0^+)$, i.e. $\widetilde{Q}_\ell(p,\epsilon|0^+) = \int_0^\infty \mathrm{d}t \mathrm{e}^{-\epsilon t}\int_0^\infty \mathrm{d}t_\ell \ \mathrm{e}^{-pt_\ell} \mathcal{P}(t_\ell,t|0^+)$, such that

$$\begin{equation} {Q_\ell }(p,|{0^ + }) = \mathop {\lim }\limits_{\varepsilon \to {0^ + }} \frac{1}{{N(\varepsilon )}}[\tilde P({0^ - },p + |{0^ + })\tilde S( - \varepsilon ,|{0^ - }) + \tilde P({0^ + },p + |{0^ + })\tilde S(\varepsilon ,|{0^ + })], \end{equation} \tag{ 57 }$$

where $\widetilde{P}(x,\epsilon|x_0)$ is given by equation (6) and $\widetilde{S}(\pm\varepsilon,\epsilon|0^\pm)$ is found to first order in ε as previously. To find the normalization factor $N(\epsilon)$, we use the normalization condition,

$$\begin{equation} \int_0^t \mathcal{P}(t_\ell,t|0^+)\mathrm{d}t_\ell = 1-S(0^-,t|0^+), \end{equation} \tag{ 58 }$$

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(\varepsilon) = \varepsilon/D$, giving

$$\begin{equation} \widetilde{Q}_\ell(p,\epsilon|0^+) = \frac{\kappa}{\sqrt{\epsilon(p+\epsilon)}(\sqrt{D \epsilon}+\kappa)}, \end{equation} \tag{ 59 }$$

using standard inverse Laplace transform relations [60], we obtain

$$\begin{equation} \mathcal{P}(t_\ell,t|0^+) = \frac{\kappa^2}{D}\mathscr{F}_\ell\left(\frac{\kappa^2}{D}t_\ell,\frac{\kappa^2}{D}(t-t_\ell)\right), \end{equation} \tag{ 60 }$$

where

$$\begin{equation} \mathscr{F}_\ell(\tau_1,\tau_2) = \mathscr{f}^+(\tau_1,\tau_2) \end{equation} \tag{ 61 }$$

for $\mathscr{f}^+(\tau_1,\tau_2)$ defined in equation (50). The striking feature here is that the scaling function that fully describes $\mathcal{P}(t_\ell,t|0^+)$ is also part of the scaling function describing the residence time density, $\mathcal{P}(t_r,t|0^+)$. Comparing the short time limit, $t\lt\lt D/\kappa^2$,

$$\begin{equation} \mathcal{P}(t_\ell,t|0^+)\simeq \frac{\kappa}{\sqrt{\pi D t_\ell}}, \end{equation} \tag{ 62 }$$

and the long time limit, $t\gt\gt D/\kappa^2$ for $t_\ell\lt\lt t$ ,

$$\begin{equation} \mathcal{P}(t_\ell,t|0^+)\simeq \frac{1}{\pi\sqrt{t_\ell(t-t_\ell)}}, \end{equation} \tag{ 63 }$$

to that of $\mathcal{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_\ell \lt\lt 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_\ell$. This breaks down as $t_\ell$ gets larger, causing the different dependencies of the distributions. In this case equation (63) is not valid for $t-t_\ell \to 0$ as $\mathscr{f}^+(\tau_1,0) = 1/\sqrt{\pi \tau_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_\ell = 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_\ell = 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.

Zoom In Zoom Out Reset image size

In summary, we have investigated the EVS 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, $\mathcal{P}(M,t_m,t|0^\pm)$. 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\to 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\lt\lt D/\kappa^2$, $\mathcal{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 managing to go through the barrier for short times. For $\mathcal{P}(M,t_m,t|0^-)$ we find a stronger approximation than just the reflecting barrier distribution, which is $\delta(M)\delta(t_m)$, which accounts for the very few particles that do cross the barrier.

We have also investigated the respective marginal distributions, $\mathcal{P}(M,t|0^\pm)$ and $\mathcal{P}(t_m,t|0^\pm)$. The presence of the barrier has a large impact on the monotonicity of the distribution, $\mathcal{P}(M,t|0^\pm)$, where for certain permeabilities a maximum appears away from the origin. At long times, $t\gt\gt D/\kappa^2$, the distribution is still dependent on κ. For $\mathcal{P}(t_m,t|0^\pm)$ 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\to 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_\ell(t)$. Using Feynman–Kac theory we have calculated $\mathcal{P}(t_r,t|0^\pm)$ 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^\pm$ are equivalent under the transformation $t_r\to t-t_r$, and in the long time limit we recover the Arcsine distribution.

Finally, we have studied $\mathcal{P}(t_\ell,t|0^\pm)$, which is equivalent for either initial condition $x_0 = 0^+$ or $x_0 = 0^-$, since $t_\ell(t)$ is a crossing event. Taking into account up and down crossings we have found $\mathcal{P}(t_\ell,t|0^\pm)$ in terms of known functions. Interestingly the scaling function that describes this distribution happens to be part of the scaling function which describes $\mathcal{P}(t_r,t|0^+)$, where the peak at the origin of both distributions have the same dependence, because for $t_\ell\lt\lt t$ the first crossing very likely corresponds to the last crossing, which implies $t_\ell = t_r$. As $t_\ell$ becomes larger this is no longer the case, and we do not observe a divergence as $t_\ell \to 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^\pm$ 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.

T K and L G acknowledge funding from, respectively, an Engineering and Physical Sciences Research Council (EPSRC) DTP student grant and the Biotechnology and Biological Sciences Research Council (BBSRC) Grant No. BB/T012196/1 and NERC Grant No. NE/W00545X/1. This work was carried out using the computational facilities of the Advanced Computing Research Centre, University of Bristol—www.bristol.ac.uk/acrc/.

All data that support the findings of this study are included within the article (and any supplementary files).

Here we show how the solution of equation (41), $Q_r(p,t|x_0)$, can be represented in terms the Green's function, $\mathcal{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,

$$\begin{align} Q_r(p,t|x_0)& = \int_{-\infty}^{\infty} \mathrm{d}y \mathcal{G}(y,p,t|x_0) Q_r(p,0|y)\\ \nonumber &\quad - \frac{D^2}{\kappa} \int_0^t \mathrm{d}t^{^{\prime}} \int_{-\infty}^{\infty} \mathrm{d}y \mathcal{G}(y,p,t-t^{^{\prime}}|x_0) \delta^{^{\prime}}(x_0) \partial_{x_0}Q_r(p,t^{^{\prime}}|0). \end{align} \tag{ A.1 }$$

Using $Q_r(p,0|y) = 1$ with $\mathcal{Q}(p,t|x_0) = \int_{-\infty}^\infty \mathcal{G}(y,p,t|x_0)$ and Laplace transforming, $t\to \epsilon$, we have

$$\begin{equation} \widetilde{Q}_r(p,\epsilon|x_0) = \widetilde{\mathcal{Q}}(p,\epsilon|x_0)+\frac{D^2}{\kappa} \partial_x \widetilde{\mathcal{G}}(0,p,\epsilon|x_0) \partial_{x_0} \widetilde{Q}_r(0,\epsilon|0). \end{equation} \tag{ A.2 }$$

Then by taking the derivative of both sides of equation (A.2) with respect to x₀ and setting $x_0 = 0$, we find

$$\begin{equation} \partial_{x_0} \widetilde{Q}_r(0,\epsilon|0) = \frac{\partial_{x_0}\widetilde{\mathcal{Q}}(p,\epsilon|0)}{1-\frac{D^2}{\kappa}\partial^2_{x,x_0}\widetilde{\mathcal{G}}(0,p,\epsilon|0)}, \end{equation} \tag{ A.3 }$$

then after inserting equation (A.3) into (A.2) we obtain (42).

To perform the Laplace inversion of $\widetilde{Q}_m(M,p,\epsilon|0^+)$, we write equation (12) as

$$\begin{equation} \widetilde{Q}_m(M,p,\epsilon|0^+) = \frac{1}{\kappa}\widetilde{\mathscr{G}}^+\left(\frac{\kappa}{D}M,\frac{D}{\kappa^2}(p+\epsilon)\right)\widetilde{\mathscr{H}}\left(\frac{\kappa}{D}M,\frac{D}{\kappa^2}\epsilon\right), \end{equation} \tag{ B.1 }$$

then $\mathcal{P}(M,t_m,t|0^+)$ is given by the two Bromwich integrals,

$$\begin{align} \mathcal{P}(M,t_m,t|0^+)& = \frac{\kappa^3}{D^2} \frac{1}{(2\pi i)^2} \int_{\gamma_1 -i \infty}^{\gamma_1+i\infty} \mathrm{d}s_1 \mathrm{e}^{\frac{\kappa^2}{D}t_m s_1 } \widetilde{\mathscr{G}}^+\left(\frac{\kappa}{D}M,s_1\right)\\\nonumber & \quad \times \int_{\gamma_2 -i \infty}^{\gamma_2+i\infty} \mathrm{d}s_2 \mathrm{e}^{\frac{\kappa^2}{D}(t-t_m) s_2} \widetilde{\mathscr{H}}\left(\frac{\kappa}{D}M,s_2\right), \end{align} \tag{ B.2 }$$

where γ₁ and γ₂ are greater than the real part of all singularities of $\widetilde{\mathscr{G}}^+\left(\frac{\kappa}{D}M,s_1\right)$ and $\widetilde{\mathscr{H}}\left(\frac{\kappa}{D}M,s_2\right)$, respectively and $\widetilde{\mathscr{G}}^+(y,s_1)$ and $\widetilde{\mathscr{H}}(y,s_2)$ are given by,

$$\begin{equation} \widetilde{\mathscr{G}}^+(y,s_1) = \frac{\mathrm{e}^{y\sqrt{s_1}}(\sqrt{s_1}+1)}{\sqrt{s_1}+\mathrm{e}^{2y\sqrt{s_1}}(\sqrt{s_1}+2)} \end{equation} \tag{ B.3 }$$

and

$$\begin{equation} \widetilde{\mathscr{H}}(y,s_2) = \frac{2\left(-\sqrt{s_2}+\mathrm{e}^{2y\sqrt{s_2}}(\sqrt{s_2}+2)\right)}{\sqrt{s_2}\left(\sqrt{s_2}+\mathrm{e}^{2 y\sqrt{s_2}}(\sqrt{s_2}+2)\right)}. \end{equation} \tag{ B.4 }$$

To find $\mathcal{P}(M,t_m,t|0^+)$ we require the following Laplace inversions $\mathcal{L}^{-1}_{s_1\to\tau_1}\{\widetilde{\mathscr{G}}^+(y,s_1) \}$ and $\mathcal{L}^{-1}_{s_2\to\tau_2}\{\widetilde{\mathscr{H}}(y,s_2) \}$, then $\mathcal{P}(M,t_m,t|0^-)$ only requires the Laplace inversion, $\mathcal{L}^{-1}_{s_1\to\tau_1}\{\widetilde{\mathscr{G}}^-(y,s_1) \}$, where $\widetilde{\mathscr{G}}^-(y,s_1) = (\sqrt{s_1}+1)^{-1}\widetilde{\mathscr{G}}^+(y,s_1)$.

B.1. Laplace inversion of $\widetilde{\mathscr{G}}^\pm(y,s_1)$

From the definition of $\widetilde{\mathscr{G}}^+(y,s_1)$ in equation (B.3) we see that $\widetilde{\mathscr{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 $\widetilde{\mathscr{G}}^+(y,s_1)$ is analytic inside the contour, C, in figure B1, meaning that $\oint_{C} \mathrm{e}^{s_1 \tau_1} \widetilde{\mathscr{G}}^+(y,s_1) \mathrm{d}s_1 = 0$. Then by taking $R\to \infty$ and α → 0 for the contour C the Laplace inversion is given by,

$$\begin{equation} \mathcal{L}^{-1}_{s_1\to \tau_1}\left\{\widetilde{\mathscr{G}}^+(y,s_1)\right\} = \frac{1}{2\pi i} \int_{C_1} \mathrm{e}^{s_1 \tau_1} \widetilde{\mathscr{G}}^+(y,s_1) \mathrm{d}s = \frac{-1}{2\pi i}\left(\int_{C_3}+\int_{C_5}\right) \mathrm{e}^{s_1 \tau_1} \widetilde{\mathscr{G}}^+(y,s_1) \mathrm{d}s_1 \end{equation} \tag{ B.5 }$$

because in the limit $R\to\infty$ the contributions from C₂ and C₆ vanish, and for α → 0 the contribution from C₄ is zero. Therefore, all we require is finding the integrals $\int_{C_3}$ and $\int_{C_5}$.

Zoom In Zoom Out Reset image size

For $\int_{C_3}$ we let $s_1 = z\mathrm{e}^{i \pi}$, giving

$$\begin{equation} \int_{C_3} \mathrm{e}^{s_1 \tau_1} \widetilde{\mathscr{G}}^+(y,s_1) \mathrm{d}s_1 = -\int_{R}^{\alpha} \mathrm{e}^{- \tau_1 z } \frac{\mathrm{e}^{i y\sqrt{z}}(i\sqrt{z}+1)}{i\sqrt{z}+\mathrm{e}^{2iy\sqrt{z}}(i\sqrt{z}+2)}\mathrm{d}z, \end{equation} \tag{ B.6 }$$

and for $\int_{C_5}$ we let $s_1 = z\mathrm{e}^{-i \pi}$, leading to

$$\begin{equation} \int_{C_5} \mathrm{e}^{s_1 \tau_1} \widetilde{\mathscr{G}}^+(y,s_1) \mathrm{d}s_1 = -\int_{\alpha}^{R} \mathrm{e}^{- \tau_1 z } \frac{\mathrm{e}^{-i y\sqrt{z}}(-i\sqrt{z}+1)}{-i\sqrt{z}+\mathrm{e}^{-2iy\sqrt{z}}(-i\sqrt{z}+2)}\mathrm{d}z. \end{equation} \tag{ B.7 }$$

Taking $R\to\infty$ and α → 0, substituting back into equation (B.5) and converting the complex exponentials to trigonometric functions, we obtain $\mathscr{G}^+(y,\tau_1)$ as defined in equation (16).

Similarly, since $\widetilde{\mathscr{G}}^-(y,s_1) = (\sqrt{s_1}+1)^{-1}\widetilde{\mathscr{G}}^+(y,s_1)$, one can see that by altering the above calculations it leads to the definition of $\mathscr{G}^-(y,\tau_1)$ in equation (26).

B.2. Laplace Inversion of $\widetilde{\mathscr{H}}(y,s_2)$

From equation (B.4) one can see that $\widetilde{\mathscr{H}}(y,s_2)$ has a branch point at $s_2 = 0$. Proceeding similarly to the above case, we use the contour in figure B1 and obtain

$$\begin{equation} \mathcal{L}^{-1}_{s_2\to \tau_2}\left\{\widetilde{\mathscr{H}}(y,s_2)\right\} = \frac{1}{2\pi i} \int_{C_1} \mathrm{e}^{s \tau_2} \widetilde{\mathscr{H}}(y,s_2) \mathrm{d}s = \frac{-1}{2\pi i}\left(\int_{C_3}+\int_{C_5}\right) \mathrm{e}^{s_2 \tau_2} \widetilde{\mathscr{H}}(y,s_2) \mathrm{d}s_2, \end{equation} \tag{ B.8 }$$

because in the limit $R\to\infty$ and α → 0 the contributions from C₂, C₆ and C₄ vanish, and we have

$$\begin{equation} \int_{C_3} \mathrm{e}^{s_2 \tau_2} \widetilde{\mathscr{H}}(y,s_2) \mathrm{d}s_2 = -\int_{R}^{\alpha} \mathrm{e}^{- \tau_2 z } \frac{2\left(-i\sqrt{z}+\mathrm{e}^{2i\sqrt{z}}(i\sqrt{z}+2)\right)}{i\sqrt{z}\left(i\sqrt{z}+\mathrm{e}^{2 i\sqrt{z}}(i\sqrt{z}+2)\right)}\mathrm{d}z, \end{equation} \tag{ B.9 }$$

and

$$\begin{equation} \int_{C_5} \mathrm{e}^{s_2 \tau_2} \widetilde{\mathscr{H}}(y,s_2) \mathrm{d}s_2 = -\int_{\alpha}^{R} \mathrm{e}^{- \tau_2 z } \frac{2\left(i\sqrt{z}+\mathrm{e}^{-2i\sqrt{z}}(-i\sqrt{z}+2)\right)}{-i\sqrt{z}\left(-i\sqrt{z}+\mathrm{e}^{-2 i\sqrt{z}}(-i\sqrt{z}+2)\right)}\mathrm{d}z. \end{equation} \tag{ B.10 }$$

After taking $R\to\infty$, α → 0 and substituting these expressions into equation (B.8) one can see that we recover $\mathscr{H}(y,\tau_2)$ in equation (17).

The marginal over t_m of $\mathcal{P}(M,t_m,t|0^+)$ corresponds to $\widetilde{Q}_m(M,0,\epsilon|0^+)$ and from equations (B.2)–(B.4), we are looking for the Laplace inversion, $\mathcal{L}^{-1}_{s \to \tau} \left\{\widetilde{\mathcal{I}}^+(y,s)\right\}$, where

$$\begin{equation} \widetilde{\mathcal{I}}^+(y,s) = \frac{2\mathrm{e}^{y\sqrt{s}}(\sqrt{s}+1)(-\sqrt{s}+\mathrm{e}^{2y\sqrt{s}}(\sqrt{s}+2))}{\sqrt{s}(\sqrt{s}+\mathrm{e}^{2y\sqrt{s}}(\sqrt{s}+2))^2}. \end{equation} \tag{ C.1 }$$

Once again we have a branch point at the origin and thus use the contour C in figure B1 to perform the Laplace inversion. The contributions from C₂, C₆ and C₄ vanish, leaving us to calculate $\int_{C_3}$ and $\int_{C_5}$. By using the substitutions, $s = z\mathrm{e}^{i\pi}$ and $s = z\mathrm{e}^{-i\pi}$ for C₃ and C₅ respectively, whilst taking $R\to \infty$ and α → 0, we have

$$\begin{equation} \mathcal{L}^{-1}_{s\to\tau}\left\{\widetilde{\mathcal{I}}^+(y,s) \right\} = \frac{-1}{2\pi i}\left(\int_{C_3}+\int_{C_5}\right) \mathrm{e}^{s \tau} \widetilde{\mathcal{I}}^+(y,s) \mathrm{d}s = \frac{1}{\pi}\int_0^\infty \frac{\mathrm{e}^{-\tau z}}{\sqrt{z}}j^+(y,z)h^2(y,z)\mathrm{d}z, \end{equation} \tag{ C.2 }$$

where

$$\begin{equation} j^+(y,z) = (5z+4)\cos(y\sqrt{z})-z\cos(3y\sqrt{z})+8\sqrt{z}\sin^3(y\sqrt{z}). \end{equation} \tag{ C.3 }$$

Similarly, since $\widetilde{\mathcal{I}}^-(y,s) = (\sqrt{s}+1)^{-1}\widetilde{\mathcal{I}}^+(y,s)+(\sqrt{s}+s)^{-1}\delta(y)$, we find

$$\begin{align} \mathcal{L}^{-1}_{s\to\tau}\left\{\widetilde{\mathcal{I}}^-(y,s) \right\}& = \frac{-1}{2\pi i}\left(\int_{C_3}+\int_{C_5}\right) \mathrm{e}^{s \tau} \widetilde{\mathcal{I}}^-(y,s) \mathrm{d}s\\ \nonumber \end{align} \tag{ C.4 }$$

$$\begin{align} & = \frac{1}{\pi}\int_0^\infty \frac{\mathrm{e}^{-\tau z}}{\sqrt{z}}j^-(y,z)h^2(y,z)\mathrm{d}z + \mathrm{e}^\tau\mathrm{erfc}(\tau)\delta(y), \end{align}$$

where,

$$\begin{equation} \quad j^-(y,z) = (4-z)\cos(y\sqrt{z})-3z\cos(3y\sqrt{z})+2\sqrt{z}\left[z+(z-2)\cos(2y\sqrt{z})\right]\sin(y\sqrt{z}). \end{equation} \tag{ C.5 }$$

Starting from equations (45) and (46), we may write

$$\begin{equation} \widetilde{Q}_r(p,\epsilon|0^\pm) = \frac{D}{\kappa^2}\widetilde{\mathscr{F}}_r^\pm\left(\frac{D}{\kappa^2}(p+\epsilon),\frac{D}{\kappa^2}\epsilon\right), \end{equation} \tag{ D.1 }$$

therefore

$$\begin{equation} \mathcal{P}(t_r,t|0^\pm) = \frac{\kappa^2}{D}\mathcal{L}^{-1}_{s_2\to \tau_2}\left\{\mathcal{L}^{-1}_{s_1\to \tau_1}\left\{\widetilde{\mathscr{F}}_r^\pm(s_1,s_2) \right\} \right\}, \end{equation} \tag{ D.2 }$$

for $\tau_1 = \kappa^2 t_r/D$ and $\tau_2 = \kappa^2 (t-t_r)/D$, where

$$\begin{equation} \widetilde{\mathscr{F}}_r^+(s_1,s_2) = \frac{s_2+\sqrt{s_1}+\sqrt{s_2}}{s_1 s_2+\sqrt{s_2} \left(\sqrt{s_1 s_2}+s_1\right)} \end{equation} \tag{ D.3 }$$

and

$$\begin{equation} \widetilde{\mathscr{F}}_r^-(s_1,s_2) = \frac{s_1+\sqrt{s_1}+\sqrt{s_2}}{s_1 s_2+\sqrt{s_2} \left(\sqrt{s_1 s_2}+s_1\right)}. \end{equation} \tag{ D.4 }$$

We now proceed to compute the Laplace inversions, $s_1\to \tau_1$ and $s_2\to\tau_2$, of equations (D.3) and (D.4).

D.1. Double Laplace inversion of $\widetilde{\mathscr{F}}_r^+(s_1,s_2)$

Let us first write equation (D.3) as

$$\begin{equation} \widetilde{\mathscr{F}}_r^+(s_1,s_2) = \frac{1}{\sqrt{s_1} \left(s_2+\sqrt{s_2}\right)+s_2}+\frac{s_2+\sqrt{s_2}}{s_1 \left(s_2+\sqrt{s_2}\right)+\sqrt{s_1} s_2}, \end{equation} \tag{ D.5 }$$

and using the following inverse Laplace transform relations [60]:

$$\begin{equation} \mathcal{L}^{-1}_{s_1\to\tau_1}\left\{\frac{1}{\sqrt{s_1}+a}\right\} = \frac{1}{\sqrt{\pi \tau_1 } }-a \mathrm{e}^{a^2 \tau_1 } \mathrm{erfc}\left(a \sqrt{\tau_1 },\right) \end{equation} \tag{ D.6 }$$

and

$$\begin{equation} \mathcal{L}^{-1}_{s_1\to\tau_1}\left\{\frac{1}{\sqrt{s_1}(\sqrt{s_1}+a)}\right\} = \mathrm{e}^{a^2 \tau_1 } \mathrm{erfc}\left(a \sqrt{\tau_1 }\right), \end{equation} \tag{ D.7 }$$

we obtain

$$\begin{align} & \mathcal{L}^{-1}_{s_1\to\tau_1}\left\{\widetilde{\mathscr{F}}_r^+(s_1,s_2)\right\} = \frac{1}{\sqrt{\pi \tau_1} \left(s_2+\sqrt{s_2}\right) } + \left(1-\frac{1}{\left(\sqrt{s_2}+1\right)^2}\right)\\ \nonumber \end{align} \tag{ D.8 }$$

$$\begin{align} & \quad \,\,\,\,\,\qquad\qquad\qquad\qquad\times e^{\frac{s_2 \tau_1}{\left(\sqrt{s_2}+1\right)^2}} \mathrm{erfc}\left(\frac{\sqrt{s_2 \tau_1}}{\sqrt{s_2}+1}\right). \end{align}$$

Taking the inverse Laplace transform of the first term on the left-hand side (LHS) of equation (D.8) gives equation (50). To take the inverse Laplace transform of the second term we use the following series representation,

$$\begin{equation} e^{z^2}\mathrm{erfc}(z) = \sum_{n = 0}^{\infty} \frac{(-1)^n z^n}{\Gamma\left(\frac{n}{2}+1\right)}, \end{equation} \tag{ D.9 }$$

where we then have

$$\begin{align} & \left(1-\frac{1}{\left(\sqrt{s_2}+1\right)^2}\right)e^{\frac{s_2 \tau_1}{\left(\sqrt{s_2}+1\right)^2}} \mathrm{erfc}\left(\frac{\sqrt{s_2 \tau_1}}{\sqrt{s_2}+1}\right)\\[8pt] \nonumber &\quad = \left(1+\frac{\widetilde{\zeta}_1(\sqrt{s_2})}{\sqrt{s_2}}\right)\sum_{n = 0}^{\infty} \frac{(-1)^n\tau_1^{n/2}}{\Gamma\left(\frac{n}{2}+1\right)} \widetilde{\zeta}_{n+1}(\sqrt{s_2}), \end{align} \tag{ D.10 }$$

where

$$\begin{equation} \widetilde{\zeta}_n(\sqrt{s_2}) = \left(\frac{\sqrt{s_2}}{1+\sqrt{s_2}}\right)^n. \end{equation} \tag{ D.11 }$$

We now proceed to calculate the inversion, $\mathcal{L}^{-1}_{s_2\to\tau_2}\left\{\widetilde{\zeta}_n(\sqrt{s_2})\right\}$, where we first find the inversion of $\mathcal{L}^{-1}_{s_2\to\tau_2}\left\{\widetilde{\zeta}_n(s_2)\right\}$. We do this by utilizing the following property,

$$\begin{equation} \mathcal{L}_{\tau \to s}\left\{\frac{\mathrm{d}^n}{\mathrm{d} \tau^n} f(\tau)\right\} = s^n \widetilde{f}(s)-\sum_{m = 1}^n s^{n-m} \lim_{\tau\to0} \frac{\mathrm{d}^{m-1}}{\mathrm{d} \tau^{m-1}}f(\tau), \end{equation} \tag{ D.12 }$$

for some arbitrary function $f(\tau)$. Using $\mathcal{L}^{-1}_{s\to \tau}\{(1+s)^{-n} \} = \tau^{n-1}\mathrm{e}^{-\tau}/\Gamma(n)$, and from the generalized product rule we have

$$\begin{equation} \frac{\mathrm{d}^n}{\mathrm{d} \tau^n}\frac{\tau^{n-1}\mathrm{e}^{-\tau}}{\Gamma(n)} = \mathrm{e}^{-\tau}\sum_{k = 0}^n (-1)^k {n\choose k} \frac{\tau^{k-1}}{\Gamma(k)} = -n {}_1F_1(n+1,2,-\tau), \end{equation} \tag{ D.13 }$$

where ${}_1F_1(a,b,z)$ is the Kummer confluent hypergeometric function [67]. From, $\frac{\mathrm{d}^{m-1}}{\mathrm{d} \tau^{m-1}}\tau^{n-1}\mathrm{e}^{-\tau}/\Gamma(n) = \tau^{n-m}{}_1F_1(n,n+m-1,-\tau)$, and since for τ → 0 this is only non-zero for m = n, we have,

$$\begin{equation} \mathcal{L}^{-1}_{s_2\to\tau_2}\left\{\widetilde{\zeta}_n(s_2)\right\} = -n{}_1F_1(n+1,2,-\tau_2)+\delta(\tau_2). \end{equation} \tag{ D.14 }$$

To find $\mathcal{L}^{-1}_{s_2\to\tau_2}\left\{\widetilde{\zeta}_n(\sqrt{s_2})\right\}$ we use the property [60]:

$$\begin{equation} \mathcal{L}^{-1}_{s\to\tau}\left\{\widetilde{f}(\sqrt{s})\right\} = \frac{1}{\sqrt{4\pi \tau^3}}\int_0^\infty u \mathrm{e}^{-\frac{u^2}{4\tau}} f(u)\mathrm{d}u, \end{equation} \tag{ D.15 }$$

and calculate the following integral as

$$\begin{align} \frac{1}{\sqrt{4 \pi \tau ^3}}\int_0^{\infty } u \mathrm{e}^{-\frac{u^2}{4 \tau }} \, _1F_1(n+1;2;-u) \, \mathrm{d}u& = \frac{1}{\sqrt{\pi \tau }}\, _2F_2\left(\frac{n}{2}+\frac{1}{2},\frac{n}{2}+1;\frac{1}{2},\frac{3}{2};\tau \right)\\[8pt] \nonumber &\quad -\frac{1}{2} (n+1) \, _2F_2\left(\frac{n}{2}+1,\frac{n}{2}+\frac{3}{2};\frac{3}{2},2;\tau \right), \end{align} \tag{ D.16 }$$

then

$$\begin{align} \mathcal{L}^{-1}_{s_2\to\tau_2}\left\{\widetilde{\zeta}_n(\sqrt{s_2})\right\}& = \frac{n}{2} (n+1) \, _2F_2\left(\frac{n}{2}+1,\frac{n}{2}+\frac{3}{2};\frac{3}{2},2;\tau_2 \right)\\[8pt]\nonumber & \quad -\frac{n}{\sqrt{\pi \tau_2 }}\, _2F_2\left(\frac{n}{2}+\frac{1}{2},\frac{n}{2}+1;\frac{1}{2},\frac{3}{2};\tau_2 \right)+\delta(\tau_2), \end{align} \tag{ D.17 }$$

where

$$\begin{equation} \, _2F_2(a,b;c,d,z) = \sum _{k = 0}^{\infty } \frac{z^k \left(a_1\right)_k \left(a_2\right)_k}{k! \left(b_1\right)_k \left(b_2\right)_k}, \end{equation} \tag{ D.18 }$$

with $(x)_n = \Gamma(x+n)/\Gamma(x)$ [67].

From equation (D.10) we also require, $\mathcal{L}^{-1}_{s_2\to\tau_2}\left\{\widetilde{\zeta}_n(\sqrt{s_2})s_2^{-1/2}\right\}$, where we use the following inversion [60]:

$$\begin{equation} \mathcal{L}^{-1}_{s\to\tau}\left\{\frac{s^{n-1}}{(s+1)^n}\right\} = {}_1F_1(n,1,-\tau) \end{equation} \tag{ D.19 }$$

with equation (D.15) to find,

$$\begin{align} \mathcal{L}^{-1}_{s_2\to\tau_2}\left\{\widetilde{\zeta}_n(\sqrt{s_2})s_2^{-1/2}\right\}& = \frac{1}{\sqrt{4 \pi \tau_2 ^3}}\int_0^{\infty } u \mathrm{e}^ {-\frac{u^2}{4 \tau_2 }} \, _1F_1(n;1;-u) \, \mathrm{d}u\\[8pt] \nonumber & = \frac{1}{\sqrt{\pi \tau_2} }\, _2F_2\left(\!\frac{n}{2},\frac{n}{2}+\frac{1}{2};\frac{1}{2},\frac{1}{2};\tau_2 \!\right)-n \, _2F_2\left(\!\frac{n}{2}+\frac{1}{2},\frac{n}{2}+1;1,\frac{3}{2};\tau_2 \!\right). \end{align} \tag{ D.20 }$$

Putting all of this together we have

$$\begin{align} \mathscr{F}_r^+(\tau_1,\tau_2) = \frac{1}{\sqrt{\pi \tau_1}}e^{\tau_2}\mathrm{erfc}(\sqrt{\tau_2})+\sum_{n = 0}^{\infty}\frac{(-1)^n \tau_1^{n/2}}{\Gamma(\frac{n}{2}+1)} \mathcal{C}_n^+(\tau_2)+e^{\tau_1}\mathrm{erfc}(\sqrt{\tau_1})\delta(\tau_2) \end{align} \tag{ D.21 }$$

where

$$\begin{align} \mathcal{C}_n^+(\tau_2)& = (n+2)\Bigg[\frac{(n+1)}{2}\, _2F_2\left(\frac{n}{2}+\frac{3}{2},\frac{n}{2}+2;\frac{3}{2},2;\tau_2 \right)\notag\\[8pt] & \quad -\, _2F_2\left(\frac{n}{2}+\frac{3}{2},\frac{n}{2}+2;1,\frac{3}{2};\tau_2 \right)\Bigg]+\frac{1}{\sqrt{\pi \tau_2}}\Bigg[\, _2F_2\left(\frac{n}{2}+1,\frac{n}{2}+\frac{3}{2};\frac{1}{2},\frac{1}{2};\tau_2 \right)\notag\\[8pt] & \quad -(n+1)\, _2F_2\left(\frac{n}{2}+1,\frac{n}{2}+\frac{3}{2};\frac{1}{2},\frac{3}{2};\tau_2 \right)\Bigg]. \end{align} \tag{ D.22 }$$

D.2. Double laplace inversion of $\widetilde{\mathscr{F}}_r^-(s_1,s_2)$

We write $\widetilde{\mathscr{F}}_r^-(s_1,s_2)$ as

$$\begin{equation} \widetilde{\mathscr{F}}_r^-(s_1,s_2) = \frac{\sqrt{s_1}+1}{\sqrt{s_1} s_2+\sqrt{s_2} \left(\sqrt{s_1}+\sqrt{s_2}\right)}+\frac{\sqrt{s_2}}{\sqrt{s_1} s_2+s_1 \left(s_2+\sqrt{s_2}\right)}, \end{equation} \tag{ D.23 }$$

then using equations (D.6) and (D.7) with (D.12), i.e.

$$\begin{equation} \mathcal{L}^{-1}_{s_1\to\tau_1}\left\{\frac{\sqrt{s_1}}{{\sqrt{s_1}+a}}\right\} = -\frac{a}{\sqrt{\pi \tau_1}}+a^2\mathrm{e}^{a^2\tau_1}\mathrm{erfc}\left(a\sqrt{\tau_1}\right) + \delta(\tau_1), \end{equation} \tag{ D.24 }$$

we obtain,

$$\begin{align} \mathcal{L}^{-1}_{s_1\to\tau_1}\left\{\widetilde{\mathscr{F}}_r^-(s_1,s_2)\right\}& = \frac{1}{\sqrt{\pi } \left(\sqrt{s_2}+1\right)^2 \sqrt{s_2 \tau_1}}+\left(1-\frac{1}{\left(\sqrt{s_2}+1\right)^2}\right)\\ \nonumber & \quad \times \frac{e^{\frac{s_2 \tau_1}{\left(\sqrt{s_2}+1\right)^2}} \mathrm{erfc}\left(\frac{\sqrt{s_2 \tau_1}}{\sqrt{s_2}+1}\right)}{\left(\sqrt{s_2}+1\right)} +\frac{1}{\sqrt{s_2}+s_2}\delta(\tau_1). \end{align} \tag{ D.25 }$$

The Laplace inversion of the first term on the RHS of equation (D.25) corresponds to equation (51). One can see the second term in equation (D.25) is only different by a factor of $(1+\sqrt{s_2})^{-1}$ to the second term on the RHS in equation (D.8), then by using equation (D.9), we have

$$\begin{align} \left(1-\frac{1}{\left(\sqrt{s_2}+1\right)^2}\right)\frac{e^{\frac{s_2 \tau_1}{\left(\sqrt{s_2}+1\right)^2}} \mathrm{erfc}\left(\frac{\sqrt{s_2 \tau_1}}{\sqrt{s_2}+1}\right)}{\left(\sqrt{s_2}+1\right)}& = \left(\frac{1}{\sqrt{s_2}}+\frac{\widetilde{\zeta}_1(s_2)}{s_2}\right)\\ \nonumber &\quad\times \sum_{n = 0}^{\infty} \frac{(-1)^n\tau_1^{n/2}}{\Gamma\left(\frac{n}{2}+1\right)} \widetilde{\zeta}_{n+2}(\sqrt{s_2}). \end{align} \tag{ D.26 }$$

While the inversion $\mathcal{L}^{-1}_{s_2\to\tau_2}\left\{\widetilde{\zeta}_n(\sqrt{s_2})s_2^{-1/2}\right\}$ is given by equation (D.20), we find $\mathcal{L}^{-1}_{s_2\to\tau_2}\left\{\widetilde{\zeta}_n(\sqrt{s_2})/s_2\right\}$ by integrating equation (D.17) over τ₂ to find,

$$\begin{align} \mathcal{L}^{-1}_{s_2\to \tau_2}\left\{\frac{\widetilde{\zeta}_n(\sqrt{s_2})}{s_2}\right\} & = \, _2F_2\left(\frac{n}{2},\frac{n}{2}+\frac{1}{2};\frac{1}{2},1;\tau _2\right)\\\nonumber &\quad -\frac{2 n \sqrt{\tau _2}}{\sqrt{\pi }} \, _2F_2\left(\frac{n}{2}+\frac{1}{2},\frac{n}{2}+1;\frac{3}{2},\frac{3}{2};\tau _2\right). \end{align} \tag{ D.27 }$$

By combining all the above we find $\mathscr{F}_r^-(\tau_1,\tau_2)$ to be,

$$\begin{align} \mathscr{F}_r^-(\tau_1,\tau_2) = \frac{2}{\pi}\sqrt{\frac{\tau_2}{\tau_1}}-\frac{2\tau_2}{\sqrt{\pi \tau_1}}\mathrm{e}^{\tau_2}\,\,\,\mathrm{erfc}(\sqrt{\tau_2})+\sum_{n = 0}^{\infty}\frac{(-1)^n \tau_1^{n/2}}{\Gamma(\frac{n}{2}+1)} \mathcal{C}_n^-(\tau_2)+\mathrm{e}^{\tau_2}\,\,\,\mathrm{erfc}(\sqrt{\tau_2})\delta(\tau_1),\nonumber\\ \end{align} \tag{ D.28 }$$

where

$$\begin{align} \mathcal{C}_n^-(\tau_2) & = \, _2F_2\left(\frac{n}{2}+\frac{3}{2},\frac{n}{2}+2;\frac{1}{2},1;\tau_2\right)+\frac{1}{\sqrt{\pi \tau_2}}\, _2F_2\left(\frac{n}{2}+1,\frac{n}{2}+\frac{3}{2};\frac{1}{2},\frac{1}{2};\tau_2\right)\\ \nonumber & \quad -\frac{2(n+3)\sqrt{\tau_2}}{\sqrt{\pi}}\, _2F_2\left(\frac{n}{2}+2,\frac{n}{2}+\frac{5}{2};\frac{3}{2},\frac{3}{2};\tau_2\right)-(n+2)\, _2F_2\left(\frac{n}{2}+\frac{3}{2},\frac{n}{2}+2;1,\frac{3}{2};\tau_2\right). \end{align} \tag{ D.29 }$$

Using $\lim_{z\to0}{}_2F_2(a,b;c,d;z) = 1$ and the definition of $\mathcal{C}^+(\tau_1,\tau_2)$ and $\mathcal{C}^-(\tau_1,\tau_2)$ in equations (D.22) and (D.29) for $\tau_1,\tau_2\to 0$ we find that, $\unicode{x210A}^\pm(\tau_1,\tau_2)$ becomes,

$$\begin{equation} \unicode{x210A}^+(\tau_1,\tau_2)\simeq \sum_{n = 1}^{\infty} \frac{(-1)^{n+1}n \tau_1^{n/2}}{\Gamma\left(1+\frac{n}{2}\right)\sqrt{\pi \tau_2}} \end{equation} \tag{ E.1 }$$

and

$$\begin{equation} \unicode{x210A}^-(\tau_1,\tau_2)\simeq \sum_{n = 0}^{\infty} \frac{(-1)^{n}\tau_1^{n/2}}{\Gamma\left(1+\frac{n}{2}\right)\sqrt{\pi \tau_2}}, \end{equation} \tag{ E.2 }$$

and for $\tau_1\to 0$ the sums in equations (E.1) and (E.2) being dominated by the first term, reduce to

$$\begin{equation} \unicode{x210A}^+(\tau_1,\tau_2)\simeq \frac{2}{\pi}\sqrt{\frac{\tau_1}{\tau_2}}, \end{equation} \tag{ E.3 }$$

and

$$\begin{equation} \unicode{x210A}^-(\tau_1,\tau_2)\simeq \frac{1}{\sqrt{\pi \tau_2}}. \end{equation} \tag{ E.4 }$$