Generalized master equation for first-passage problems in partitioned spaces

Before introducing the general framework, we formulate the problem for the case of two domains amended by an absorbing target, which serves as an illustration and a test bed of the proposed GME to calculate the FPT distribution for target search. We consider free diffusion (i) on a 1D half axis with an absorbing boundary at position R > 0, and (ii) in the 3D space between two concentric spheres of radii R < b, where the boundary of the inner sphere is absorbing and the outer one is reflecting (figures 1(a) and (b)). A domain boundary is placed at the position a > R (1D) and at radius a with b > a > R (3D), respectively, and we refer to the region between R and a as inner domain and to the remaining accessible space as outer domain. The trajectories start in the outer domain at x₀ with |x₀| > a. (For simplicity, we use the vector notation also in 1D space.) The central quantity of interest is the probability density p_FPT(t) of the random time t until the first encounter with the boundary at R.

The trajectories are dissected into partial ones whenever the boundary at |x| = a is hit. By the Markov property of diffusion, the evolution of a partial trajectory inside a domain depends on the entry point at the boundary, but not on the motion in the other domain. The partial trajectory ends when reaching another point of the domain boundary. The dwell time on a given partial trajectory (equivalently, in its respective domain) is the time the molecule spends between two consecutive events of hitting the domain boundary, which include entrance to and exit from the domain, but also merely touching the boundary. The probability density of dwell times is an FPT density itself, which is why it will be referred to as a partial FPT density in the following. In the outer domain, there are two types of trajectories: starting at x₀ and starting at the boundary |x| = a, both types end at this boundary; the corresponding partial FPT densities are ϕ₀(t) and ϕ_out(t). Trajectories in the inner domain always start at |x| = a, but either leave the domain at this boundary again or are stopped at |x| = R with partial FPT densities ϕ_in(t) and ${\phi }_{\mathtt{\text{in}}}^{\varnothing }\left(t\right)$, respectively. The specific forms of these dwell time densities for the geometries chosen here are discussed in sections 3.1 and 4.1

As a first test, we assume no physical difference between the inner and outer domains so that the boundary at |x| = a is only an artificial one. This allows us to check the proposed GME approach of dissecting trajectories and assembling the overall FPT density ϕ_FPT(t) for reaching the boundary |x| = R from the partial ones. The result must reproduce the known distribution of FPTs for reaching |x| = R directly when starting at x₀. In section 5, we elaborate on the situation of different diffusion coefficients in the two domains.

The derivation of the GME of the present problem follows ideas in reference [36] and makes use of two types of conservation laws: (1). Probability is conserved locally: net gains or losses within one state determine the change of local probability, and (2). Probability is conserved in the transitions between the states (continuity of the fluxes). The other key ingredient is an equation accounting for the renewal character of the stochastic process, linking present losses from a state with the gains at a previous time.

For the general scheme, we consider a partition of the accessible space into n domains labeled by α ∈ {1, ..., n}. The trajectory of a diffusing molecule is dissected whenever it hits a boundary between domains, and the label of the corresponding domain is assigned to each partial trajectory. This gives rise to a set of coarse-grained states and to the probabilities ρ_α (t) that the molecule is found in domain α at time t. If two domains α and β share a boundary α|β, then transitions across this boundary induce a probability flux. Accounting for the direction of the transition, j_α|β (t) denotes the flux from state β to state α at time t. Transitions from a domain to itself are possible (j_α|α ≠ 0) since at the boundary of a domain the origin of the trajectory is irrelevant by the assumed renewal property. The overall probability gain of state α is ${j}_{\alpha }^{+}={\sum }_{\beta }\;{j}_{\alpha \vert \beta }$, whereas the total loss flux reads ${j}_{\alpha }^{-}={\sum }_{\beta }\;{j}_{\beta \vert \alpha }$. Local conservation of probability [principle (i)] then implies

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}t}{\rho }_{\alpha }\left(t\right)={j}_{\alpha }^{+}\left(t\right)-{j}_{\alpha }^{-}\left(t\right)=\sum\limits _{\beta }\left[{j}_{\alpha \vert \beta }\left(t\right)-{j}_{\beta \vert \alpha }\left(t\right)\right],\end{equation} \tag{ 1 }$$

that is, the temporal change of probability in a state is the difference between the total gain and loss fluxes. Summing over α shows that the overall probability is conserved globally, as it should: (d/dt)∑_α ρ_α (t) = 0.

The continuity of the fluxes, principle (ii), determines the total gain fluxes as a fixed linear combination of the total losses,

$$\begin{equation}{j}_{\alpha }^{+}=\sum\limits _{\beta }{w}_{\alpha \vert \beta }{j}_{\beta }^{-},\end{equation} \tag{ 2 }$$

where the weights w_α|β encode the connectivity of the domains and form a stochastic n × n matrix with columns adding up to unity, ∑_α w_α|β = 1. The sum in equation (2) is actually restricted to those domains β that are adjacent to α, otherwise we can put w_α|β = 0. Each summand represents the partial gain of state α stemming from a loss of probability in state β. Thus, the probability flux across the boundary α|β directed toward α is given by the product

$$\begin{equation}{j}_{\alpha \vert \beta }={w}_{\alpha \vert \beta }{j}_{\beta }^{-}.\end{equation} \tag{ 3 }$$

The w_α|β include the transmission probability q_α|β that a partial trajectory ending at the α|β boundary is indeed continued in the domain α. The fraction of the loss flux ${j}_{\beta }^{-}$ that reaches the boundary α|β is w_α|β /q_α|β and sums to unity, ∑_α≠β w_α|β /q_α|β = 1. The flux ${j}_{\beta \vert \beta }={w}_{\beta \vert \beta }{j}_{\beta }^{-}$ describes those trajectories that hit the boundary of domain β, but return and are continued inside of β. If the transition probabilities at all boundaries of β are equal, q_α|β ≕ q_β , the fraction of such 'remainers' is

$$\begin{equation}{w}_{\beta \vert \beta }=1-\sum\limits _{\alpha \ne \beta }{w}_{\alpha \vert \beta }=1-{q}_{\beta }.\end{equation} \tag{ 4 }$$

For unbiased diffusion, q_α|β = 1/2 at all boundaries α|β, so that w_β|β = 1/2 is the probability of a trajectory to remain in the domain after reaching its boundary. In the presence of one (or more) absorbing domain ∅, exits from such a domain cannot occur, j_α|∅ = 0 and so ${j}_{\varnothing }^{-}=0$. The temporal change of ρ_∅(t) is the sought FPT density of the target search problem, ${p}_{\text{FPT}}\left(t\right)=\mathrm{d}{\rho }_{\varnothing }\left(t\right)/\mathrm{d}t={j}_{\varnothing }^{+}\left(t\right)$.

Invoking the picture of an ensemble of particles, a loss of particles from a domain at time t can only happen if the particles had been there before, either from the very beginning or through gains at an earlier time t − τ, where τ is the dwell time of a specific particle (i.e., a partial trajectory) in that domain. The loss fluxes are thus linked to the gain fluxes and obey renewal-like relations [36], which for a domain β reads schematically:

$$\begin{equation}{j}_{\beta }^{-}\left(t\right)={\int }_{0}^{t}{\phi }_{\beta }\left(\tau \right)\enspace {j}_{\beta }^{+}\left(t-\tau \right)\enspace \mathrm{d}\tau +{\phi }_{\beta }^{\left(0\right)}\left(t\right){\rho }_{\beta }^{\left(0\right)}\end{equation} \tag{ 5 }$$

with ϕ_β (τ) denoting the probability density of dwell times. The last term describes particles that were initially placed inside of the domain with probability ${\rho }_{\beta }^{\left(0\right)}={\rho }_{\beta }\left(t=0\right)$ and reach the domain boundary at time t according to the FPT density ${\phi }_{\beta }^{\left(0\right)}\left(t\right)$. Equation (5) applies only under certain conditions, e.g., for a highly symmetric domain, as shown in appendix A.

More generally, we distinguish the parts of the boundary of domain β according to their adjacent domain and consider the 'ports' α|β that allow for transitions to and from a domain α. Then, the dwell time of a partial trajectory in β starting at the boundary β|γ and stopping at α|β is statistically characterized by its FPT density ${\phi }_{\beta \vert \gamma }^{\alpha }\left(t\right)$, with α and γ taken from the set of adjacent domains. The probability of leaving the domain β through the boundary α|β after starting at β|γ is referred to as splitting probability,

$$\begin{equation}{P}_{\beta \vert \gamma }^{\alpha }={\int }_{0}^{\infty }{\phi }_{\beta \vert \gamma }^{\alpha }\left(t\right)\enspace \mathrm{d}t;\end{equation} \tag{ 6 }$$

the splitting probabilities of the same initial boundary sum up to unity, ${\sum }_{\alpha \ne \beta }\;{P}_{\beta \vert \gamma }^{\alpha }=1$

For the directed probability fluxes j_α|β through these ports, we postulate the following renewal equation as a generalization of equation (5):

$$\begin{equation}{j}_{\alpha \vert \beta }={q}_{\alpha \vert \beta }\sum\limits _{\gamma \ne \beta }{\phi }_{\beta \vert \gamma }^{\alpha }{\ast}\left[{j}_{\beta \vert \gamma }+\left({q}_{\gamma \vert \beta }^{-1}-1\right){j}_{\gamma \vert \beta }\right]+{q}_{\alpha \vert \beta }\enspace {\phi }_{\beta \vert 0}^{\alpha }\enspace {\rho }_{\beta }^{\left(0\right)}\enspace \quad \alpha \ne \beta ,\end{equation} \tag{ 7 }$$

omitting the time arguments and introducing the convolution symbol as $\left(f{\ast}g\right)\left(t\right)={\int }_{0}^{t}f\left(\tau \right)\enspace g\left(t-\tau \right)\enspace \mathrm{d}\tau$ for integrable functions f, g. In the last term of equation (7), ${\phi }_{\beta \vert 0}^{\alpha }\left(t\right)$ is the FPT density of trajectories starting initially in the interior of domain β with probability ${\rho }_{\beta }^{\left(0\right)}{:=}{\rho }_{\beta }\left(0\right)$ and hitting the boundary α|β at time t; the initial position is either fixed at some point r₀ within the domain, or ${\phi }_{\beta \vert 0}^{\alpha }\left(t\right)$ is an average over a distribution of initial positions. The multiplication of the rhs of the renewal equation by the transmission probability q_α|β takes into account that the flux j_α|β includes only those trajectories that are continued in the domain α, but the FPT densities describe only the arrival at the boundary α|β. Concerning the start point of the partial trajectories, there are two possibilities to reach the boundary β|γ, which are expressed by the brackets inside of the convolution: (i) entering the domain β by a transition from γ and (ii) touching the boundary from inside of β, i.e., remainers that were in β before time t − τ. Whereas the flux corresponding to (i) is simply j_β|γ (t − τ), the second situation requires that the flux j_γ|β (t − τ)/q_γ|β of trajectories that reached the boundary from inside is multiplied by the probability 1 − q_γ|β of not leaving to domain γ; for unbiased diffusion, $\left({q}_{\gamma \vert \beta }^{-1}-1\right)=1$. The flux ${j}_{\beta \vert \beta }={w}_{\beta \vert \beta }{j}_{\beta }^{-}$ due to self-transition is proportional to the overall flux to the neighboring domains, where the prefactor follows with ${j}_{\beta }^{-}={\sum }_{\alpha }\;{j}_{\alpha \vert \beta }$, so that

$$\begin{equation}{j}_{\beta \vert \beta }=\frac{{w}_{\beta \vert \beta }}{1-{w}_{\beta \vert \beta }}\sum\limits _{\alpha \ne \beta }{j}_{\alpha \vert \beta };\end{equation} \tag{ 8 }$$

the prefactor is unity if w_β|β = 1/2 as for unbiased diffusion.

The renewal equation equation (7) simplifies considerably if, due to symmetries of the domain β, all boundaries are equivalent with respect to the partial FPT densities; examples are a slab-like domain delimited by two parallel, infinite planes, and a domain that forms a regular simplex. Then, there are only two types of FPT densities, namely ${\phi }_{\beta }^{\left(d\right)}={\phi }_{\beta \vert \gamma }^{\alpha }$ for trajectories connecting distinct boundaries and ${\phi }_{\beta }^{\left(s\right)}={\phi }_{\beta \vert \alpha }^{\alpha }$ for self-transitions (α = γ). Summation of equation (7) over the adjacent domains α of β yields equation (5) for the total loss flux ${j}_{\beta }^{-}$ with the FPT density for reaching some part of the boundary of β:

$$\begin{equation}{\phi }_{\beta }={\phi }_{\beta }^{\left(s\right)}+\left({z}_{\beta }-1\right){\phi }_{\beta }^{\left(d\right)},\end{equation} \tag{ 9 }$$

where the coordination number z_β counts the adjacent domains of β (see appendix A). Further, the loss flux is equally distributed over all boundaries, i.e., w_α|β = q_β /z_β for α ≠ β and w_β|β = 1 − q_β .

The set of linear equation (5) together with equation (2) for β = 1, ..., n can be solved for the n loss fluxes ${j}_{\beta }^{-}$, provided that the weight matrix (w_α|β ) is known a priori; the FPT densities are considered model parameters. Alternatively, equations (7) and (8) specify a set of 2m linear equations that can be solved for the 2m directed fluxes j_α|β , assuming that the system contains m boundaries (ports) between domains.

In the two-domain model of section 2.1, the state of the molecule is described by the probabilities ρ_in(t) and ρ_out(t) that at time t it is found in the inner and outer domain, respectively. It is amended by the probability ρ_∅(t) that the trajectory was stopped at time t or earlier: the molecule has reached the target and was removed from the system, e.g., by a chemical reaction. The transitions between the states induce probability fluxes j_α|β with α, β ∈ {in, out}, which lead to the total gain (+) and loss (−) fluxes, ${j}_{\mathtt{\text{in}}}^{{\pm}}$ and ${j}_{\mathtt{\text{out}}}^{{\pm}}$, of the in and out states, respectively. In addition, there is a flux j_∅|in indicating the trajectories that traversed the in state before being stopped at the target (∅). In the subsequent notation, the flux j_∅|in is treated as an additional loss from the in state and is not included in the symbol ${j}_{\mathtt{\text{in}}}^{-}$. Figure 1(c) shows a transition graph of the model along with the loss fluxes.

For an unbiased diffusion, a partial trajectory ending at the domain boundary |x| = a is continued in either the inner or the outer domain with equal probability 1/2. By the Markov property of diffusion, this applies for partial trajectories irrespective of which side of the boundary they belong to (in or out state). Thus, ${j}_{\mathtt{\text{in}}\vert \mathtt{\text{out}}}={j}_{\mathtt{\text{out}}\vert \mathtt{\text{out}}}=\frac{1}{2}{j}_{\mathtt{\text{out}}}^{-}$, and we will not distinguish between j_in|out and j_out|out in the following; correspondingly, ${j}_{\mathtt{\text{in}}\vert \mathtt{\text{in}}}={j}_{\mathtt{\text{out}}\vert \mathtt{\text{in}}}=\frac{1}{2}{j}_{\mathtt{\text{in}}}^{-}$. Therefore, continuity of the fluxes in the transitions requires that (see equation (2))

$$\begin{equation}{j}_{\mathtt{\text{out}}}^{+}\left(t\right)=\frac{1}{2}\left[{j}_{\mathtt{\text{in}}}^{-}\left(t\right)+{j}_{\mathtt{\text{out}}}^{-}\left(t\right)\right]\quad \text{and}\quad {j}_{\mathtt{\text{in}}}^{+}\left(t\right)=\frac{1}{2}\left[{j}_{\mathtt{\text{out}}}^{-}\left(t\right)+{j}_{\mathtt{\text{in}}}^{-}\left(t\right)\right],\end{equation} \tag{ 10 }$$

where the terms on the right represent the incoming arrows of the state given on the left as depicted in figure 1(c); it follows that ${j}_{\mathtt{\text{out}}}^{+}\left(t\right)={j}_{\mathtt{\text{in}}}^{+}\left(t\right)$. Next, local balance in each state demands (equation (1)):

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}t}{\rho }_{\mathtt{\text{out}}}\left(t\right)={j}_{\mathtt{\text{out}}}^{+}\left(t\right)-{j}_{\mathtt{\text{out}}}^{-}\left(t\right),\end{equation} \tag{ 11a }$$

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}t}{\rho }_{\mathtt{\text{in}}}\left(t\right)={j}_{\mathtt{\text{in}}}^{+}\left(t\right)-{j}_{\mathtt{\text{in}}}^{-}\left(t\right)-{j}_{\varnothing \vert \mathtt{\text{in}}}\left(t\right),\end{equation} \tag{ 11b }$$

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}t}{\rho }_{\varnothing }\left(t\right)={j}_{\varnothing \vert \mathtt{\text{in}}}\left(t\right).\end{equation} \tag{ 11c }$$

One confirms readily that the overall probability, including that of the absorbed particles, is conserved: $\left(\mathrm{d}/\mathrm{d}t\right)\left[{\rho }_{\mathtt{\text{in}}}\left(t\right)+{\rho }_{\mathtt{\text{out}}}\left(t\right)+{\rho }_{\varnothing }\left(t\right)\right]=0$. The quantity ρ_∅(t) in equation (11c) collects the trajectories that have stopped up to time t, and its change is thus equal to the sought FPT density of the target search problem, dρ_∅(t)/dt = p_FPT(t). So the task is to compute the flux j_∅|in(t) onto the boundary |x| = R.

Finally, the loss fluxes make a recursion to the gain fluxes and obey a set of renewal-like relations:

$$\begin{equation}{j}_{\mathtt{\text{out}}}^{-}\left(t\right)={\phi }_{0}\left(t\right)+{\int }_{0}^{t}{\phi }_{\mathtt{\text{out}}}\left(t-{t}^{\prime }\right)\enspace {j}_{\mathtt{\text{out}}}^{+}\left({t}^{\prime }\right)\enspace \mathrm{d}{t}^{\prime },\end{equation} \tag{ 12a }$$

$$\begin{equation}{j}_{\mathtt{\text{in}}}^{-}\left(t\right)={\int }_{0}^{t}{\phi }_{\mathtt{\text{in}}}\left(t-{t}^{\prime }\right)\enspace {j}_{\mathtt{\text{in}}}^{+}\left({t}^{\prime }\right)\enspace \mathrm{d}{t}^{\prime },\end{equation} \tag{ 12b }$$

$$\begin{equation}{j}_{\varnothing \vert \mathtt{\text{in}}}\left(t\right)={\int }_{0}^{t}{\phi }_{\mathtt{\text{in}}}^{\varnothing }\left(t-{t}^{\prime }\right)\enspace {j}_{\mathtt{\text{in}}}^{+}\left({t}^{\prime }\right)\enspace \mathrm{d}{t}^{\prime }.\end{equation} \tag{ 12c }$$

These equations follow from equation (5), which is applicable because in the two-domain problem merely two types of partial FPT densities occur: the outer domain permits transitions only to the in state through the boundary at |x| = a and thus the only FPT density in the domain out is ${\phi }_{\mathtt{\text{out}}}{:=}{\phi }_{\mathtt{\text{out}}\vert \mathtt{\text{in}}}^{\mathtt{\text{in}}}$ for trajectories from this boundary to itself, along with the FPT density ϕ₀ for the initial part of the trajectories starting at r₀ in the interior of out. Partial trajectories in the inner domain in can either end at the boundary |x| = a (flux ${j}_{\mathtt{\text{in}}}^{-}$) or at |x| = R (absorption, j_∅|in), see figure 1(c). Further, there are no gains to in from the absorbing boundary (q_∅|in = 1), and we have only the self-transition term corresponding to ${\phi }_{\mathtt{\text{in}}}={\phi }_{\mathtt{\text{in}}\vert \mathtt{\text{out}}}^{\mathtt{\text{out}}}$ in equation (12b). The loss j_∅|in to the absorbed state is a 'distinct' contribution expressed by the FPT density ${\phi }_{\mathtt{\text{in}}}^{\varnothing }{:=}{\phi }_{\mathtt{\text{in}}\vert \mathtt{\text{out}}}^{\varnothing }$. Note that ϕ_in and ${\phi }_{\mathtt{\text{in}}}^{\varnothing }$ are no proper FPT densities in the sense that they do not normalize, the full density of dwell times in the state in is given by their sum, which satisfies ${\int }_{0}^{\infty }\left[{\phi }_{\mathtt{\text{in}}}\left(t\right)+{\phi }_{\mathtt{\text{in}}}^{\varnothing }\left(t\right)\right]\enspace \mathrm{d}t=1$. This fact expresses the splitting of probabilities to follow one or the other type of partial trajectory (i.e., to survive or to be stopped at the end of the current step).

The system of equations (11a)–(11c) and (12a)–(12c), together with equation (10), forms the GME of the first-passage problem with two domains. Using equation (10), the gain fluxes ${j}_{\mathtt{\text{in}}}^{+}$ and ${j}_{\mathtt{\text{out}}}^{+}$ can be expressed in terms of loss fluxes, which leaves us with a system of six linear integro-differential equations in six variables (three loss fluxes and three occupation probabilities). This type of equation system is conveniently solved in Laplace domain.

The Laplace transform $\tilde {f}{:=}\mathcal{L}\left[f\right]$ of a measurable function f(t) on ${\mathbb{R}}_{{\geqslant}0}$ is defined as [34]

$$\begin{equation}\tilde {f}\left(u\right)=\mathcal{L}\left[f\right]\left(u\right){:=}{\int }_{0}^{\infty }{\text{e}}^{-ut}f\left(t\right)\enspace \mathrm{d}t,\end{equation} \tag{ 13 }$$

where the frequency u > 0 is the Laplace variable conjugate to t. For a convolution $\left(f{\ast}g\right)\left(t\right)={\int }_{0}^{t}f\left(t-{t}^{\prime }\right)g\left({t}^{\prime }\right)\mathrm{d}{t}^{\prime }$ it holds $\mathcal{L}\left[f{\ast}g\right]=\mathcal{L}\left[f\right]\mathcal{L}\left[g\right]$, and for the time derivative one has $\mathcal{L}\left[\dot {f}\right]\left(u\right)=u\enspace \mathcal{L}\left[f\right]\left(u\right)-f\left(0\right)$.

Laplace transformation of the self-consistent, linear system of integro-differential equations (11a)–(11c) and (12a)–(12c), yields a closed set of linear, algebraic equations in the probabilities and fluxes with the partial FPT densities as coefficients. Substituting ${j}_{\mathtt{\text{in}}}^{+}$ and ${j}_{\mathtt{\text{out}}}^{+}$ in equations (12a)–(12c) by means of equations (11a) and (11b), we obtain:

$$\begin{equation}{\tilde {j}}_{\mathtt{\text{out}}}^{-}\left(u\right)={\tilde {\phi }}_{0}\left(u\right)+{\tilde {\phi }}_{\mathtt{\text{out}}}\left(u\right)\left[u{\tilde {\rho }}_{\mathtt{\text{out}}}\left(u\right)-1+{\tilde {j}}_{\mathtt{\text{out}}}^{-}\left(u\right)\right],\end{equation} \tag{ 14a }$$

$$\begin{equation}{\tilde {j}}_{\mathtt{\text{in}}}^{-}\left(u\right)={\tilde {\phi }}_{\mathtt{\text{in}}}\left(u\right)\left[u{\tilde {\rho }}_{\mathtt{\text{in}}}\left(u\right)+{\tilde {j}}_{\mathtt{\text{in}}}^{-}\left(u\right)+{\tilde {j}}_{\varnothing \vert \mathtt{\text{in}}}\left(u\right)\right],\end{equation} \tag{ 14b }$$

$$\begin{equation}{\tilde {j}}_{\varnothing \vert \mathtt{\text{in}}}\left(u\right)={\tilde {\phi }}_{\mathtt{\text{in}}}^{\varnothing }\left(u\right)\left[u{\tilde {\rho }}_{\mathtt{\text{in}}}\left(u\right)+{\tilde {j}}_{\mathtt{\text{in}}}^{-}\left(u\right)+{\tilde {j}}_{\varnothing \vert \mathtt{\text{in}}}\left(u\right)\right],\end{equation} \tag{ 14c }$$

where we made use of the initial conditions ρ_out(0) = 1 and ρ_in(0) = ρ_∅(0) = 0. Solving for the loss fluxes, we have:

$$\begin{equation}{\tilde {j}}_{\mathtt{\text{out}}}^{-}\left(u\right)=\frac{{\tilde {\phi }}_{0}\left(u\right)-{\tilde {\phi }}_{\mathtt{\text{out}}}\left(u\right)}{1-{\tilde {\phi }}_{\mathtt{\text{out}}}\left(u\right)}+\frac{u{\tilde {\phi }}_{\mathtt{\text{out}}}\left(u\right)}{1-{\tilde {\phi }}_{\mathtt{\text{out}}}\left(u\right)}{\tilde {\rho }}_{\mathtt{\text{out}}}\left(u\right)\end{equation} \tag{ 15a }$$

$$\begin{equation}{\tilde {j}}_{\mathtt{\text{in}}}^{-}\left(u\right)=\frac{{\tilde {\phi }}_{\mathtt{\text{in}}}\left(u\right)}{1-{\tilde {\phi }}_{\mathtt{\text{in}}}\left(u\right)}\left[u{\tilde {\rho }}_{\mathtt{\text{in}}}\left(u\right)+{\tilde {j}}_{\varnothing \vert \mathtt{\text{in}}}\right]\end{equation} \tag{ 15b }$$

$$\begin{equation}{\tilde {j}}_{\varnothing \vert \mathtt{\text{in}}}\left(u\right)=\frac{{\tilde {\phi }}_{\mathtt{\text{in}}}\left(u\right)}{1-{\tilde {\phi }}_{\mathtt{\text{in}}}^{\varnothing }\left(u\right)}\left[u{\tilde {\rho }}_{\mathtt{\text{in}}}\left(u\right)+{\tilde {j}}_{\mathtt{\text{in}}}^{-}\right].\end{equation} \tag{ 15c }$$

The system of equations in Laplace domain is completed by

$$\begin{equation}u{\tilde {\rho }}_{\mathtt{\text{out}}}\left(u\right)-1=\frac{1}{2}\left[{\tilde {j}}_{\mathtt{\text{in}}}^{-}\left(u\right)-{\tilde {j}}_{\mathtt{\text{out}}}^{-}\left(u\right)\right]\end{equation} \tag{ 16a }$$

$$\begin{equation}u{\tilde {\rho }}_{\mathtt{\text{in}}}\left(u\right)=\frac{1}{2}\left[{\tilde {j}}_{\mathtt{\text{out}}}^{-}\left(u\right)-{\tilde {j}}_{\mathtt{\text{in}}}^{-}\left(u\right)\right]-{\tilde {j}}_{\varnothing \vert \mathtt{\text{in}}}\left(u\right)\end{equation} \tag{ 16b }$$

$$\begin{equation}u{\tilde {\rho }}_{\varnothing }\left(u\right)={\tilde {j}}_{\varnothing \vert \mathtt{\text{in}}}\left(u\right)\end{equation} \tag{ 16c }$$

from equations (11a)–(11c). Solving the linear system for ${\tilde {j}}_{\varnothing \vert \mathtt{\text{in}}}\left(u\right)$ provides us with an explicit expression for the sought FPT density in Laplace domain, ${\tilde {p}}_{\text{FPT}}\left(u\right)={\tilde {j}}_{\varnothing \vert \mathtt{\text{in}}}\left(u\right)$, which is fully specified by the partial FPT densities ${\tilde {\phi }}_{x}\left(u\right)$:

$$\begin{equation}{\tilde {p}}_{\text{FPT}}\left(u\right)=\frac{2{\tilde {\phi }}_{0}\left(u\right){\tilde {\phi }}_{\varnothing }\left(u\right)}{{\tilde {\phi }}_{\mathtt{\text{in}}}\left(u\right)\left[{\tilde {\phi }}_{\varnothing }\left(u\right)-2\right]-2{\tilde {\phi }}_{\mathtt{\text{out}}}\left(u\right)+4};\end{equation} \tag{ 17 }$$

see equation (B1) of the appendix for the solution for all densities and fluxes. As a by-product, our approach also provides the overall sojourn time in a certain state α ∈ {in, out} up to a time t simply by integrating ρ_α (t) over time, which amounts to calculating ${\tilde {\rho }}_{\alpha }\left(u\right)/u$ in the Laplace domain.

The actual FPT density ϕ_FPT(t) in time domain can be obtained from a numerical Laplace backtransform. The procedure is understood best by switching to the characteristic function of the FPTs, given by the one-sided Fourier transform $\hat{\phi }\left(\omega \right){=\int }_{0}^{\infty }{\text{e}}^{\text{i}\omega t}\phi \left(t\right)\mathrm{d}t$, which is well defined for all frequencies $\omega \in \mathbb{R}$ since ϕ(t) is a probability density and thus integrable. It allows for the analytic continuation to the upper complex plane and is connected to the Laplace transform by $\tilde {\phi }\left(u\right)=\hat{\phi }\left(\text{i}u\right)$. The backtransform is uniquely given by a cosine transform:

$$\begin{equation}\phi \left(t\right)=\frac{2}{\pi }{\int }_{0}^{\infty } \mathrm{cos}\left(\omega t\right)\text{Re}\enspace \hat{\phi }\left(\omega \right)\enspace \mathrm{d}\omega ,\end{equation} \tag{ 18 }$$

using that $\text{Re}\enspace \hat{\phi }\left(\omega \right)$ is an even function in ω. For the robust numerical evaluation of the Fourier integral, we used a modified Filon quadrature as developed recently for the back and forth transformation between time correlation functions and their dissipation spectra [39].

The coefficients in equations (15a)–(15c) become singular if any of the FPT densities ${\tilde {\phi }}_{x}\left(u\right)=1$, i.e., if ϕ_x (t) = δ₊(t) is the density of the Dirac measure on the positive reals, ${\mathbb{R}}_{{\geqslant}0}$, supported at t = 0. Unfortunately, we are facing this problem for ϕ_in(t) and ϕ_out(t) with the setup of figure 1: partial trajectories starting at the boundary |x| = a return to that boundary immediately, almost surely, the starting and ending points are the same. This issue is familiar from the first-return problem for diffusion in the continuum, it does not arise for random walks on a lattice.

As a regularization of these singularities, we require a minimum length ɛ > 0 of the partial trajectories in the calculation of partial FPT distributions and we take the limit ɛ → 0 for the final result of the FPT distribution. This can be accomplished by slightly shifting the boundary between in and out for trajectories leaving the in state. Technically, we introduce an auxiliary boundary at |x| = a + ɛ, with a + ɛ < b in the 3D case (figures 2(a) and (b)). This auxiliary boundary is transparent for partial trajectories in the out state, i.e., they start in the outer domain and end at the boundary |x| = a. Partial trajectories starting at |x| = a are assigned to the in state, they either end at |x| = a + ɛ or are stopped at the absorbing boundary |x| = R, whereas the boundary a is invisible to the partial trajectories in the in-state. In the former case, the following partial trajectory belongs to the out state, as do all trajectories starting at a + ɛ, and it ends at |x| = a. Thus, the out state is followed by an in state again and so forth. With the different states being assigned according to the different starting points, a partial trajectory cannot be successed by a partial trajectory in the same state, as opposed to the situation in section 2.3 (Reflections at the boundary |x| = b in the 3D case are included in the partial FPT density of the out state and do not subdivide the partial trajectory further.) The transport properties (e.g., diffusion coefficient) along a partial trajectory are those of the assigned state (or domain), which leads to a hybrid character of the region a < |x| < a + ɛ. Clearly, this interpretation is an approximation to the original diffusion problem, which is restored in the limit ɛ → 0.

Zoom In Zoom Out Reset image size

Our definition of the occupation probabilities ρ_α (t) of the states α ∈ {in, out, ∅} remains unchanged. And as before, the probability fluxes ${j}_{\alpha }^{{\pm}}\left(t\right)$ denote gains (+) and losses (−) of the state α; the loss j_∅|in from the in state to the absorbed state ∅ is not included in ${j}_{\mathtt{\text{in}}}^{-}\left(t\right)$. The resulting transition graph of the amended two-domain GME is given in figure 2. In the terminology of section 2.2, the transmission probabilities at the in|out boundary are q_in|out = q_out|in = 1, i.e., there are no transitions from a state to itself, w_in|in = w_out|out = 0, and the transitions between in and out are thus strictly alternating. The loss fluxes do not split and flux continuity (equation (2)) demands:

$$\begin{equation}{j}_{\mathtt{\text{in}}}^{+}\left(t\right)={j}_{\mathtt{\text{out}}}^{-}\left(t\right)\quad \text{and}\quad {j}_{\mathtt{\text{out}}}^{+}\left(t\right)={j}_{\mathtt{\text{in}}}^{-}\left(t\right),\end{equation} \tag{ 19 }$$

which replaces equation (10). For the local balance in each state, the same equations (11a)–(11c) as previously hold. Also the equations (12a)–(12c) for the fluxes apply without modifications. This set of equations constitutes the GME of the regularized two-domain model.

In the Laplace domain, the three expressions for the loss fluxes, equations (15a)–(15c), carry over as well since their derivation did not rely on equation (10). The linear system is completed by three equations for the occupation probabilities, which follow from equations (11a)–(11c) and (19):

$$\begin{equation}u{\tilde {\rho }}_{\mathtt{\text{out}}}\left(u\right)-1={\tilde {j}}_{\mathtt{\text{in}}}^{-}\left(u\right)-{j}_{\mathtt{\text{out}}}^{-}\left(u\right),\end{equation} \tag{ 20a }$$

$$\begin{equation}u{\tilde {\rho }}_{\mathtt{\text{in}}}\left(u\right)={\tilde {j}}_{\mathtt{\text{out}}}^{-}\left(u\right)-{\tilde {j}}_{\mathtt{\text{in}}}^{-}\left(u\right)-{\tilde {j}}_{\varnothing \vert \mathtt{\text{in}}}\left(u\right),\end{equation} \tag{ 20b }$$

$$\begin{equation}u{\tilde {\rho }}_{\varnothing }\left(u\right)={\tilde {j}}_{\varnothing \vert \mathtt{\text{in}}}\left(u\right).\end{equation} \tag{ 20c }$$

Solving the linear system in six variables, we find the desired FPT density ${\tilde {p}}_{\text{FPT}}\left(u\right)={\tilde {j}}_{\varnothing \vert \mathtt{\text{in}}}\left(u\right)$ in terms of the partial FPT densities:

$$\begin{equation}{\tilde {p}}_{\text{FPT}}\left(u\right)=\frac{{\tilde {\phi }}_{0}\left(u\right)\enspace {\tilde {\phi }}_{\mathtt{\text{in}}}^{\varnothing }}{1-{\tilde {\phi }}_{\mathtt{\text{in}}}\left(u\right)\enspace {\tilde {\phi }}_{\mathtt{\text{out}}}\left(u\right)},\end{equation} \tag{ 21 }$$

which is a main result of this work. Note that the partial FPT densities (except ${\tilde {\phi }}_{0}\left(u\right)$) implicitly depend on the regularization parameter ɛ. The complete solution for all probabilities and fluxes is given in equation (B2) of the appendix.

As mentioned above, the dwell time probability densities ϕ_α (t) are themselves FPT densities, namely of the first passage from their entrance to the exit from the respective regions. We will refer to these regions as 'outer' and 'inner' regions, respectively (see figure 2). These partial FPT densities can be obtained from the (backwards) Fokker–Planck equation for the probability density ψ(x, t) of the particle position,

$$\begin{equation}\frac{\partial }{\partial t}\psi \left(\mathbf{x},t\right)=D\frac{{\partial }^{2}}{\partial {\mathbf{x}}^{2}}\psi \left(\mathbf{x},t\right),\end{equation} \tag{ 22 }$$

with an initial condition and the respective boundary conditions. More precisely, the setup where particles reaching some point in space for the first time are immediately removed from the system, corresponds to the setup in terms of position probability density with an absorbing boundary at that point. The flux into that point x_ξ will be the FPT density to visit that point for the first time,

$$\begin{equation}{\phi }^{\xi }\left({\mathbf{x}}_{\xi },{\mathbf{x}}_{{\Omega}},t;{\mathbf{x}}_{\mathbf{0}}\right)=-D\frac{\partial }{\partial \mathbf{x}}{\left.\psi \left(\mathbf{x},t\right)\right\vert }_{{\mathbf{x}}_{\xi }},\end{equation} \tag{ 23 }$$

where the particles started at x = x₀ and x_Ω is a point at the boundary.

For each region we will solve the respective PDEs in Laplace domain where they take the form of ordinary differential equations linear in x. The frequency u will be kept as a variable, although it takes the role of a parameter during calculations. In a first step, we will solve the pertinent initial-boundary-value-problems more generally for arbitrary starting points x₀ or r₀ and lower and upper boundaries at x_ℓ , x_u or r_ℓ , r_u , respectively and customize them later. A more detailed exposition of the procedure than we are able to give here can be found in [40].

Specializing to 1D spaces, we have the following equation governing the probability to find a particle at some x

$$\begin{equation}D\frac{{\partial }^{2}}{\partial {x}^{2}}\tilde {\psi }\left(x,u\right)-u\tilde {\psi }\left(x,u\right)=-\delta \left(x-{x}_{0}\right),\end{equation} \tag{ 24 }$$

where the right side represents the initial condition of all particles starting at x₀. The fundamental system is $\left\{\mathrm{exp}\left(\sqrt{u/D}\enspace x\right),\mathrm{exp}\left(-\sqrt{u/D}\enspace x\right)\right\},$ i.e. our solutions will be a linear combination of these two exponentials in x and $\sqrt{u/D}$.

3.1.1. Outer region

In the outer region we have an absorbing lower boundary at some x_ℓ and a reflecting upper one at x_u :

$$\begin{equation}\tilde {\psi }\left({x}_{\ell },u\right)=0,\end{equation} \tag{ 25 }$$

$$\begin{equation}\frac{\partial }{\partial x}{\left.\tilde {\psi }\left(x,u\right)\right\vert }_{x={x}_{u}}=0,\end{equation} \tag{ 26 }$$

for which

$$\begin{equation}\tilde {\psi }\left(x,u\right)=\begin{cases}\frac{1}{\sqrt{uD}}\frac{\mathrm{cosh}\left(\sqrt{u/D}\left({x}_{u}-{x}_{0}\right)\right)}{\mathrm{cosh}\left(\sqrt{u/D}\left({x}_{u}-{x}_{\ell }\right)\right)}\enspace \mathrm{sinh}\left(\sqrt{u/D}\left(x-{x}_{\ell }\right)\right),\quad \hspace{77pt}x{< }{x}_{0},\quad \hfill \\ \frac{1}{\sqrt{uD}}\frac{\mathrm{sinh}\left(\sqrt{u/D}\left({x}_{0}-{x}_{\ell }\right)\right)}{\enspace }\mathrm{cosh}\left(\sqrt{u/D}\left({x}_{u}-{x}_{\ell }\right)\right)\mathrm{cosh}\left(\sqrt{u/D}\left({x}_{u}-x\right)\right),\quad x{ >}{x}_{0},\quad \hfill \end{cases}\end{equation} \tag{ 27 }$$

solves equation (24). The fluxes onto the boundaries give us the FPT densities for a particle to reach the respective boundary. Thus with ${\tilde {\phi }}^{\xi }\left({x}_{\ell },{x}_{u},u;{x}_{0}\right)={\pm}D\frac{\partial }{\partial x}{\left.\tilde {\psi }\left(x,u\right)\right\vert }_{x={x}_{\xi }}$ (+ for the flux onto the lower and − for the flux onto the upper boundary) we have

$$\begin{equation}{\tilde {\phi }}_{\mathtt{\text{out}}}^{\ell }\left({x}_{\ell },{x}_{u},u;{x}_{0}\right)=\frac{\mathrm{cosh}\left(\sqrt{u/D}\left({x}_{u}-{x}_{0}\right)\right)}{\mathrm{cosh}\left(\sqrt{u/D}\left({x}_{u}-{x}_{\ell }\right)\right)},\end{equation} \tag{ 28 }$$

$$\begin{equation}{\tilde {\phi }}_{\mathtt{\text{out}}}^{u}\left({x}_{\ell },{x}_{u},u;{x}_{0}\right)=0.\end{equation} \tag{ 29 }$$

The cumulative FPT probability to reach x_ℓ for t → ∞ is obtained by putting u → 0:

$$\begin{equation}{\int }_{0}^{\infty }{\phi }_{\mathtt{\text{out}}}^{\ell }\left({x}_{\ell },{x}_{u},t;{x}_{0}\right)\mathrm{d}t=\underset{u\to 0}{\mathrm{lim}}\;{\tilde {\phi }}_{\mathtt{\text{out}}}^{\ell }\left({x}_{\ell },{x}_{u},u;{x}_{0}\right)=1,\end{equation} \tag{ 30 }$$

as expected. To see what happens in an infinite system, we let x_u → ∞:

$$\begin{equation}\underset{{x}_{u}\to \infty }{\mathrm{lim}}\enspace {\tilde {\phi }}_{\mathtt{\text{out}}}^{\ell }\left({x}_{\ell },{x}_{u},u;{x}_{0}\right)=\mathrm{exp}\left(-\sqrt{u/D}\left({x}_{0}-{x}_{\ell }\right)\right),\end{equation} \tag{ 31 }$$

which is the known density for first passage to a point x_ℓ on an infinite domain. An expansion in small u yields a non-analytic expression which gives rise to long time tail. (In fact the exponential of a square root of u gives the one sided Lévy stable law L_α of parameter α = 1/2 in t, thus the tail is ∝ t^−3/2.)

3.1.2. Inner region

Here we consider the boundary conditions

$$\begin{equation}\tilde {\psi }\left({x}_{\ell },u\right)=0,\end{equation} \tag{ 32 }$$

$$\begin{equation}\tilde {\psi }\left({x}_{u},u\right)=0,\end{equation} \tag{ 33 }$$

for which the solution to equation (24) yields

$$\begin{equation}\tilde {\psi }\left(x,u\right)=\begin{cases}\frac{1}{\sqrt{uD}}\frac{\mathrm{sinh}\left(\sqrt{u/D}\left({x}_{u}-{x}_{0}\right)\right)}{\enspace }\mathrm{sinh}\left(\sqrt{u/D}\left({x}_{u}-{x}_{\ell }\right)\right)\mathrm{sinh}\left(\sqrt{u/D}\left(x-{x}_{\ell }\right)\right),\quad x{< }{x}_{0},\quad \hfill \\ \frac{1}{\sqrt{uD}}\frac{\mathrm{sinh}\left(\sqrt{u/D}\left({x}_{0}-{x}_{\ell }\right)\right)}{\enspace }\mathrm{sinh}\left(\sqrt{u/D}\left({x}_{u}-{x}_{\ell }\right)\right)\mathrm{sinh}\left(\sqrt{u/D}\left({x}_{u}-x\right)\right),\quad x{ >}{x}_{0}.\quad \hfill \end{cases}\end{equation} \tag{ 34 }$$

Again, we compute the fluxes onto the boundaries [cf equation (23)]:

$$\begin{equation}{\tilde {\phi }}_{\mathtt{\text{in}}}^{\ell }\left({x}_{\ell },{x}_{u},u;{x}_{0}\right)=\frac{\mathrm{sinh}\left(\sqrt{u/D}\left({x}_{u}-{x}_{0}\right)\right)}{\mathrm{sinh}\left(\sqrt{u/D}\left({x}_{u}-{x}_{\ell }\right)\right)},\end{equation} \tag{ 35 }$$

$$\begin{equation}{\tilde {\phi }}_{\mathtt{\text{in}}}^{u}\left({x}_{\ell },{x}_{u},u;{x}_{0}\right)=\frac{\mathrm{sinh}\left(\sqrt{u/D}\left({x}_{0}-{x}_{\ell }\right)\right)}{\mathrm{sinh}\left(\sqrt{u/D}\left({x}_{u}-{x}_{\ell }\right)\right)}.\end{equation} \tag{ 36 }$$

These fluxes are the FPT densities to the respective boundaries. The splitting probabilities, i.e. the cumulative probabilities to leave the system via the respective boundary, are obtained by sending u → 0, which corresponds to t → ∞ in time domain:

$$\begin{equation}{P}_{{x}_{\ell }}\left({x}_{0}\right)={\int }_{0}^{\infty }{\phi }_{\mathtt{\text{in}}}^{\ell }\left({x}_{\ell },{x}_{u},t;{x}_{0}\right)\mathrm{d}t=\frac{{x}_{u}-{x}_{0}}{{x}_{u}-{x}_{\ell }},\end{equation} \tag{ 37 }$$

$$\begin{equation}{P}_{{x}_{u}}\left({x}_{0}\right)={\int }_{0}^{\infty }{\phi }_{\mathtt{\text{in}}}^{u}\left({x}_{\ell },{x}_{u},t;{x}_{0}\right)\mathrm{d}t=\frac{{x}_{0}-{x}_{\ell }}{{x}_{u}-{x}_{\ell }}.\end{equation} \tag{ 38 }$$

Sending the outer boundary to infinity,

$$\begin{equation}\underset{{x}_{u}\to \infty }{\mathrm{lim}}\enspace {\tilde {\phi }}_{\mathtt{\text{in}}}^{\ell }\left({x}_{\ell },{x}_{u},u;{x}_{0}\right)=\mathrm{exp}\left(-\sqrt{u/D}\left({x}_{0}-{x}_{\ell }\right)\right),\end{equation} \tag{ 39 }$$

we recover the known FPT density on an infinite domain.

3.1.3. Scaling form and full solution for the FPT density

Let now, closer in accordance with the original problem (figure 1(a)), start the partial trajectory at a small distance ɛ to its end point a, thus x₀ = a + ɛ, x_ℓ = a for ${\tilde {\phi }}_{\mathtt{\text{out}}}$ and x₀ = a, x_u = a + ɛ for ${\tilde {\phi }}_{\mathtt{\text{in}}}$. The lower boundary for the ${\phi }_{\mathtt{\text{in}}}^{u}$ be R and the upper boundary of the ${\phi }_{\mathtt{\text{out}}}^{\ell }$ be b. Furthermore, let us introduce the scaled Laplace variable s ≕ uɛ²/D. Then, we have

$$\begin{equation}{\tilde {\phi }}_{\mathtt{\text{in}}}^{u}\left(\varepsilon ,s\right)=\frac{\mathrm{sinh}\left(\sqrt{s}\left(R-a\right)/\varepsilon \right)}{\mathrm{sinh}\left(\sqrt{s}\left(\left(R-a\right)/\varepsilon -1\right)\right)},\end{equation} \tag{ 40 }$$

$$\begin{equation}{\tilde {\phi }}_{\mathtt{\text{out}}}^{\ell }\left(\varepsilon ,s\right)=\frac{\mathrm{cosh}\left(\sqrt{s}\left(\left(a-b\right)/\varepsilon +1\right)\right)}{\mathrm{sinh}\left(\sqrt{s}\left(a-b\right)/\varepsilon \right)}.\end{equation} \tag{ 41 }$$

The FPT densities of the partial trajectories attain their scaling forms for ɛ → 0:

$$\begin{equation}\underset{\varepsilon \to 0}{\mathrm{lim}}\enspace {\tilde {\phi }}_{\mathtt{\text{in}}}^{u}\left(\varepsilon ,s\right)=\mathrm{exp}\left(-\sqrt{s}\right),\end{equation} \tag{ 42 }$$

and the same for ${\mathrm{lim}}_{\varepsilon \to 0}\enspace {\tilde {\phi }}_{\mathtt{\text{out}}}^{\ell }\left(\varepsilon ,s\right)$, which permits an analytic inversion of the Laplace transform:

$$\begin{equation}\frac{1}{2\sqrt{\pi }{\left({t}^{{\ast}}\right)}^{3/2}}\enspace \mathrm{exp}\left(-\frac{1}{4{t}^{{\ast}}}\right)\end{equation} \tag{ 43 }$$

with t* := tD/ɛ² the scaled time variable conjugate to s.

The partial FPT densities for the inner and outer domains are depicted in the left panels of figure 3. The closer to the absorbing boundary or target a particle starts, the shorter the time of the peak position and the higher the peak value, indicating that the particles are more rapidly absorbed. Thus, moving the starting position closer to the target corresponds to a limiting procedure for the partial FPT densities, e.g., lim_ɛ→0  ϕ_in(ɛ, t) = δ₊(t), converging to a singular peak at t = 0: if the particle starts at the position of the target, all probability is absorbed immediately within zero time. The construction of an auxiliary ɛ-shell around the boundary at a, see section 2.5, circumvents the difficulties arising from such singular FPT densities. At intermediate times, we find the expected t^−3/2 power law decays for the particles exploring the outer domain yet without hitting the outer confinement. The FPT densities at large times decay rapidly due to the confinement.

Zoom In Zoom Out Reset image size

For the inner domain, there is another reflecting boundary at |x| = a. A sharp exponential cutoff sets in at times t ≈ a²/D (figures 3(a) and (b)), which we attribute to partial trajectories that end as soon as they reach the outer boundary and do not contribute to the FPT density at longer times. The outer domain has a reflecting boundary at |x| = b, which also results in an exponential cutoff of the tail (figures 3(c) and (d)). However, the trajectories are continued and move on in their attempts to reach the target, which they will eventually do, but at a later time. This shifts and smears out the cutoff at large times. Moreover, due to the conservation of probability, the reflected trajectories, which in the case of an unbounded domain would have contributed into the power law tail, are responsible for the small shoulder of the FPT density at large times. The right panels of the figures show the partial FPT densities scaled with respect to the distance of the starting point to the target ɛ. Small values of ɛ are equivalent to a large domain size b ≫ a, the particle feels the confinement at a later time and the differences between reflecting or absorbing outer boundary vanish.

With the partial FPT densities calculated above, we are ready to write down the analytical expression for the Laplace transform of the full FPT density for a particle starting at x₀ and being absorbed at R by substituting

$$\begin{align*}\hfill {\tilde {\phi }}_{0}\left(u\right)& ={\phi }_{\mathtt{\text{out}}}^{\ell }\left({x}_{\ell }=a,{x}_{u}=b,u;{x}_{0}\right),\hfill \\ \hfill {\tilde {\phi }}_{\mathtt{\text{in}}}^{\varnothing }\left(u\right)& ={\phi }_{\mathtt{\text{in}}}^{\ell }\left({x}_{\ell }=R,{x}_{u}=a+\varepsilon ,u;{x}_{0}=a\right),\hfill \\ \hfill {\tilde {\phi }}_{\mathtt{\text{in}}}\left(u\right)& ={\phi }_{\mathtt{\text{in}}}^{u}\left({x}_{\ell }=R,{x}_{u}=a+\varepsilon ,u;{x}_{0}=a\right)\enspace ,\quad \text{and}\hfill \\ \hfill {\tilde {\phi }}_{\mathtt{\text{out}}}\left(a,u\right)& ={\phi }_{\mathtt{\text{out}}}^{\ell }\left({x}_{\ell }=a,{x}_{u}=b,u;{x}_{0}=a+\varepsilon \right)\hfill \end{align*}$$

into equation (21).

As an ultimate test for the validity of the method we compare our GME solution for the FPT density for the homogeneous system (i.e. particles behave the same way in inner and outer region) with the FPT density for particles starting at r₀ and ending at first encounter with R (i.e. for non-dissected trajectories). We find perfect agreement between the two of them,

$$\begin{equation}{\tilde {p}}_{\text{FPT}}\left(u\right)={\tilde {\phi }}_{0}^{\ell }\left({x}_{\ell }=R,{x}_{u}=b,u;{x}_{0}\right)=\frac{\mathrm{cosh}\left(\sqrt{u/D}\left(b-{x}_{0}\right)\right)}{\mathrm{cosh}\left(\sqrt{u/D}\left(R-b\right)\right)}.\end{equation} \tag{ 44 }$$

Observing that in this case, the ɛ-dependence vanishes for the GME solution, taking the limit ɛ → 0 is not needed anymore. Also the dependence of a vanishes as it should in the completely homogeneous case. For an infinitely large outer domain, b → ∞, we have:

$$\begin{equation}{\tilde {p}}_{\text{FPT}}^{\infty }\left(u\right)=\mathrm{exp}\left(-\sqrt{u/D}\left({x}_{0}-R\right)\right).\end{equation} \tag{ 45 }$$

In time domain, this limiting form corresponds to

$$\begin{equation}{p}_{\text{FPT}}^{\infty }\left(t\right)=\frac{\left({x}_{0}-R\right)}{2\sqrt{\pi D}{t}^{3/2}}\enspace \mathrm{exp}\left(-\frac{{\left({x}_{0}-R\right)}^{2}}{4Dt}\right),\end{equation} \tag{ 46 }$$

which is the expected result for first passage on an infinite domain. As a last step, the inverse Laplace transform of expression (44) is calculated numerically and depicted in figure 4. It is indeed qualitatively the same picture as already for the partial FPT densities with reflecting outer boundary: for small times, it fits the Lévy–Smirnow law, equation (46), at large times, there is an exponential cutoff, and at FPTs smaller than the cutoff time, a small elevation relative to the t^−3/2 decay is visible, which collects the probability in the truncated tail.

Zoom In Zoom Out Reset image size

Analogously to our calculations in section 3.1, we now will calculate the partial FPT densities from the fluxes onto the boundaries of the respective PDEs for the radial symmetric setup (figure 1(b)). Thus in this case, equation (22) becomes in Laplace domain:

$$\begin{equation}\frac{1}{{r}^{2}}\frac{\partial }{\partial r}{r}^{2}D\frac{\partial }{\partial r}\tilde {\psi }\left(r,u\right)-u\tilde {\psi }\left(r,u\right)=-\frac{\delta \left(r-{r}_{0}\right)}{4\pi {r}_{0}}.\end{equation} \tag{ 47 }$$

By substituting $\tilde {\eta }\left(r,u\right){:=}\tilde {\psi }\left(r,u\right)/r$, the equation for $\tilde {\eta }$ attains a similar form as the one in the 1D case:

$$\begin{equation}D\frac{{\partial }^{2}}{\partial {r}^{2}}\tilde {\eta }\left(r\right)-u\tilde {\eta }\left(r\right)=-\frac{1}{4\pi {r}_{0}}\delta \left(r-{r}_{0}\right),\end{equation} \tag{ 48 }$$

with the initial condition specified on the right-hand side (all trajectories start initially from the radius r₀). The fundamental system is again $\left\{\mathrm{exp}\left(\sqrt{u/D}\enspace r\right),\mathrm{exp}\left(-\sqrt{u/D}\enspace r\right)\right\}\enspace .$

4.1.1. Outer region

For the outer region, we have the transformed boundary conditions (absorbing at r_ℓ , reflecting at r_u ):

$$\begin{equation}\tilde {\eta }\left({r}_{\ell },u\right)=0,\end{equation} \tag{ 49 }$$

$$\begin{equation}\frac{\partial }{\partial r}{\left.\tilde {\eta }\left(r,u\right)\right\vert }_{r={r}_{u}}=\frac{1}{{r}_{u}}\eta \left({r}_{u},u\right),\end{equation} \tag{ 50 }$$

so that the solution to equation (48) is given piecewise for r ≶ r₀ by

$$\begin{align}\hfill \tilde {\eta }\left(r{< }{r}_{0},u\right)=& \frac{1}{4\pi \sqrt{uD}{r}_{0}}\frac{{r}_{u}\sqrt{\frac{u}{D}}\enspace \mathrm{cosh}\left(\sqrt{u/D}\left({r}_{0}-{r}_{u}\right)\right)+\mathrm{sinh}\left(\sqrt{u/D}\left({r}_{0}-{r}_{u}\right)\right)}{{r}_{u}\sqrt{u/D}\enspace \mathrm{cosh}\left(\sqrt{u/D}\left({r}_{\ell }-{r}_{u}\right)\right)+\mathrm{sinh}\left(\sqrt{u/D}\left({r}_{\ell }-{r}_{u}\right)\right)}\hfill \\ \hfill & {\times}\mathrm{sinh}\left(\sqrt{u/D}\left(r-{r}_{\ell }\right)\right)\hfill \end{align} \tag{ 51 }$$

and

$$\begin{align}\hfill \tilde {\eta }\left(r{ >}{r}_{0},u\right)=& \;\frac{1}{4\pi \sqrt{uD}{r}_{0}}\frac{\mathrm{sinh}\left(\sqrt{u/D}\left({r}_{0}-{r}_{\ell }\right)\right)}{{r}_{u}\sqrt{u/D}\enspace \mathrm{cosh}\left(\sqrt{u/D}\left({r}_{\ell }-{r}_{u}\right)\right)+\mathrm{sinh}\left(\sqrt{u/D}\left({r}_{\ell }-{r}_{u}\right)\right)}\hfill \\ \hfill & {\times}\left[{r}_{u}\sqrt{u/D}\enspace \mathrm{cosh}\left(\sqrt{u/D}\left(r-{r}_{u}\right)\right)+\mathrm{sinh}\left(\sqrt{u/D}\left(r-{r}_{u}\right)\right)\right].\hfill \end{align} \tag{ 52 }$$

With ${\tilde {\psi }}^{\prime }={\tilde {\eta }}^{\prime }/r-\tilde {\eta }/{r}^{2}$ and ${\tilde {\phi }}^{\xi }\left({r}_{\ell },{r}_{u},u;{r}_{0}\right)={\pm}4D\pi {r}_{\xi }^{2}\frac{\partial }{\partial r}{\left.\tilde {\psi }\left(r,u\right)\right\vert }_{r={r}_{\xi }}$ (+ for the lower, − for the upper boundary) the fluxes onto the boundaries are thus

$$\begin{equation}{\tilde {\phi }}_{\mathtt{\text{out}}}^{\ell }\left({r}_{\ell },{r}_{u},u;{r}_{0}\right)=\frac{{r}_{\ell }}{{r}_{0}}\frac{{r}_{u}\sqrt{u/D}\enspace \mathrm{cosh}\left(\sqrt{u/D}\left({r}_{0}-{r}_{u}\right)\right)+\mathrm{sinh}\left(\sqrt{u/D}\left({r}_{0}-{r}_{u}\right)\right)}{{r}_{u}\sqrt{u/D}\enspace \mathrm{cosh}\left(\sqrt{u/D}\left({r}_{\ell }-{r}_{u}\right)\right)+\mathrm{sinh}\left(\sqrt{u/D}\left({r}_{\ell }-{r}_{u}\right)\right)},\end{equation} \tag{ 53 }$$

$$\begin{equation}{\tilde {\phi }}_{\mathtt{\text{out}}}^{u}\left({r}_{\ell },{r}_{u},u;{r}_{0}\right)=0.\end{equation} \tag{ 54 }$$

${\tilde {\phi }}_{\mathtt{\text{out}}}^{\ell }\left({r}_{\ell },{r}_{u},u;{r}_{0}\right)$ is the density for first passage from r₀ onto r_ℓ .

The integral over time is calculated as the limit u → 0, which yields

$$\begin{align}\hfill \underset{u\to 0}{\mathrm{lim}}\enspace {\tilde {\phi }}_{\mathtt{\text{out}}}^{\ell }\left({r}_{\ell },{r}_{u},u;{r}_{0}\right)& =\underset{v\to 0}{\mathrm{lim}}\left[\frac{{r}_{\ell }}{{r}_{0}}\frac{{r}_{u}v\enspace \mathrm{cosh}\left(v\left({r}_{0}-{r}_{u}\right)\right)+\mathrm{sinh}\left(v\left({r}_{0}-{r}_{u}\right)\right)}{\mathrm{cosh}\left(v\left({r}_{\ell }-{r}_{u}\right)\right)+\mathrm{sinh}\left(v\left({r}_{\ell }-{r}_{u}\right)\right)}\right]\hfill \\ \hfill & =\frac{{r}_{\ell }}{{r}_{0}}\frac{{r}_{0}}{{r}_{\ell }}=1\hfill \end{align} \tag{ 55 }$$

and confirms that with probability 1 the particles ultimately reach the boundary r_ℓ , as expected for a finite domain with the boundary at r_ℓ as the only exit. We may consider the transition to an infinite domain by taking the limit of large r_u :

$$\begin{equation}\underset{{r}_{u}\to \infty }{\mathrm{lim}}\enspace {\tilde {\phi }}_{\mathtt{\text{out}}}^{\ell }\left({r}_{\ell },{r}_{u},u;{r}_{0}\right)=\frac{{r}_{\ell }}{{r}_{0}}\enspace \mathrm{exp}\left(-\sqrt{u/D}\left({r}_{0}-{r}_{\ell }\right)\right),\end{equation} \tag{ 56 }$$

where again, as in the 1D case, the $\sqrt{u}$-term generates the long time tail. Note that the limit u → 0 in the infinite domain is less than unity, indicating the transience of diffusion in three dimensions.

4.1.2. Inner region

For the inner region, we have the boundary conditions:

$$\begin{equation}\tilde {\eta }\left({r}_{\ell },u\right)=0,\end{equation} \tag{ 57 }$$

$$\begin{equation}\tilde {\eta }\left({r}_{u},u\right)=0,\end{equation} \tag{ 58 }$$

which results in

$$\begin{equation}\tilde {\eta }\left(r,u\right)=\begin{cases}\frac{1}{4\pi \sqrt{uD}{r}_{0}}\frac{\mathrm{sinh}\left(\sqrt{u/D}\left({r}_{u}-{r}_{0}\right)\right)}{\mathrm{sinh}\left(\sqrt{u/D}\left({r}_{u}-{r}_{\ell }\right)\right)}\enspace \mathrm{sinh}\left(\sqrt{u/D}\left(r-{r}_{\ell }\right)\right),\quad \hfill & r{< }{r}_{0}\hfill \\ \frac{1}{4\pi \sqrt{uD}{r}_{0}}\frac{\mathrm{sinh}\left(\sqrt{u/D}\left({r}_{0}-{r}_{\ell }\right)\right)}{\mathrm{sinh}\left(\sqrt{u/D}\left({r}_{u}-{r}_{\ell }\right)\right)}\mathrm{sinh}\left(\sqrt{u/D}\left({r}_{u}-r\right)\right),\quad \hfill & r{ >}{r}_{0}.\hfill \end{cases}\end{equation} \tag{ 59 }$$

We find for the fluxes onto the boundaries:

$$\begin{equation}{\tilde {\phi }}_{\mathtt{\text{in}}}^{\ell }\left({r}_{\ell },{r}_{u},u;{r}_{0}\right)=\frac{{r}_{\ell }}{{r}_{0}}\frac{\mathrm{sinh}\left(\sqrt{u/D}\left({r}_{u}-{r}_{0}\right)\right)}{\mathrm{sinh}\left(\sqrt{u/D}\left({r}_{u}-{r}_{\ell }\right)\right)}\end{equation} \tag{ 60 }$$

$$\begin{equation}{\tilde {\phi }}_{\mathtt{\text{in}}}^{u}\left({r}_{\ell },{r}_{u},u;{r}_{0}\right)=\frac{{r}_{u}}{{r}_{0}}\frac{\mathrm{sinh}\left(\sqrt{u/D}\left({r}_{0}-{r}_{\ell }\right)\right)}{\mathrm{sinh}\left(\sqrt{u/D}\left({r}_{u}-{r}_{\ell }\right)\right)}.\end{equation} \tag{ 61 }$$

The time integrals, calculated as small-u limits, provide the splitting probabilities:

$$\begin{equation}\underset{u\to 0}{\mathrm{lim}}\enspace {\tilde {\phi }}_{\mathtt{\text{in}}}^{\ell }\left({r}_{\ell },{r}_{u},u;{r}_{0}\right)=\frac{{r}_{\ell }}{{r}_{0}}\frac{{r}_{u}-{r}_{0}}{{r}_{u}-{r}_{\ell }},\end{equation} \tag{ 62 }$$

$$\begin{equation}\underset{u\to 0}{\mathrm{lim}}\enspace {\tilde {\phi }}_{\mathtt{\text{in}}}^{u}\left({r}_{\ell },{r}_{u},u;{r}_{0}\right)=\frac{{r}_{u}}{{r}_{0}}\frac{{r}_{0}-{r}_{\ell }}{{r}_{u}-{r}_{\ell }},\end{equation} \tag{ 63 }$$

and indeed the sum of both is 1 as it should. Taking the limit r_u → ∞ yields the same cumulative FP probability for an infinite system as in the case of a reflecting upper boundary (equation (56)),

$$\begin{equation}\underset{{r}_{u}\to \infty }{\mathrm{lim}}\enspace {\tilde {\phi }}_{\mathtt{\text{in}}}^{\ell }\left({r}_{\ell },{r}_{u},u;{r}_{0}\right)=\frac{{r}_{\ell }}{{r}_{0}}\enspace \mathrm{exp}\left(-\sqrt{u/D}\enspace \left({r}_{0}-{r}_{\ell }\right)\right).\end{equation} \tag{ 64 }$$

4.1.3. Scaling form

Again, according to figure 1(b), we let the partial trajectory start at a small distance ɛ to its end point a, thus r₀ = a + ɛ, r_ℓ = a for ${\tilde {\phi }}_{\mathtt{\text{out}}}$ and r₀ = a, r_u = a + ɛ for ${\tilde {\phi }}_{\mathtt{\text{in}}}$. The inner radius for the ϕ_in be R and the outer radius of the ϕ_out be b. With the scaled Laplace variable s ≕ uɛ²/D, we have

$$\begin{equation}{\tilde {\phi }}_{\mathtt{\text{in}}}^{u}\left(\varepsilon ,s\right)=\frac{a+\varepsilon }{a}\frac{\mathrm{sinh}\left(\sqrt{s}\left(R-a\right)/\varepsilon \right)}{\mathrm{sinh}\left(\sqrt{s}\left[\left(R-a\right)/\varepsilon -1\right]\right)},\end{equation} \tag{ 65 }$$

$$\begin{equation}{\tilde {\phi }}_{\mathtt{\text{out}}}^{\ell }\left(\varepsilon ,s\right)=\frac{a}{a+\varepsilon }\frac{b\sqrt{s}/\varepsilon \enspace \mathrm{cosh}\left(\sqrt{s}\left[\left(a-b\right)/\varepsilon +1\right]\right)+\mathrm{sinh}\left(\sqrt{s}\left[\left(a-b\right)/\varepsilon +1\right]\right)}{b\sqrt{s}/\varepsilon \enspace \mathrm{cosh}\left(\sqrt{s}\left(a-b\right)/\varepsilon \right)+\mathrm{sinh}\left(\sqrt{s}\left(a-b\right)/\varepsilon \right)}.\end{equation} \tag{ 66 }$$

The limits as ɛ → 0 are equal to the scaling forms of the 1D case (equations (42) and (43)).

Figure 5 shows the partial FPT densities for the inner and outer domains (panels (a), (c)). Similarly to the 1D case we find again the peak at small times, which shifts to the left and gets the higher (and narrower) the closer to the target the particle starts. Intermediate times are again governed by a t^−3/2 power law decay. At large times, the cutoff sets in, again relatively sharply at t ≈ a²/D for the inner domain with its absorbing outer boundary at |x| = a, and more blurred out near t ≈ b²/D for the outer domain with its reflecting boundary at |x| = b. In the 3D case, the small shoulder at large times is more pronounced as compared to the 1D case. This is a well-known effect and is due to the compact exploration of space in low dimensions, where the first passage is more dominated by the direct trajectories [24]. A smaller large-time shoulder indicates that the geometry plays a minor role in low dimensions than in 3D or higher. Data collapse of the partial FPT densities is demonstrated again by rescaling with respect to the distance ɛ of the starting point to the target (figures 5(b) and (d)).

Zoom In Zoom Out Reset image size

Again we use the partial FPT densities to obtain the Laplace transform of the full FPT densities for a particle starting at r₀ and being absorbed at R by substituting

$$\begin{equation}{\tilde {\phi }}_{0}\left(u\right)={\phi }_{\mathtt{\text{out}}}^{\ell }\left({r}_{\ell }=a,{r}_{u}=b,u;{r}_{0}\right),\end{equation} \tag{ 67a }$$

$$\begin{equation}{\tilde {\phi }}_{\mathtt{\text{in}}}^{\varnothing }\left(u\right)={\phi }_{\mathtt{\text{in}}}^{\ell }\left({r}_{\ell }=R,{r}_{u}=a+\varepsilon ,u;{r}_{0}=a\right),\end{equation} \tag{ 67b }$$

$$\begin{equation}{\tilde {\phi }}_{\mathtt{\text{in}}}\left(u\right)={\phi }_{\mathtt{\text{in}}}^{u}\left({r}_{\ell }=R,{r}_{u}=a+\varepsilon ,u;{r}_{0}=a\right)\quad \text{and}\end{equation} \tag{ 67c }$$

$$\begin{equation}{\tilde {\phi }}_{\mathtt{\text{out}}}\left(a,u\right)={\phi }_{\mathtt{\text{out}}}^{\ell }\left({r}_{\ell }=a,{r}_{u}=b,u;{r}_{0}=a+\varepsilon \right)\end{equation} \tag{ 67d }$$

into equation (21).

The coincidence of our solution for the FPT density for the homogeneous system obtained via the GME approach with the time density for first passage to R starting from r₀, again underlines its vanishing dependence on ɛ and a,

$$\begin{align}\hfill {\tilde {p}}_{\text{FPT}}\left(u\right)& ={\tilde {\phi }}_{0}^{\ell }\left({r}_{\ell }=R,{r}_{u}=b,u;{r}_{0}\right)\hfill \\ \hfill & =\frac{R}{{r}_{0}}\frac{b\sqrt{u/D}\enspace \mathrm{cosh}\left(\sqrt{u/D}\left({r}_{0}-b\right)\right)+\mathrm{sinh}\left(\sqrt{u/D}\left({r}_{0}-b\right)\right)}{b\sqrt{u/D}\enspace \mathrm{cosh}\left(\sqrt{u/D}\left(R-b\right)\right)},\hfill \end{align} \tag{ 68 }$$

as expected. For an infinite outer domain, b → ∞, we have:

$$\begin{equation}{\tilde {p}}_{\text{FPT}}^{\infty }\left(u\right)={\tilde {\phi }}_{0}^{\infty }\left({r}_{\ell }=R,u;{r}_{0}\right)=\frac{R}{{r}_{0}}\enspace \mathrm{exp}\left(-\sqrt{u/D}\left({r}_{0}-R\right)\right),\end{equation} \tag{ 69 }$$

which is translated to the time domain as

$$\begin{equation}{p}_{\text{FPT}}^{\infty }\left(t\right)=\frac{R}{{r}_{0}}\frac{\left({r}_{0}-R\right)}{2\sqrt{\pi D}{t}^{3/2}}\enspace \mathrm{exp}\left(-\frac{{\left({r}_{0}-R\right)}^{2}}{4Dt}\right).\end{equation} \tag{ 70 }$$

The numerical Laplace inversion of equation (68) yields the final FTP densities p_FPT(t), shown in figure 6 and deviating qualitatively from the 1D case (figure 4): the plateau region is much more pronounced in the 3D case. Notice that the limiting FPT density for infinite domains (equation (70)) does not normalize due to its prefactor, while the FPT densities in confined domains do. The plateau in the FPT profile for reflecting boundaries does not only compensate for those particles that would have contributed to the power law tail in the case of an unbounded domain, but in addition it compensates also for the particles that would have got lost forever due to the transient character of diffusion in dimensions higher than two.

Zoom In Zoom Out Reset image size