A First Look at First-Passage Processes

The first-passage probability is defined as the probability that a diffusing particle or a random walk first reaches a given site (or set of sites) at a specified time. Typical examples of first-passage processes include: fluorescence quenching, in which light emission by a fluorescent molecule stops when it reacts with a quencher; integrate-and-fire neurons, in which a neuron fires when a fluctuating voltage level first reaches a specified level; and the execution of buy/sell orders when a stock price first reaches a threshold. To appreciate why first-passage phenomena might be relevant practically, consider the following example. You are an investor who buys stock in a company at a price of $100. Suppose that this price fluctuates daily by $\pm\$1$. You will sell if the stock price reaches $200 and if the stock price reaches $0, the company has gone bankrupt and you’ve lost all your investment. What it the probability of doubling your investment or losing your entire investment? How long will it take before one of these two events occurs? These are the types of questions that are the purview of first-passage phenomena. Much of the material covered here is discussed more detail in this monograph [1], and in other general reviews and texts on probability theory and stochastic processes [2, 3, 4, 5].

Let’s start by deriving the formal relation between the first-passage probability and the familiar occupation probability. For concreteness, consider a random walk in discrete space and discrete time. We define $P(\mathbf{r},t)$ as the occupation probability; this is the probability that a random walk is at site $\mathbf{r}$ at time $t$ when it starts at the origin. Similarly, let $F(\mathbf{r},t)$ be the first-passage probability, namely, the probability that the random walk first visits $\mathbf{r}$ at time $t$ with the same initial condition. Clearly $F(\mathbf{r},t)$ decays more rapidly in time than $P(\mathbf{r},t)$ because once a random walk reaches $\mathbf{r}$, there can be no further contribution to $F(\mathbf{r},t)$, although the same walk may still contribute to $P(\mathbf{r},t)$.

It is convenient to write $P(\mathbf{r},t)$ in terms of $F(\mathbf{r},t)$ and then invert this relation to find $F(\mathbf{r},t)$. For a random walk to be located at $\mathbf{r}$ at time $t$, the walk must first reach $\mathbf{r}$ at some earlier time step $t^{\prime}$ and then return to $\mathbf{r}$ after $t-t^{\prime}$ additional steps (Fig. 1). This connection between $F(\mathbf{r},t)$ and $P(\mathbf{r},t)$ may be expressed as the convolution

$$P(\mathbf{r},t)=\delta_{\mathbf{r},0}\delta_{t,0}+\sum_{t^{\prime}\leq t}F(\mathbf{r},t^{\prime})\,P(0,t-t^{\prime})\,,$$

(2.1)

where $\delta_{t,0}$ is the Kronecker delta function. This equation expresses the fact that if a random walk is at $\mathbf{r}$ at time $t$, it must have first reached $\mathbf{r}$ at some earlier time $t^{\prime}$ (which could even be $t$). If the walk reached $\mathbf{r}$ at a time earlier than $t$, then it must return to $\mathbf{r}$ (and any number of such returns could occur) in the remaining time $t-t^{\prime}$. The probability for this set of events is expressed by $P(0,t-t^{\prime})$. The delta function term accounts for the initial condition that the walk starts at $\mathbf{r}=0$.

The above convolution is most conveniently solved by introducing the generating functions,

$\displaystyle P(\mathbf{r},z)=\sum_{t=0}^{\infty}P(\mathbf{r},t)z^{t}\,,\qquad F(\mathbf{r},z)=\sum_{t=0}^{\infty}F(\mathbf{r},t)z^{t}\,.$

For a random walk in continuous time, we would merely replace the sum over discrete time in Eq. (2.1) by an integral and then use the Laplace transform. However, the asymptotic results would be identical. To solve for the first-passage probability, we multiply Eq. (2.1) by $z^{t}$ and sum over all $t$. We thereby find that the generating functions for $P$ and $F$ are related by

$$P(\mathbf{r},z)=\delta_{\mathbf{r},0}+F(\mathbf{r},z)P(0,z)\,.$$

(2.2)

Thus we obtain the fundamental connection between the generating functions

$$F(\mathbf{r},z)=\begin{cases}{\displaystyle\frac{P(\mathbf{r},z)}{P(0,z)}}&\mathbf{r}\neq 0\,,\\[11.38109pt] {\displaystyle 1-\frac{1}{P(0,z)}},&\mathbf{r}=0\,.\end{cases}$$

(2.3)

The important point is that the first-passage probability can be determined solely from the occupation probability. Many profound results about random walks in infinite space can be obtained from the fundamental relation (2.3) (see, e.g., [6, 7, 8, 9]). Our focus here will be on random walks or diffusion in confined geometries that reflect important physical constraints.

Suppose that a diffusing particle starts at $x_{0}>0$ on the infinite half line $[0,\infty]$ and is absorbed when it reaches the origin. Does the particle ever reach the origin? If so, when does this particle first reach the origin? To answer these questions, we have to solve the diffusion equation for the concentration $c(x,t)$, subject to the initial condition $c(x,t\!=\!0)=\delta(x-x_{0})$, and the boundary condition $c(x\!=\!0,t\!>\!0)=0$; the latter enforces the absorption of the particle when it reaches the origin.

A standard approach to solve this problem is to first take the Laplace transform of the diffusion equation and then solve for the Green’s function of this transformed equation. Then one inverts the Laplace transform of the Green’s function to obtain the concentration $c(x,t)$ in the time domain. The diffusive flux to the origin gives the probability that the particle reaches the origin at time $t$. Because of the absorbing boundary condition, when the particle does reach the origin, it is removed from the system. Thus the diffusive flux corresponds to the probability for the particle to reach the origin for the first time—namely, the first-passage probability to the origin.

A more fun way to solve this problem is by invoking the familiar image method from electrostatics. In this method, a diffusing particle that starts at $x_{0}>0$ and is subject to an absorbing boundary condition at $x=0$ is equivalent to removing the boundary altogether and introducing an image “antiparticle” that is initially at $x=-x_{0}$. Because the concentration is clearly equal to zero at the origin by symmetry, this image antiparticle effectively imposes the absorbing boundary condition $c(x\!=\!0,t)=0$. The initial particle and the image antiparticle both diffuse freely on $[-\infty,\infty]$ and their superposition gives a resultant concentration $c(x,t)$ for $x>0$ that solves the original problem.

Hence the concentration for a diffusing particle on the positive half-line is the sum of a Gaussian centered at $x_{0}$ and an anti-Gaussian centered at $-x_{0}$:

$$c(x,t)=\frac{1}{\sqrt{4\pi Dt}}\left[e^{-(x-x_{0})^{2}/4Dt}-e^{-(x+x_{0})^{2}/4Dt}\right].$$

(3.1)

This concentration profile has a linear dependence on $x$ near the origin and a Gaussian tail for $x/\sqrt{Dt}\gg 1$, as illustrated in Fig. 2. Because the initial condition is normalized, the first-passage probability to the origin at time $t$ is just the diffusive flux to this point. From the above expression for $c(x,t)$, we find

$\displaystyle F(0,t)=+D\frac{\partial c(x,t)}{\partial x}\bigg{|}_{x=0}=\frac{x_{0}}{\sqrt{4\pi Dt^{3}}}\,\,e^{-x_{0}^{2}/4Dt}\qquad t\to\infty\,.$

(3.2)

This fundamental and simple formula has a number of striking implications:

1. 1.

   The particle is sure to reach to the origin because $\int_{0}^{\infty}F(0,t)\,dt=1$. That is, eventual absorption is certain.

2. 2.

   The average time for the particle to reach the origin is infinite! This fact arises because the first-passage probability has the long-time algebraic tail $F(0,t)\to x_{0}/\sqrt{4\pi Dt^{3}}$ for $t\gg x_{0}^{2}/D$. This dichotomy between hitting the origin with certainty but taking an infinite average time to do so underlies many of the intriguing features of one-dimensional diffusion.

3. 3.

   The typical time to reach the origin is finite. We can define the term typical time in a precise way as follows: As a preliminary, define the typical position of the particle, $x_{T}$, as

   $\displaystyle\int_{x_{T}}^{\infty}c(x,t)\,dx=\tfrac{1}{2}\,.$

   (3.3)

   That is, one half of the total probability lies in the range beyond $x_{T}$ and one half lies in the range $[0,x_{T}]$. Substituting the concentration profile (3.1) into the above integral and performing the integral leads to the transcendental equation

   $\displaystyle\text{erfc}\left(\frac{x_{T}-x_{0}}{\sqrt{4Dt}}\right)-\text{erfc}\left(\frac{x_{T}+x_{0}}{\sqrt{4Dt}}\right)=1\,.$

   (3.4)

   This equation can only be solved numerically and the salient feature is that $x_{T}$ monotonically decreases with time and reaches zero at a time that is roughly $1.1\,x_{0}^{2}/D$; this defines the typical hitting time.

4. 4.

   Even though the average time to reach the origin is infinite, the number of times the origin is reached in a time $t$ is proportional to $\sqrt{t}$. This result follows directly from the probability distribution of a freely diffusing particle. At large times, the bulk of Gaussian distribution for a particle that starts at $x_{0}$ will spread past the origin. Each site that is within the Gaussian envelope will have been visited of the order of $\sqrt{t}$ times. This last fact will be relevant for our discussion of the reaction rate of the sphere in Sec. 6.

We now turn to the first-passage properties in a finite interval. The reason for focusing on the finite interval is that the basic first-passage questions in this geometry have many profound implications. Moreover, the interval geometry is sufficiently simple that many results can be readily derived. Let us begin by outlining the basic questions that we will address. We consider a diffusing particle that starts at some point $x_{0}$ within the interval $[0,L]$, with absorbing boundary conditions at both ends of the interval. Eventually the particle is absorbed, and our goal is to characterize the time dependence of this absorption. Basic first-passage questions include:

1. 1.

   What is the time dependence of the survival probability $S(t)$? This is the probability that a diffusing particle does not touch either absorbing boundary before time $t$.

2. 2.

   What is the time dependence of the first-passage, or exit probabilities, to either 0 or to $L$ as a function of $x_{0}$? Integrating these probabilities over all time gives the eventual hitting, or splitting probability to a specified boundary. What is the dependence of the splitting probability to 0 or to $L$ as a function of the starting position?

3. 3.

   What is the average exit time, that is, the average time until the particle hits either of the absorbing boundaries as a function of starting position? What are the conditional exit times, that is, the average time to hit a specified boundary (without ever touching the other boundary) as a function of the starting position?

To answer the first question, we need to solve the diffusion equation in the interval, subject to the initial condition that a particle starts at $x_{0}$, and with absorbing boundary conditions at both ends. This is a standard exercise and the result for the concentration is

$$c(x,t)=\sum_{n\geq 1}^{\infty}\frac{2}{L}\sin\frac{n\pi x_{0}}{L}\,\sin\frac{n\pi x}{L}\;e^{-n^{2}\pi^{2}\,Dt/L^{2}}~{}.$$

(4.1)

Since the large-$n$ eigenmodes decay more rapidly in time, the most slowly decaying eigenmode dominates in the long-time limit. As a result, the survival probability, which is the spatial integral of the concentration over the interval, asymptotically decays as

$$S(t)\sim e^{-D\pi^{2}t/L^{2}}\,.$$

(4.2)

For answering the second question, it is mathematically simpler to work in the Laplace transform domain. Applying the Laplace transform to the diffusion equation recasts it as

$$sc(x,s)-c(x,t=0)=Dc^{\prime\prime}(x,s)\,,$$

(4.3)

where the prime denotes differentiation with respect to $x$. The argument $s$ indicates that $c(x,s)$ is the Laplace transform. Within the standard Green’s function approach, the homogeneous equation in each subdomain $(0,x_{0})$ and $(x_{0},L)$ has elementary solutions of the form $c(x,s)=A\exp(x\sqrt{{s/D}})+B\exp(-x\sqrt{{s/D}})$, with the constants $A$ and $B$ determined by the boundary conditions. Because the absorbing boundary condition at $x=0$ mandates an antisymmetric combination of exponentials and because the form of the Green’s function as $x\to L$ must be identical to that as $x\to 0$, we can be immediately write

$\displaystyle\begin{split}c_{<}(x,s)&=A\sinh\left(\sqrt{\frac{s}{D}}\,x\right)\hskip 36.135ptx<x_{0}\,,\\[5.69054pt] c_{>}(x,s)&=B\sinh\left(\sqrt{\frac{s}{D}}(L-x)\right)\quad x>x_{0}\,,\end{split}$

(4.4)

for the subdomain Green’s functions $c_{<}$ and $c_{>}$, where $A$ and $B$ are constants.

We now impose the continuity condition $c_{<}(x_{0},s)=c_{>}(x_{0},s)$ and the jump condition that is obtained by integrating Eq. (4.3) over an infinitesimal interval that includes $x_{0}$:

$$c^{\prime}(x,s)\big{|}_{x=x_{0}^{+}}-c^{\prime}(x,s)\big{|}_{x=x_{0}^{-}}=-1/D\,,$$

to finally obtain

$$c(x,s)=\frac{\displaystyle\sinh\left(\sqrt{\frac{s}{D}}\,x_{<}\right)\sinh\left(\sqrt{\frac{s}{D}}(L-x_{>})\right)}{\displaystyle\sqrt{sD}\sinh\left(\sqrt{\frac{s}{D}}\,L\right)}\,.$$

(4.5)

From this Green’s function, the Laplace transform of the fluxes at $x=0$ and at $x=L$ are

	$\displaystyle j_{-}(s|x_{0})$	$\displaystyle\equiv+D\frac{\partial c(x,s)}{\partial x}\Bigg{|}_{x=0}={\displaystyle\sinh\left({\sqrt{\frac{s}{D}}}(L-x_{0})\right)}\bigg{/}{\displaystyle\sinh\left({\sqrt{\frac{s}{D}}}\,L\right)}\,.$		(4.6a)
	$\displaystyle j_{+}(s|x_{0})$	$\displaystyle\equiv-D\frac{\partial c(x,s)}{\partial x}\Bigg{|}_{x=L}={\displaystyle\sinh\left(\sqrt{\frac{s}{D}}\,x_{0}\right)}\bigg{/}{\displaystyle\sinh\left(\sqrt{\frac{s}{D}}\,L\right)}\,.$		(4.6b)

The subsidiary argument $x_{0}$ in $j$ emphasizes that the flux depends on the initial particle position. Since the initial condition is normalized, the magnitude of the flux to each boundary is identical to the respective first-passage probabilities.

For $s=0$, these Laplace transforms are just the time-integrated first-passage probabilities to 0 and at $L$. These quantities therefore coincide with the respective splitting probabilities, ${\mathcal{E}}_{-}(x_{0})$ and ${\mathcal{E}}_{+}(x_{0})$, namely, the probabilities to eventually hit the left and the right ends of the interval as a function of the initial position $x_{0}$:

$\displaystyle\begin{split}{\mathcal{E}}_{-}(x_{0})&=j_{-}(s\!=\!0|x_{0})=1-\frac{x_{0}}{L}\,,\\[5.69054pt] {\mathcal{E}}_{+}(x_{0})&=j_{+}(s\!=\!0|x_{0})=\frac{x_{0}}{L}\,.\end{split}$

(4.7)

Thus the splitting probabilities are given by an amazingly simple formula—the probability of reaching one endpoint is just the fractional distance to the other endpoint!

It is instructive to also derive these splitting probabilities by the backward Kolmogorov approach [1, 2]. The word backward reflects the feature that the initial condition becomes the dependent variable, rather than the current position of the particle. As we shall see, this method provides a powerful tool for determining first-passage properties. Physically, we obtain the eventual hitting probability ${\mathcal{E}}_{+}(x_{0})$ to the right boundary by summing the probabilities for all paths that start at $x_{0}$ and reach $L$ without touching $0$. Thus

$${\mathcal{E}}_{+}(x_{0})=\sum_{\text{paths}}\Pi_{x_{0}\to L}\,,$$

(4.8)

where $\Pi_{x_{0}\to L}$ denotes the probability of a path from $x_{0}$ to $L$ that avoids $0$. As illustrated in Fig. 3, the sum over all such paths can be decomposed into the outcome after one step and the sum over all path remainders from the intermediate point $x^{\prime}$ to $L$. This gives

$\displaystyle\mathcal{E}_{+}(x_{0})=\tfrac{1}{2}\sum_{\text{paths}^{\prime}}\Pi_{x_{0}+\delta x\to L}+\tfrac{1}{2}\,\sum_{\text{paths}^{\prime\prime}}\Pi_{x_{0}-\delta x\to L}=\tfrac{1}{2}[{\mathcal{E}}_{+}(x_{0}\!+\!\delta x)+{\mathcal{E}}_{+}(x_{0}\!-\!\delta x)]\,.$

(4.9)

Here $\delta x$ is the length of a single random-walk step and paths^′ and paths^′′ indicate, respectively, all paths that start at $x_{0}\pm\delta x$ and reach $L$ without touching 0.

Equation (4.9) reduces to $\Delta^{(2)}{\mathcal{E}}_{\pm}(x_{0})=0$, where $\Delta^{(2)}$ is the discrete second-difference operator, $\Delta^{(2)}f(x)\equiv f(x-\delta x)-2f(x)+f(x+\delta x)$. This difference equation is subject to the boundary conditions ${\mathcal{E}}_{+}(0)=0$, ${\mathcal{E}}_{+}(L)=1$. The solution is simply

$\displaystyle{\mathcal{E}}_{+}(x_{0})=\frac{x}{L}\,,$

(4.10)

and correspondingly ${\mathcal{E}}_{-}(x)=1-\frac{x}{L}$.

It is worth mentioning that there is an even simpler way to determine the exit probabilities by the martingale method [10]. This solution relies on the fact that the motion of the random walk in the interval is a “fair game” at any time. Thus the average position of the walk is time independent. Formally, a martingale is a process in which the average value of a random variable at time $t+\delta t$ equals the average value of this variable at time $t$.

For the present example, at $t=0$, the average position of the particle is $\langle x\rangle=x_{0}$. At infinite time, the walk is either at the left end or the right end of the interval, with respective probabilities ${\mathcal{E}}_{+}(x_{0})$ or ${\mathcal{E}}_{-}(x_{0})$. Thus the average position at infinite time is $\langle x\rangle=0\times{\mathcal{E}}_{-}(x_{0})+L\times{\mathcal{E}}_{+}(x_{0})$. Since the initial average position equals the final average position, we immediately recover (4.10).

For a biased random walk with a probability $p$ of hopping to the right and probability $q$ of hopping to the left, the analog of Eq. (4.9) for the splitting probability is

$\displaystyle{\mathcal{E}}(x_{0})=p{\mathcal{E}}_{+}(x_{0}\!+\!\delta x)+q{\mathcal{E}}_{+}(x_{0}\!-\!\delta x)\,,$

(4.11)

with solution

$\displaystyle{\mathcal{E}}_{+}(x_{0})=\frac{1-e^{-vx_{0}/D}}{1-e^{-vL/D}}\equiv\frac{1-e^{-u_{0}P\!e}}{1-e^{-P\!e}}\,,$

(4.12)

where $v=(p-q)\delta x$, $D=\delta x^{2}/2$, $u_{0}=x_{0}/L$, and $P\!e=vL/D$ is the Péclet number, which is a dimensionless measure of the influence of the bias relative to diffusive fluctuations. When the Péclet number is large, the exit probability clearly reflects the strong influence of the bias (Fig. 4).

We now extend the backward Kolmogorov approach to determine the average exit time from the finite interval. We distinguish between the unconditional average exit time, namely, the average time for a particle to reach either end of the interval, and the conditional average exit time, namely, the average time for a particle to reach, say, the right end of the interval without ever touching the other end.

In close analogy with (4.8), the unconditional exit time satisfies

$$t(x_{0})=\sum_{\text{paths}}(\Pi t)_{x_{0}\to\pm L}\,.$$

(4.13)

That is, to compute the average unconditional exit time, we take the time for a path to go from $x_{0}$ to $\pm L$ times the probability of this path and sum over all possible paths. Using the same decomposition that which led to (4.9), the unconditional exit time satisfies

$\displaystyle t(x_{0})$

$\displaystyle=\tfrac{1}{2}\!\!\sum_{\text{paths}^{\prime}}\!\!\Pi_{x_{0}+\delta x\to\pm L}(t_{x_{0}+dx\to\pm L}\!+\!dt)+\tfrac{1}{2}\!\!\sum_{\text{paths}^{\prime\prime}}\!\!\Pi_{x_{0}-\delta x\to\pm L}(t_{x_{0}-dx\to\pm L}\!+\!dt)\,.$

(4.14)

Notice that the term

$\displaystyle\tfrac{1}{2}\!\!\sum_{\text{paths}^{\prime}}\!\!\Pi_{x_{0}+\delta x\to\pm L}\;t_{x_{0}+dx\to\pm L}$

has exactly the same form as (4.13), so that the above expression is merely $\frac{1}{2}t(x_{0}+\delta x)$. A similar identification holds for the analogous term

$\displaystyle\tfrac{1}{2}\!\!\sum_{\text{paths}^{\prime\prime}}\!\!\Pi_{x_{0}-\delta x\to\pm L}\;t_{x_{0}-dx\to\pm L}\,.$

Finally, the terms multiplying the factor $dt$ in (4.14) just gives the probability of all possible paths from $x_{0}$ to either end of the interval; this probability is clearly equal to 1. Thus we have

$\displaystyle t(x_{0})=\tfrac{1}{2}[t(x_{0}\!+\!\delta x)+t(x_{0}\!-\!\delta x)]+dt\,.$

(4.15)

In the continuum limit, we expand (4.15) in a Taylor series to second order in $\delta x$ to give

$\displaystyle\frac{\delta x^{2}}{2}t^{\prime\prime}=-dt\,.$

Now identifying $\delta x^{2}/(2dt)$ as the diffusion coefficient $D$, the equation for the unconditional exit time reduces to $D\,t^{\prime\prime}=-1$, subject to the boundary conditions $t(0)=t(L)=0$. The solution is (see Fig. 5)

$$t(x_{0})=\frac{x_{0}(L-x_{0})}{2D}\,.$$

(4.16)

Notice that this exit time is of the order of $L$ for a particle that starts a distance of the order of one from an absorbing boundary and of the order of $L^{2}$ for a particle that starts near the middle of the interval.

Now let’s calculate the conditional exit time to the right boundary when starting from $x_{0}$, $t_{+}(x_{0})$. By definition this conditional time is

$\displaystyle t_{+}(x_{0})=\frac{\sum_{\text{paths}_{+}}(\Pi t)_{x_{0}\to L}}{\sum_{\text{paths}_{+}}\Pi_{x_{0}\to L}}=\frac{\sum_{\text{paths}_{+}}(\Pi t)_{x_{0}\to L}}{{\mathcal{E}}_{+}(x_{0})}\,.$

(4.17)

That is, the conditional exit time is the time for a path to start at $x_{0}$ and reach $L$ without touching 0 multiplies by the probability of this path, summed over all such allowable paths. Here the subscript + on the word paths indicates that only paths that go from $x_{0}$ to $L$ without touching 0 are included. Since the total probability of all these restricted paths is less than 1, we need to divide by this total probability, ${\mathcal{E}}_{+}(x_{0})$, to obtain the properly normalized conditional exit time. Thus, for the quantity ${\mathcal{E}}_{+}(x_{0})t_{+}(x_{0})$, we have

	$\displaystyle{\mathcal{E}}_{+}(x_{0})t_{+}(x_{0})$	$\displaystyle=\sum_{\text{paths}_{+}}(\Pi t)_{x_{0}\to L}$
		$\displaystyle=\tfrac{1}{2}\!\!\sum_{\text{paths}^{\prime}_{+}}\Pi(t_{x_{0}+dx\to L}+dt)+\tfrac{1}{2}\!\!\sum_{\text{paths}^{\prime\prime}_{+}}\Pi(t_{x_{0}-dx\to L}+dt)$
		$\displaystyle=\tfrac{1}{2}\left[{\mathcal{E}}_{+}(x_{0}\!+\!dx)t_{+}(x_{0}\!+\!dx)\right]+\tfrac{1}{2}\left[{\mathcal{E}}_{+}(x_{0}\!-\!dx)t_{+}(x_{0}\!-\!dx)\right]$
		$\displaystyle\hskip 158.99377pt+{\mathcal{E}}_{+}(x_{0})\,dt\,.$		(4.18)

In the continuum limit, the above equation reduces to $D({\mathcal{E}}_{+}t_{+})^{\prime\prime}=-{\mathcal{E}}_{+}$, with solution

$\displaystyle t_{+}(x_{0})=\frac{L^{2}-x_{0}^{2}}{6D}\,.$

(4.19)

From this result, we also immediately find $t_{-}(x_{0})=t_{+}(L-x_{0})$. The dependences of all the exit times on the starting position are illustrated in Fig. 5.

One of the alluring features of first-passage processes is its intimate connection to electrostatics. By this connection, one can recast an electrostatic problem in a given geometry as a first-passage problem in the same geometry. With this perspective, it is possible to solve seemingly difficult first-passage problems in a simple way by this electrostatic connection.

To illustrate the basic principle, consider the following general problem. Suppose that a diffusing particle starts at some point $\mathbf{r}_{0}$ inside an arbitrary bounded domain. At the boundary of this domain, the particle is absorbed. Eventually, all of the initial probability is absorbed on the boundary and we ask: what is the exit probability at some arbitrary point $\mathbf{r}_{B}$ on the domain boundary? Formally, we have to solve the diffusion equation in this domain, subject to the appropriate initial and boundary conditions:

$\displaystyle\begin{cases}\displaystyle{\frac{\partial c(\mathbf{r},t)}{\partial t}}=D\nabla^{2}c(\mathbf{r},t)\\[5.69054pt] c(\mathbf{r},t=0)=\delta(\mathbf{r}-\mathbf{r}_{0})\\[2.84526pt] c(\mathbf{r},t)\big{|}_{\mathbf{r}_{B}}=0\end{cases}$

(5.1)

Then the exit probability at the boundary point $\mathbf{r}_{B}$ is given by

$\displaystyle\mathcal{E}(\mathbf{r}_{B})=-D\int_{0}^{\infty}\frac{\partial c(\mathbf{r},t)}{\partial\hat{\mathbf{n}}}\bigg{|}_{\mathbf{r}_{B}}\,dt\,,$

(5.2)

where $\hat{\mathbf{n}}$ is the outward normal to the surface of the domain at $\mathbf{r}_{B}$.

Let’s look critically at this calculation. We are attempting to solve a partial differential equation in some domain (which may well be difficult), then take the result of this calculation and integrate over all time. That is, we really don’t need that exit probability at all times, but merely the time integral of the exit probability. This observation suggests that it will be useful to take the original problem (5.1) and integrate it over all time. To simplify what emerges, we also define the time integrated concentration, $\mathcal{C}(\mathbf{r})\equiv\int_{0}^{\infty}c(\mathbf{r},t)\,dt$. Performing this time integration on Eq. (5.1) leads to

$\displaystyle\begin{cases}-\delta(\mathbf{r}-\mathbf{r}_{0})=D\nabla^{2}\mathcal{C}(\mathbf{r})\\ \mathcal{C}(\mathbf{r})\big{|}_{\mathbf{r}_{B}}=0\end{cases}$

The delta function on the left-hand side is what remains when we integrate the time derivative in (5.1) over all time. At $t=\infty$ the concentration is zero, while at $t=0$, we merely have the initial condition. But notice that in terms of the time-integrated concentration, the exit probability may be written as

$\displaystyle\mathcal{E}(\mathbf{r}_{B})=-D\frac{\partial\mathcal{C}(\mathbf{r},t)}{\partial\hat{\mathbf{n}}}\bigg{|}_{\mathbf{r}_{B}}\,.$

(5.3)

Thus we arrive at the fundamental result:

	$\displaystyle\mathcal{E}(\mathbf{r}_{B})$	$\displaystyle=\text{the\ electric\ field\ at\ }\mathbf{r}_{B}\text{\ on the surface of a grounded conductor,}$
		$\displaystyle\hskip 14.45377pt\text{when a point charge of magnitude\ }1/(D\Omega_{d})\text{\ is placed at\ }\mathbf{r}_{0}\,.$

Here $\Omega_{d}$ is the surface area of a $d$-dimensional sphere; this factor is needed to convert the prefactor $D$ in the diffusion equation to the correct prefactor in the Laplace equation. Thus a given first-passage problem can be expressed as an equivalent electrostatic problem in the same geometry.

What is the probability $H(r)$ that a diffusing particle eventually hits a sphere of radius $a$, when the particle starts at a distance $r>a$ from the origin? One can determine this hitting probability in the standard way by solving the diffusion equation exterior to the sphere, computing the flux to the sphere, and then integrating over all time. However, it is much simpler to use the connection between first passage and electrostatics. Indeed, by a direct extension of Eq. (4.9) to three dimensions, the hitting probability satisfies (here for a discrete random walk in Cartesian coordinates for simplicity)

	$\displaystyle H(x,y)$	$\displaystyle=\tfrac{1}{6}\big{[}H(x+1,y,z)+H(x-1,y,z)+H(x,y+1,z)$
		$\displaystyle\hskip 36.135pt+H(x,y-1,z)+H(x,y,z+1)+H(x,y,z-1)\big{]}\,.$		(6.1)

Let us now take the continuum limit so that we can work in spherical coordinates. Then the above discrete difference equation becomes $\nabla^{2}H=0$, subject to the boundary conditions $H(r\!=\!a)=1$ and $H(\infty)=0$. The solution is

$\displaystyle H(r)=\frac{a}{r}\,.$

(6.2)

Amazingly simple!

Now let’s treat a related problem that is fundamental in chemical kinetics. Suppose that there is an initially uniform concentration of particles exterior to an absorbing sphere of radius $a$. A fundamental kinetic characteristic of the absorbing sphere is its reaction rate $k$, namely, the efficiency at which this sphere captures particles. Formally, $k$ is defined as the number of particles absorbed per unit time divided by the initial concentration. This normalization ensures that the reaction rate is a quantity that is intrinsic to the system. By dimensional analysis, the reaction rate as defined above has units of (length)${}^{d}/$time. Moreover the reaction rate can only be a function of intrinsic parameters of the system, namely, the diffusion coefficient $D$ and the sphere radius $a$. Since the units of $D$ are length²/time, we infer that $k$ must have units of $D\,a^{d-2}$. Thus the reaction rate $k=A\,D\,a^{d-2}$, where $A$ is a constant of the order of 1. This example shows the power of dimensional analysis in obtaining a non-trivial physical property of multi-particle system.

This simple result has some surprising implications. First, in three dimensions, the reaction rate is linear in $a$; it is not proportional to the cross-sectional area of the sphere. Second for $d<2$, the reaction rate increases when the radius of the absorbing sphere is decreased! This nonsensical result indicates that there is a basic problem with classic chemical kinetics for $d<2$. This pathology arises because the concentration field exterior to the absorber never reaches a steady state for $d<2$. Instead, the absorption rate of the sphere is time dependent. In contrast, for $d>2$, a steady state concentration field does arise. Once this steady state is reached, it is a trivial exercise to compute the steady-state density by solving the Laplace equation rather than the diffusion equation, and thereby obtain the steady-state flux to the sphere. From these steps, one finds the exact reaction rate

$\displaystyle k=4\pi Da\,.$

(6.3)

Armed with the basic result for the reaction rate, we now turn to a much more profound problem that is fundamental to living systems, namely, how many receptors should there be on the surface of a cell? One might expect that much of the cell surface should be covered by receptors so that its detection efficiency is high. On the other hand, one can imagine that receptors are complex and evolutionarily expensive machines. Based on cost considerations only, it would be advantageous for a cell to minimize the number of receptors. What is the appropriate balance between these two competing attributes? As a first step to address this question, we want to compute the reaction rate of a cell that is sensitive to its environment only at the locations of the receptors. We model the cell as a sphere of radius $a$ in which most of the surface is reflecting. However, on the sphere surface there are also $N$ circular domains of radius $s$ that are absorbing (Fig. 6). We view these $N$ absorbing circles as the receptors on the cell surface. What is the reaction rate of this toy model of a cell? If most of the sphere surface is reflecting, one might anticipate that the reaction efficiency of the cell will be poor. Surprisingly, the reaction efficiency of the sphere with $N$ absorbing receptors is almost as good as a perfectly absorbing sphere, even when the area fraction covered by the receptors is vanishingly small! This realization is the brilliant insight of the article by Berg and Purcell that was far ahead of its time [11]. Here I outline their argument.

The first step in their argument relies on the feature that if a diffusing particle hits the surface of a sphere, it will hit again many times before diffusing away; this point was discussed above in Sec. 3. There are two types of subsequent hitting events: (i) the particle rises a distance less than $s$ above the surface before hitting it again, and (ii) the particle rises a distance greater than $s$ above the surface before hitting it again (Fig. 6). In the former case, if the particle initially misses a receptor, it will likely miss upon the second encounter. In the latter case, if the particle misses a receptor initially, we have no information about whether the second encounter will hit or miss a receptor. Thus the rise distance $s$ demarcates the regime of dependent subsequent hits and independent subsequent hits. It is the latter events that are relevant for estimating the reaction rate. We thus use the height $s$ as the criterion for determining the number of times that a particle independently hits the cell surface before diffusing away.

When a particle is a distance $s$ above the cell surface, the probability that it eventually hits the cell again is, using Eq. (6.2),

	$\displaystyle p_{s}=\frac{a}{a+s}\,.$		(6.4a)
Thus the probability that the diffusing particle independently hits the cell $n$ times before diffusing away is $p_{s}^{n}(1-p_{s})$. Correspondingly, the average number of independent hits to the surface is
	$\displaystyle\langle n\rangle=\sum_{n>0}np_{s}^{n}\;(1-p_{s})=(1-p_{s})^{-1}=\frac{a+s}{s}\sim\frac{a}{s}\,.$		(6.4b)

The probability for a diffusing particle to not land on a receptor in a single independent hitting event is

$\displaystyle\beta=1-\frac{N\pi s^{2}}{4\pi a^{2}}\,,$

(6.4c)

namely, the area fraction of the surface that is not covered by receptors. Thus the probability $p_{\rm escape}$ that a diffusing particle that reaches the cell but never lands on a receptor is the probability that the particle always misses a receptor in each of its independent hitting attempts. This is

$\displaystyle p_{\rm escape}=\sum_{n>0}(\beta p_{s})^{n}(1-p_{s})=\frac{1-p_{s}}{1-\beta p_{s}}=\frac{4a}{4a+Ns}\,.$

(6.4d)

The probability that a diffusing particle ultimately hits a receptor thus is

$\displaystyle p_{\rm hit}=1-p_{\rm escape}=\frac{Ns}{4a+Ns}\,.$

(6.4e)

The final result for the reaction rate of a cell of radius $a$ that is covered by $N$ receptors, each of radius $s$, is

$\displaystyle k=4\pi Da\times\frac{Ns}{4a+Ns}\equiv 4\pi Da\times\text{efficiency}\,.$

(6.5)

To understand the implication of this result, let’s use some numbers that typify a cell: $a=1$ micron, $s=10$ nanometers, and $N=10^{4}$. The area fraction covered by the receptors is roughly $10^{-3}$, but the absorption efficiency of the cell is roughly 2/3! Evidently, Mother Nature is very smart to not waste resources on endowing a cell with too many receptors.

We now turn to the first-passage properties of diffusion in a two-dimensional wedge domain with absorption when the particle hits the wedge boundary. One of our motivations for studying this system is that first passage in the wedge can be mapped onto a simple diffusive capture problem in one dimension that we’ll treat in the next section. We will obtain first-passage properties in the wedge geometry both by direct solution of the diffusion equation and also, for two dimensions, in a more aesthetically pleasing fashion by conformal transformation techniques, in conjunction with the electrostatic formulation. We also present a heuristic extension of the electrostatic approach that allows us to infer, with little additional computational effort, time-dependent first-passage properties in the wedge from corresponding time-integrated properties.

What is the time dependence of the survival probability of a diffusing lamb that is hunted by $N$ diffusing lions? We define this survival probability as $S_{N}(t)$. This toy problem is most interesting in a one-dimensional geometry where all the lions are located to one side of the lamb. It is known [14] that this survival probability asymptotically decays algebraically with time,

$\displaystyle S_{N}(t)\sim t^{-\beta_{N}}\,,$

(8.1)

and the goal is to compute the decay exponent $\beta_{N}$. As we shall discuss, the decay exponent is known for $N=1$ and $N=2$ only: $\beta_{1}=\frac{1}{2}$, and $\beta_{2}=\frac{3}{4}$. For $N\geq 3$, $\beta_{N}$ grows slowly with $N$, and numerical simulations give $\beta_{3}\approx 0.91$, $\beta_{4}\approx 1.03$, and $\beta_{10}\approx 1.4$. The focus of this section is to derive $\beta_{2}=\frac{3}{4}$ by a simple geometric approach and to develop some analytical understanding of the dependence of $\beta_{N}$ on $N$ for $N\to\infty$.

Let us begin by treating a lamb that starts at $x_{0}>0$ and a single lion that starts at $x=0$. For simplicity, the diffusivities of the lamb and the lion are assumed to both equal $D$. The separation between the lamb and the lion thus diffuses with diffusion coefficient $2D$. When this separation reaches zero, the lamb has been eaten. This problem is just the classic first-passage problem on the positive infinite line, except that the diffusion coefficient is $2D$. The probability that the lamb survives until time $t$ is the same as the first-passage time being greater than $t$. This probability therefore is

$\displaystyle S_{1}(t)$

$\displaystyle=\int_{t}^{\infty}\!\frac{x_{0}}{\sqrt{8\pi Dt^{\prime 3}}}\;e^{-x_{0}^{2}/8Dt^{\prime}}\,dt^{\prime}=\text{erf}\biggl{(}\frac{x_{0}}{\sqrt{8Dt}}\biggr{)}\sim\frac{x_{0}}{\sqrt{8\pi Dt}}\,.$

(8.2)

Thus the survival probability of the lamb asymptotically decays as $t^{-1/2}$. While the lamb is sure to die, its average lifetime is infinite. Thus a single diffusing lion is not a particularly good hunter and it might starve before eating the lamb.

What happens when there are $N=2$ lions? We again assume that the lions start from the origin while the lamb starts at $x_{0}>0$, and that the diffusivities of all particles are the same. Let us label the positions of the lions as $x_{1}$ and $x_{2}$, and the position of the lamb as $x_{3}$ The lamb survives up to time $t$ if the conditions $x_{3}>x_{2}$ and $x_{3}>x_{1}$ always hold. We can give an insightful geometric interpretation of this problem by viewing the motion of the three particles on the line as equivalent to the motion of a single effective particle in three dimensions with coordinates $(x_{1},x_{2},x_{3})$. The constraints $x_{3}>x_{2}$ and $x_{3}>x_{1}$ mean that the effective particle in three-space remains to the left of the plane $x_{2}=x_{3}$ and behind the plane $x_{1}=x_{3}$ (Fig. 7(a)). This allowed region is a wedge of opening angle $\Theta$ that is defined by the intersection of these two planes. If the particle hits one of the planes, then one of the lions has eaten the lamb.