The distribution of the time elapsed before a random variable reaches a threshold value for the first time, called the first passage time (FPT) distribution, is a fundamental characteristic of stochastic processes. Here, by solving the standard macroscopic diffusion equation, we derive analytical expressions for the FPT distribution of a diffusing particle hitting a spherical object in two dimensions (2D) and three dimensions (3D) in the course of unrestricted diffusion in open space. In addition, we calculate, analytically, the angular dependence of the FPT, known as the hit distribution. The analytical results are also compared to simulations of the motions of a random walker on a discrete lattice. This topic could be of wide pedagogical interest because the FPT is important not only in physics but also in chemistry, biology, medicine, agriculture, engineering, and finance. Additionally, the central equations often appear in physics and engineering with only trivial variations, making the solution techniques widely applicable.

Topics
Heat transfer, Integral transforms, Students, Education, Engineers, Careers and professions, Scholarly publishing, Journal, Literature search, Stochastic processes

First passage processes, which are concerned with a threshold being reached for the first time by a stochastic variable, are ubiquitous.1 The most well known example is a diffusing particle reaching a threshold position, but other examples include the formation of astronomical objects,2 chemical reactions,3 biological processes like neuron firing or moth mating,4–7 rupture,8 timing precision in intracellular events,9 random search,10–12 disease spreading,13 radiotherapy planning,14 economics and finance,15,16 psychology17 as well as reliability theory.18 A fundamental characteristic of such processes is the distribution of the time t elapsed before the threshold is reached for the first time, the so-called first passage (FPT) distribution ψ ( t )⁠. Analytical expressions are readily available in textbooks for ψ ( t ) in one-dimensional (1D) diffusion processes1,2 and widely used.19,20 However, the counterparts in higher dimensions, i.e., two-dimensional (2D) and three-dimensional (3D), seem not readily available in the physics literature. Arguably, these higher-dimensional counterparts could be more useful, as the real world is rarely 1D, even under simplifications. Recent years have seen a growing interest in these topics and analytical expressions for ψ ( t )⁠, and its angular decomposition in higher dimensions is, therefore, in high demand. The calculation of ψ ( t ) in higher dimensions is not extremely difficult and can be performed using techniques well known in heat conduction theory.21,22 However, despite an extensive literature search, we have not found a publication containing systematic and detailed accounts of the results and derivation. We present them here in order to make them more accessible to all who need them.

Our main results include the analytical expressions for the FPT and its angular dependence (i.e., the hitting location distribution) in 2D and 3D diffusion processes in open space. The geometry is shown in Fig. 1. The target considered here is a circle in the 2D case or a sphere in the 3D case. Such processes are typical in, for example, intracellular matter transport,6,23,24 e.g., the transport of proteins between the cell cytoplasm and the interior of nucleus, and the collection of photons arriving from outer space or the absorption of photons by the Chlorophyll antenna of the photosynthetic apparatus.25 

Fig. 1.
Geometry for the first passage processes. A diffusing particle or a random walker starts at position x 0, which is x 0 = r for 1D, ( r , 0 ) for 2D, and ( 0 , 0 , r ) for 3D and disperses from there. The target is a segment of length 2a in 1D, a circle of radius a in 2D, and a sphere of radius a in 3D. In all cases, the target is centered at the origin. In the case of a random walker, the space-time is discretized as usual: the space is a lattice, and time is counted in steps.
View largeDownload slide

Geometry for the first passage processes. A diffusing particle or a random walker starts at position x 0⁠, which is x 0 = r for 1D, ( r , 0 ) for 2D, and ( 0 , 0 , r ) for 3D and disperses from there. The target is a segment of length 2a in 1D, a circle of radius a in 2D, and a sphere of radius a in 3D. In all cases, the target is centered at the origin. In the case of a random walker, the space-time is discretized as usual: the space is a lattice, and time is counted in steps.

Fig. 1.
Geometry for the first passage processes. A diffusing particle or a random walker starts at position x 0, which is x 0 = r for 1D, ( r , 0 ) for 2D, and ( 0 , 0 , r ) for 3D and disperses from there. The target is a segment of length 2a in 1D, a circle of radius a in 2D, and a sphere of radius a in 3D. In all cases, the target is centered at the origin. In the case of a random walker, the space-time is discretized as usual: the space is a lattice, and time is counted in steps.
View largeDownload slide

Geometry for the first passage processes. A diffusing particle or a random walker starts at position x 0⁠, which is x 0 = r for 1D, ( r , 0 ) for 2D, and ( 0 , 0 , r ) for 3D and disperses from there. The target is a segment of length 2a in 1D, a circle of radius a in 2D, and a sphere of radius a in 3D. In all cases, the target is centered at the origin. In the case of a random walker, the space-time is discretized as usual: the space is a lattice, and time is counted in steps.

Close modal

The remainder of this paper is organized as follows. In Sec. II, the standard diffusion equation is solved to establish ψ ( t ) for diffusion processes in 2D and 3D. We also revisit the result for 1D. Our derivations follow a similar method to a closely related problem on heat conduction in a circular cylinder.22 In Sec. III, we derive the angular dependence of the FPT distribution and, therefore, obtain the hit distribution. The paper is concluded in Sec. IV. Two appendices are provided. In  Appendix A, the inverse Laplace transform used for numerical calculations is discussed in detail to add to the pedagogical value of the paper. In  Appendix B, simulations are conducted of random walkers moving on simple lattices in 2D and 3D. The analytical results are shown to be well reproduced by the simulations in the limit where the discrete nature of the lattices becomes irrelevant, thereby verifying the analytical results.

The quantity that serves as a starting point for us is the probability density of a diffusing particle to be found at position x from its initial position x 0 in a time t, which we denote by G ( x , x 0 , t )⁠. For t → 0⁠, the particle has not moved, meaning that lim t → 0 G ( x , x 0 , t ) = δ d ( x − x 0 )⁠, where d is the spatial dimension and δd denotes the d-dimensional Dirac function. The subsequent evolution of G ( x , x 0 , t ) is governed by the diffusion equation,
(1)
(2)
where D is the diffusion constant of the particle and J is the probability flux density. Equation (1) is also called the equation of continuity, which expresses the conservation of total probability, while Eq. (2) is known as Fick's law.
We shall find G ( x , x 0 , t ) by solving the aforementioned diffusion equation. It is convenient to perform a Laplace transform, defined as f ̃ ( s ) = ∫ 0 ∞ d t   f ( t ) e − s t for a general function f(t), to the diffusion equation. After the transform, Eq. (1) becomes1 
(3)
Here, G ̃ is the Laplace transform of G. Note that the right-hand side of this equation arises as the initial value of G, which appears as a point source term in Eq. (3).
We aim to find the FPT distribution ψ ( r , t )⁠, whose Laplace transform is, as we denote, ψ ̃ ( r , s )⁠. Physically, ψ ( r , t ) d t gives the probability that the particle initially reaches a target centered a distance r away from its starting position x 0 during the time interval [ t , t + d t ]⁠. As aforementioned, the target is a d-dimensional sphere of radius a < r, as shown in Fig. 1. To get ψ ( t )⁠, one may solve Eq. (3) subjected to the following boundary conditions:
(4)
(5)
where ∂ Ω denotes the surface of the target. Equation (5) defines an absorbing boundary: any particles reaching the surface are instantly removed so that no particles can diffuse outward after having reached that location. Now ψ ( r , t ) amounts to the total flux into the target and is given by an integral of the flux density J over the target surface. We then find
(6)
Here, d s is the surface element (inwardly directed) of the target and J ̃ is the Laplace transform of J. In what follows, we solve Eqs. (3)–(5) and obtain ψ ̃ for 2D and 3D.
Before we do this, however, it is instructive to revisit the result for 1D. Equation (3) for d = 1 can be rewritten as
(7)
This is a linear ordinary differential equation with an inhomogeneous term on the right-hand side. The general solution can, therefore, be written as the sum of a special solution and a general solution to the homogeneous counterpart of the equation. One may verify by direct substitution that
(8)
Here, the first term represents the aforesaid special solution, which describes the propagation of a particle from its initial position in an infinite unbounded space (i.e., without the target), and the remaining terms give the general solution to the homogeneous part of the equation, with A and B being coefficients to be determined by boundary conditions. Now let us place the center of the target, which is just a segment of length 2a as shown in Fig. 1, at x = 0 and release the particle at x = x 0 = r > a⁠. Then, no particles shall be found for x < a due to the absorbing boundary. Requiring that G ̃ vanishes for x → ∞⁠, we find B = – 1. Requiring that G ̃ vanishes for x = a gives ( 1 + A ) e s / D ( x 0 − a ) − e s / D ( a − x 0 ) = 0⁠. Hence
(9)
Now the flux density is
(10)
Noting that the inward normal of the target surface is to the negative x-axis, we obtain
(11)
To get ψ ( r , t )⁠, we need to do an inverse Laplace transform, ψ ( r , t ) = ∫ c − i ∞ c + i ∞ ( d s / 2 π i ) e s t ψ ̃ ( r , s )⁠, where c is some positive number. This integral can be done analytically in various ways.1 Here, we have opted to do it by means of contour integration, since this technique is also applicable to higher dimensions with the same contour; see  Appendix A for details. We get
(12)
which is the celebrated result for 1D. One may note that here ψ ( r , t ) is just the product of the average diffusion speed ( r − a ) / t with the probability density G ( r , a , t )⁠.
In the 2D case, a particle is released at x 0 = ( r > a , 0 )⁠, and the target is a circle of radius a centered about the origin, see Fig. 1. The symmetry of the geometry suggests that polar coordinates will be useful. We, thus, write x = ( x , y ) = R ( cos   θ ,   sin   θ )⁠, where R is the length of x and θ is the polar angle (i.e., the angle with the x-axis) lying between − π and π. With ∂ x 2 = ∂ R 2 + ( 1 / R ) ∂ R + ( 1 / R 2 ) ∂ θ 2⁠, Eq. (3) now takes on the following form:
(13)
where c = s / D and P ( R , θ , s ) = D G ̃ ( x , x 0 , s )⁠. Here, we have also used that δ 2 ( x − x 0 ) = ( 1 / r ) δ ( R − r ) δ ( θ )⁠. This follows from the facts that
∫ d 2 x · δ 2 ( x − x 0 ) · f ( x ) = ∫ d θ ∫ d R R   · 1 R δ ( R − r ) δ ( θ ) · f ( x ) ,
valid for arbitrary function f ( x )⁠, and that ( 1 / R ) δ ( R − r ) = ( 1 / r ) δ ( R − r )⁠.
To make progress, we convert Eq. (13) into an ordinary differential equation involving R only. To this end, we note that δ ( θ ) is an even function and, hence, can be expanded as a cosine series,
(14)
Now, as dictated by the linearity of Eq. (13), P ( R , θ , s ) can also be expanded as a cosine series,
(15)
where P n ( R , s ) are the coefficients. The flux density along the radial direction, given by ( x / R ) · J ̃ = − ( x / R ) · ∂ x P = − ∂ R P⁠, now obtains as
(16)
where J n ( R , s ) = − ∂ R P n ( R , s )⁠. Integrating this over the entire target yields
(17)
Here, the minus sign accounts for the fact that the inward normal of the circle is oriented opposite to the radial unit vector x / R⁠. In Eq. (17), the integration averages out all the contributions from higher-order harmonics. However, the latter contains information regarding the distribution of the hit location, which is to be discussed in Sec. III.
Substituting the expressions (14) and (15) into (13) and equating the corresponding terms in the series, we find that P n ( R , s ) satisfies the following ordinary differential equation:
(18)
or, multiplying this by R2
(19)
It may be pointed out that, with some trivial modifications, this equation also describes the propagation of electromagnetic waves and matter waves in quantum mechanics. Its homogeneous part is known as the Bessel equation, to which the general solutions are linear combinations of modified Bessel functions.26 Since the inhomogeneous term vanishes everywhere except for R = r, we may write down the solutions separately for R > r and R < r. Namely,
(20)
Here, A, B, C, and D are coefficients, while In and Kn are the nth order modified Bessel functions of the first and second kind, respectively. The coefficients are determined by the following conditions:
  1. Continuity at R = r gives
    (21)

  2. In the limit R → ∞ ,   I n ( c R ) diverges, whereas K n ( c R ) vanishes. By requiring P ( R , θ , s ) to vanish in this limit, one finds A = 0.

  3. Integrating Eq. (19) over a tiny section including R = r results in a discontinuity in the radial flux density, which reflects on the fact that the particle can diffuse either toward the target or away from it. Explicitly,
    P ′ n ( r + 0 + , s ) − P ′ n ( r − 0 + , s ) = − 1 / r .
    Here, P ′ n ( R , s ) = ∂ R P n ( R , s ) and 0+ denotes the positive infinitesimal. Explicitly, this gives
    (22)

    where the prime indicates the derivative to R.

  4. At the absorbing boundary at R = a, P n ( a , s ) = 0⁠, which yields
    (23)
These conditions allow us to fully determine the coefficients. Using the identity that
(24)
and after some algebra, we obtain
(25)
By using the aforementioned identity again, one may show that
(26)
With this, we arrive at
(27)
Now an inverse Laplace transform to this expression gives ψ ( r , t )⁠. This is done by the technique of contour integration in  Appendix A, where we find
(28)
Here, J0 and Y0 are the zeroth-order Bessel functions of the first and the second kind, respectively. The integral in Eq. (28) can be readily carried out numerically.
For r ≃ a⁠, a simplification occurs to Eq. (28). Writing Y 0 ( r x / D ) ≈ Y 0 ( a x / D ) + ( r − a ) x / D Y ′ 0 ( r x / D ) and similarly for J 0 ( r x / D ) and then using the property that J 0 ( z ) Y ′ 0 ( z ) − J ′ 0 ( z ) Y 0 ( z ) = ( 2 / π z )⁠, we obtain
(29)
This roughly estimates the distribution of the time it takes for the particle to return for the first time to the location where it starts. Physically, for r ∼ a⁠, the target looks flat to the particle, and, hence, the 1D result [c.f. Eq. (12)] is recovered.
In the 3D case, the particle is released at x 0 = ( 0 , 0 , r > a )⁠, and the target is a sphere of radius a centered about the origin, see Fig. 1. We then use spherical coordinates and write x = ( x , y , z ) = R ( sin   θ   cos   ϕ ,   sin   θ   sin   ϕ ,   cos   θ )⁠, where θ ∈ [ 0 , π ] and ϕ ∈ [ 0 , 2 π ) are angles. The setup is rotationally symmetric about the z-axis. Thus, we expect that G ̃ ( x , x 0 , s ) does not depend on the azimuthal angle ϕ⁠. As in the 2D case, we introduce P ( R , θ , s ) = D G ̃ ( x , x 0 , s )⁠, which satisfies
(30)
where c = s / D as before, and we have used ∂ x 2 = ( 1 / R 2 ) ∂ R R 2 ∂ R + ( 1 / R 2 ) ( 1 / sin   θ ) ∂ θ   sin   θ ∂ θ and
(31)
This relation follows from the fact that
∫ d 3 x · δ 3 ( x − x 0 ) · f ( x ) = 2 π ∫ d θ   sin   θ ∫ d R R 2 · 1 2 π R 2   sin   θ · f ( x ) δ ( R − r ) δ ( θ )
valid for arbitrary function f ( x )⁠. Here, the factor 2 π arises from the integration over ϕ due to rotational symmetry.
Since θ ∈ [ 0 , π ]⁠, there is a one-to-one correspondence between q = cos   θ and θ. It turns out to be more convenient to work with q instead of θ directly. We shall hereafter write p ( R , q , s ) = P ( R , θ , s )⁠. Now, the source term can be written as
(32)
where we have used δ ( q − 1 ) = δ ( θ ) /   sin   θ⁠. Note that26,
(33)
where p l ( q ) are Legendre polynomials.26 In accord with the linearity of Eq. (30), we expand
(34)
where P l ( R , s ) are the coefficients satisfying an ordinary differential equation (see later). Now, the radial flux density reads − ∂ R p ( R , q , s )⁠, which is integrated over the target surface to get
(35)
Here, we have used the following identity:26 
(36)
Again, the angular information was erased by the integration and will be recovered in Sec. III.
Substituting the expansions Eqs. (33) and (34) in Eq. (30) and using (32),
(37)
In obtaining this equation, we have also used26 
(38)
Multiplying Eq. (35) by R2, we obtain
(39)
Note that this equation is also encountered when solving the Schrödinger's equation for a particle moving in a central potential in quantum mechanics. In particular, for l = 0
(40)
This equation can readily be solved. Let P 0 = u / R⁠. Then, ∂ R ( R 2 ∂ R P 0 ) = R ∂ R 2 u⁠, and the equation becomes
(41)
or, dividing it by R
(42)
This is basically the same equation we have solved for d = 1. The solution is then obtained as
(43)
which is plugged in Eq. (35) to yield
(44)
After an inverse Laplace transform,
(45)
which is the FPT distribution in 3D. The dependence on t is exactly the same as in 1D. However, the probability that the target could be reached at all, given by ∫ 0 ∞ ψ ( r , t ) d t⁠, is unity in 1D but less in 3D. Thus, in 3D, the diffusing particle may never find the target, unlike in 1 and 2D where diffusion is recurrent. This accords with the fact that in higher dimensions, the particle has a much bigger space to explore. It may be mentioned that the time-integrated probability can be computed in any dimension in terms of elementary functions by analogy with electrostatics.1 

While we have so far been focused on the FPT ψ ( t )⁠, which only gives an angularly averaged picture of the FPT in 2D and 3D, the results obtained in Sec. II also allow us to get the distribution of the hitting location, which reveals the angular dependence (i.e., the hit distribution) of the FPT and, hence, may be more desirable in some applications. These aspects are discussed in this section, where it is shown that the target is hit by the particle mostly at the heading part of the surface at short times but gradually spread out over the whole surface.

Let us begin with the 2D case and decompose the corresponding FPT distribution as
ψ ( r , t ) = ∫ 0 2 π d θ ψ ( r , θ , t ) ,
where dtd θ ψ ( r , θ , t ) gives the probability of the circle being hit where the polar angle lies in [ θ , θ + d θ ] for the first time by the particle in the time interval [ t , t + d t ]⁠. The Laplace transform of ψ ( r , θ , t )⁠, denoted by ψ ̃ ( r , θ , s )⁠, corresponds to the flux through the arc spanning from θ to θ + d θ normalized by d θ⁠, i.e., − a J ( a , θ , s ) with J being the radial flux density given by Eq. (16). Using Eq. (26),
(46)
Performing an inverse Laplace transform to this expression by means of the method described in  Appendix A, one obtains
(47)
where ψ n ( r , t ) has the same form as Eq. (28) except that the Bessel functions Y0 and J0 are replaced by Yn and Jn, respectively.
Similarly, in 3D, we decompose
ψ ( r , t ) = ∫ 0 π d θ   ψ ( r , θ , t ) ,
where dtd θ ψ ( r , θ , t ) now gives the probability of the sphere being hit within a ribbon spanning from θ to θ + d θ for the first time by the particle in the time interval [ t , t + d t ]⁠. The Laplace transform of ψ ( r , θ , t )⁠, which corresponds to the flux through the ribbon normalized by the angular span, reads
(48)
where P ′ l ( R , s ) = ∂ R P l ( R , s )⁠, and p l ( cos   θ ) are Legendre polynomials defined before in Sec. II B. In Sec. II B, we have shown that Pl satisfies the following equation:
(49)
The solution to this equation, whose homogeneous part is called the spherical Bessel function,26 can be found by writing down the general solutions for R > r and R < r and then joining them at R = r, in an analogous manner to how we solved Eq. (18). For R > r,
P l ( R , s ) = A i l ( c R ) + B k l ( c R ) ,
while for R < r
P l ( R , s ) = C i l ( c R ) + D k l ( c R ) .
Here, il and kl are the lth order modified spherical Bessel functions of the first and the second kind, respectively.26 The coefficients A, B, C, and D are determined by the boundary conditions,
P l ( R → ∞ , s ) = 0 ,   P l ( a , s ) = 0 ,   P l ( r + 0 + , s ) = P l ( r − 0 + , s ) ,
and the discontinuity at r
  P ′ l ( r + 0 + , s ) − P ′ l ( r − 0 + , s ) = − 1 2 π r 2 ,
where P ′ l = ∂ R P l⁠, and 0+ denotes the positive infinitesimal as before. Using the Wronskian
i ′ l ( x ) k l ( x ) − i l ( x ) k ′ l ( x ) = π 2 x 2 ,
we find
(50)
Hence, using the Wronskian again,
(51)
It follows that
(52)
The inverse Laplace transform of this expression can be done in the same manner as described in  Appendix A. It suffices to note that the inverse Laplace transform of k l ( s / D r ) / k l ( s / D a )⁠, denoted by ψ l ( r , t )⁠, has the same form as Eq. (28) except that the Y0 and J0 therein are replaced by spherical Bessel functions yl and jl, respectively. See that p l ( 1 ) = 1 regardless of l.

In Figs. 2 and 3, we display the hit distribution ψ ( r , θ , t ) and the cumulative hit distribution ∫ 0 θ d θ ′   ψ ( r , θ ′ , t ) (both normalized by the total FPT) for 2D and 3D, respectively. When calculating the cumulative hit distribution for 3D, we have used the formula that p l ( x ) = ( 1 / 2 l + 1 ) d / d x ( p l + 1 ( x ) − p l − 1 ( x ) )⁠. For 2D, as seen on the left panel in Fig. 2, at short times, most particles hit the target near the shortest point, and the distribution is symmetric about θ = 0. As time goes by, particles start to hit the target from elsewhere, and the distribution gets spread out. Such trend is also exhibited in the cumulative distribution as shown on the right panel. For 3D, as shown in Fig. 3, at short times, most particles also hit the target near the closest point, but the distribution is not peaked at θ = 0. Rather, ψ ( r , θ , t ) vanishes at the pole θ = 0; this is so because the area of the ribbon is proportional to sin   ( θ ) [see the text above Eq. (48)] and so is the flux through it [see also Eq. (52)]. For the same reason, the hit distribution also vanishes at the other pole θ = π⁠. This feature persists at all times. So, unlike the 2D case, ψ ( r , θ , t ) for 3D does not peak at θ = 0. Instead, as the time elapses, the position of the peak shifts from small θ toward the equator at θ = π / 2⁠. Eventually, the peak dwells at the equator. This is so because at long times, the flux density becomes evenly distributed over the entire target, and the hit distribution is solely determined by the ribbon area, which maximizes at the equator.

Fig. 2.
Hit distribution ψ ( r , θ , t ) (left panel) and cumulative hit distribution ∫ 0 θ d θ ′ ψ ( r , θ ′ , t ) (right panel), both normalized by ψ ( r , t ) = ∫ 0 2 π d θ ′ ψ ( r , θ ′ , t ), for 2D at various instants. Units of time and length are D / a 2 and a, respectively. r = 10 a.
View largeDownload slide

Hit distribution ψ ( r , θ , t ) (left panel) and cumulative hit distribution ∫ 0 θ d θ ′ ψ ( r , θ ′ , t ) (right panel), both normalized by ψ ( r , t ) = ∫ 0 2 π d θ ′ ψ ( r , θ ′ , t )⁠, for 2D at various instants. Units of time and length are D / a 2 and a, respectively. r = 10 a⁠.

Fig. 2.
Hit distribution ψ ( r , θ , t ) (left panel) and cumulative hit distribution ∫ 0 θ d θ ′ ψ ( r , θ ′ , t ) (right panel), both normalized by ψ ( r , t ) = ∫ 0 2 π d θ ′ ψ ( r , θ ′ , t ), for 2D at various instants. Units of time and length are D / a 2 and a, respectively. r = 10 a.
View largeDownload slide

Hit distribution ψ ( r , θ , t ) (left panel) and cumulative hit distribution ∫ 0 θ d θ ′ ψ ( r , θ ′ , t ) (right panel), both normalized by ψ ( r , t ) = ∫ 0 2 π d θ ′ ψ ( r , θ ′ , t )⁠, for 2D at various instants. Units of time and length are D / a 2 and a, respectively. r = 10 a⁠.

Close modal
Fig. 3.
Hit distribution ψ ( r , θ , t ) (left panel) and cumulative hit distribution ∫ 0 θ d θ ′ ψ ( r , θ ′ , t ) (right panel), both normalized by ψ ( r , t ) = ∫ 0 π d θ ′ ψ ( r , θ ′ , t ), for 3D at various instants. Units of time and length are D / a 2 and a, respectively. r = 10 a.
View largeDownload slide

Hit distribution ψ ( r , θ , t ) (left panel) and cumulative hit distribution ∫ 0 θ d θ ′ ψ ( r , θ ′ , t ) (right panel), both normalized by ψ ( r , t ) = ∫ 0 π d θ ′ ψ ( r , θ ′ , t )⁠, for 3D at various instants. Units of time and length are D / a 2 and a, respectively. r = 10 a⁠.

Fig. 3.
Hit distribution ψ ( r , θ , t ) (left panel) and cumulative hit distribution ∫ 0 θ d θ ′ ψ ( r , θ ′ , t ) (right panel), both normalized by ψ ( r , t ) = ∫ 0 π d θ ′ ψ ( r , θ ′ , t ), for 3D at various instants. Units of time and length are D / a 2 and a, respectively. r = 10 a.
View largeDownload slide

Hit distribution ψ ( r , θ , t ) (left panel) and cumulative hit distribution ∫ 0 θ d θ ′ ψ ( r , θ ′ , t ) (right panel), both normalized by ψ ( r , t ) = ∫ 0 π d θ ′ ψ ( r , θ ′ , t )⁠, for 3D at various instants. Units of time and length are D / a 2 and a, respectively. r = 10 a⁠.

Close modal

To summarize, we have derived analytical expressions for the FPT distribution of a symmetric target (circle or sphere) being hit by a diffusing particle in 2D and 3D and the corresponding hitting location distribution as well. Our results were derived with the particle released outside the target region. However, it is easy to derive the results with the particle released inside the region by the same method. We shall present these results elsewhere in connection with first passage phenomena subjected to resetting, a topic that has attracted lots of attention in recent years.11,27

The authors thank A. E. Lindsay and A. J. Bernoff for some help with the calculations. The support of the Supercomputing Wales project, which is part-funded by the European Regional Development Fund via Welsh Government, is gratefully acknowledged. The work was developed out of an undergraduate project by A. Clarkson with H.-Y. Deng.

The authors have no conflicts of interest to declare.

The data and codes that support the findings of this study are available from the corresponding author upon reasonable request.

We first do the ILT of Eq. (11). This is done as follows. With l = r − a and c being some positive number,
(A1)
where T = t D / l 2 is the dimensionless time. To evaluate the integral on the right-hand-side, let us consider the following integral in the complex s-plane:
(A2)
Here, C is a contour shown in Fig. 4, where we have chosen Arg ( s ) ∈ ( − π , π ]⁠, so that the branch cut of s lies on the negative real axis. It is composed of six pieces, C 1 , … , C 6⁠,
Fig. 4.
The contour in the complex-s plane used for doing the ILT.
View largeDownload slide

The contour in the complex-s plane used for doing the ILT.

Fig. 4.
The contour in the complex-s plane used for doing the ILT.
View largeDownload slide

The contour in the complex-s plane used for doing the ILT.

Close modal

  • C1: the line s ∈ [ c − i R , c + i R ] for some large value of R.

  • C2: the circular arc of radius R from the top of C1 to just above the negative real axis, plus the top straight segment joining the arc to C1.

  • C3: the line just above the negative real axis between [ − R , − ϵ ] for some small ϵ.

  • C4: the circular arc of radius ϵ about the origin.

  • C5: the line just below the negative real axis between [ − ϵ , − R ]⁠.

  • C6: the circular arc of radius R from just below the negative real axis to the bottom of C1, plus the bottom straight segment joining the arc to C1.

Thanks to Jordan's lemma,26 the integral vanishes on C2 and C6 in the limit R → ∞⁠. One may also show that the integral vanishes on C4 for ϵ → 0⁠. Note that the integrand has no poles inside C, and, hence, the integral (A2) vanishes according to Cauchy's theorem, yielding
(A3)
On C3, we parameterize by s = e i π x and have
(A4)
On C5, we parameterize by s = e − i π x and have
(A5)
Combining the aforementioned results,
(A6)
Dividing this by 2 π i⁠, we get the integral on the r.h.s. of Eq. (A1), i.e.,
(A7)
With x = u 2⁠, the integral over x becomes
(A8)
To do this integral, let us introduce
I ( μ ) = ∫ − ∞ ∞ d u   e − T u 2   cos   ( μ u ) = ∫ − ∞ ∞ d u   e i μ u − T u 2 = π T   exp   ( − μ 2 4 T ) .
Then, the integral (A8) equals − I ′ ( 1 )⁠, where I ′ ( μ ) = d I / d μ⁠. After some algebra, one recovers Eq. (11).
To evaluate the ILT of Eq. (27), we follow the same procedures. Now, we look at the following integral on the same contour C,
(A9)
Again, the integrals on C2 and C6 as well as C4 vanish, for K 0 ( z ) ∼ e − z / z for large argument z and K 0 ( z ) ∼ − ln   z for small argument. On C3, we get
(A10)
On C5, we get
(A11)
By Cauchy's theorem, we finally obtain the ILT of Eq. (27) as
(A12)
which, upon using K 0 ( i k ) = − ( π / 2 ) ( Y 0 ( k ) + i J 0 ( k ) ) with J0 and Y0 being the Bessel functions of the first and second kind, respectively, immediately reduces to Eq. (28).

The correctness of Eq. (28) has been directly verified numerically. In Fig. 5, we compute the ILT of Eq. (27) numerically using the method proposed by Abate and Whitt (red line),28 which agrees nicely with the analytical result Eq. (28) (dashed green line).

Fig. 5.
FPT distribution for 2D. Parameters: D = 0.025, a = 10, and r = 15, all in arbitrary units.
View largeDownload slide

FPT distribution for 2D. Parameters: D = 0.025, a = 10, and r = 15, all in arbitrary units.

Fig. 5.
FPT distribution for 2D. Parameters: D = 0.025, a = 10, and r = 15, all in arbitrary units.
View largeDownload slide

FPT distribution for 2D. Parameters: D = 0.025, a = 10, and r = 15, all in arbitrary units.

Close modal

Here, we perform numerical simulations to verify the analytical results obtained in the main text. We simulate the motion of a random walker on a simple lattice (square for 2D and cubic for 3D) with lattice constant b. Let p be the probability that the walker moves to a neighboring site in a step, and, thus, 1 − p is the probability it stays on the site where it currently resides. Let τ be the time interval for a step. One may define ν = ( p / 2 d ) ( 1 / τ ) as the attempt rate for the walker to hop to an adjacent site. Note that 2d is the coordinate number for the lattice. For ν τ ≪ 1 and b smaller than any other lengths in the system (i.e., the continuum limit), the behaviors of the random walker are expected to be the same as a diffusing particle with diffusion constant D = ν b 2⁠.

In Figs. 6–8, we display the results for 1D, 2D, and 3D, respectively. Here, time t corresponds to the number of steps. The agreement between the analytical result, Eq. (12), and the simulation is very good even for parameters for which the analytical results were not intended (see the curve with r − a = 5 b in Fig. 6). In 2D and 3D, the agreement is still very good for all parameters at long t, but there is an obvious discrepancy at short t for r and a comparable to b. Increasing r and a improves the agreement, as is clear from Figs. 7 and 8. We have performed discrete-time random walks with p ≪ 1 to average out disparity between odd and even sites. Alternatively, one can study values averaged over odd and even time steps or perform continuous-time random walks, and similar results are expected.

Fig. 6.
FPT distribution for 1D. Smooth curves stand for analytical results and fluctuating ones for simulations (106 iterations).
View largeDownload slide

FPT distribution for 1D. Smooth curves stand for analytical results and fluctuating ones for simulations (106 iterations).

Fig. 6.
FPT distribution for 1D. Smooth curves stand for analytical results and fluctuating ones for simulations (106 iterations).
View largeDownload slide

FPT distribution for 1D. Smooth curves stand for analytical results and fluctuating ones for simulations (106 iterations).

Close modal
Fig. 7.
FPT distribution for 2D. Smooth curves stand for analytical results and fluctuating ones for simulations (106 iterations). Parameters: r = 2 a and p = 0.4.
View largeDownload slide

FPT distribution for 2D. Smooth curves stand for analytical results and fluctuating ones for simulations (106 iterations). Parameters: r = 2 a and p = 0.4.

Fig. 7.
FPT distribution for 2D. Smooth curves stand for analytical results and fluctuating ones for simulations (106 iterations). Parameters: r = 2 a and p = 0.4.
View largeDownload slide

FPT distribution for 2D. Smooth curves stand for analytical results and fluctuating ones for simulations (106 iterations). Parameters: r = 2 a and p = 0.4.

Close modal
Fig. 8.
FPT distribution for 3D. Smooth curves stand for analytical results and fluctuating ones for simulations (106 iterations). Parameters: r = 2 a and p = 0.8.
View largeDownload slide

FPT distribution for 3D. Smooth curves stand for analytical results and fluctuating ones for simulations (106 iterations). Parameters: r = 2 a and p = 0.8.

Fig. 8.
FPT distribution for 3D. Smooth curves stand for analytical results and fluctuating ones for simulations (106 iterations). Parameters: r = 2 a and p = 0.8.
View largeDownload slide

FPT distribution for 3D. Smooth curves stand for analytical results and fluctuating ones for simulations (106 iterations). Parameters: r = 2 a and p = 0.8.

Close modal

We are aware that the scripts for the simulations are likely to be of interest to college education.

1.
S.
Redner
,
A Guide to First-Passage Processes
(
Cambridge U. P
.,
Cambridge, UK
,
2001
).
Google Scholar
Crossref
Search ADS
 
2.
S.
Chandrasekhar
, “
Stochastic problems in physics and astronomy
,”
Rev. Mod. Phys.
15
(
1
),
1
–
89
(
1943
).
https://doi.org/10.1103/RevModPhys.15.1
Google Scholar
Crossref
Search ADS
 
3.
R. J.
Preston
,
M. F.
Gelin
, and
D. S.
Kosov
,
J. Chem. Phys.
154
,
114108
(
2021
).
https://doi.org/10.1063/5.0045652
Google Scholar
Crossref
Search ADS
PubMed
 
4.
G. L.
Gerstein
 and
B. B.
Mandelbrot
, “
Random walk models for the spike activity of a single neuron
,”
Biophys. J.
4
(
1
),
41
–
68
(
1964
).
https://doi.org/10.1016/S0006-3495(64)86768-0
Google Scholar
Crossref
Search ADS
PubMed
 
5.
C.
Loverdo
,
O.
Bénichou
,
M.
Moreau
, and
R.
Voituriez
, “
Enhanced reaction kinetics in biological cells
,”
Nat. Phys.
4
(
2
),
134
–
137
(
2008
).
https://doi.org/10.1038/nphys830
Google Scholar
Crossref
Search ADS
 
6.
A. E.
Lindsay
,
A. J.
Bernoff
, and
M. J.
Ward
, “
First passage statistics for the capture of a Brownian particle by a structured spherical target with multiple surface traps
,”
Multiscale Model. Simul.
15
(
1
),
74
–
109
(
2017
).
https://doi.org/10.1137/16M1077659
Google Scholar
Crossref
Search ADS
 
7.
T. L.
Stepien
,
C.
Zmurchok
,
J. B.
Hengenius
,
R. M.
Caja Rivera
,
M. R.
D'Orsogna
, and
Alan E.
Lindsay
,
Appl. Sci.
2020
(
10
),
6543
.
https://doi.org/10.3390/app10186543
Crossref
Search ADS
 
8.
Z.
Hu
,
L.
Cheng
, and
B. J.
Berne
,
J. Chem. Phys.
133
,
034105
(
2010
).
https://doi.org/10.1063/1.3456556
Google Scholar
Crossref
Search ADS
PubMed
 
9.
K. R.
Ghusinga
,
J. J.
Dennehy
, and
A.
Singh
, “
First-passage time approach to controlling noise in the timing of intracellular events
,”
Proc. Natl. Acad. Sci. U. S. A.
114
(
4
),
693
–
698
(
2017
).
https://doi.org/10.1073/pnas.1609012114
Google Scholar
Crossref
Search ADS
PubMed
 
10.
G. M.
Viswanathan
 et al,
The Physics of Foraging: An Introduction to Random Searches and Biological Encounters
(
Cambridge U. P
.,
Cambridge
,
2011
).
Google Scholar
Crossref
Search ADS
 
11.
M. R.
Evans
 and
S. N.
Majumdar
, “
Diffusion with stochastic resetting
,”
Phys. Rev. Lett.
106
,
16
,
160601
(
2011
);
https://doi.org/10.1103/PhysRevLett.106.160601
Crossref
Search ADS
[PubMed]
Google Scholar
 
M. R.
Evans
 and
S. N.
Majumdar
,
J. Phys. A
44
,
435001
(
2011
).
https://doi.org/10.1088/1751-8113/44/43/435001
Crossref
Search ADS
Google Scholar
 
12.
F.
Rojo
,
C. E.
Budde
, and
H. S.
Wio
, “
Optimal intermittent search strategies
,”
J. Phys. A
42
(
12
),
125002
(
2009
);
https://doi.org/10.1088/1751-8113/42/12/125002
Crossref
Search ADS
Google Scholar
 
F.
Rojo
 et al,
J. Phys. A
43
,
345001
(
2009
).
https://doi.org/10.1088/1751-8113/43/34/345001
Crossref
Search ADS
Google Scholar
 
13.
A. L.
Lloyd
 and
R. M.
May
, “
How viruses spread among computers and people
,”
Science
292
(
5520
),
1316
–
1317
(
2001
).
https://doi.org/10.1126/science.1061076
Google Scholar
Crossref
Search ADS
PubMed
 
14.
P. G.
Hufton
,
E.
Buckingham-Jeffery
, and
T.
Galla
, “
First-passage times and normal tissue complication probabilities in the limit of large populations
,”
Sci. Rep.
10
(
1
),
8786
(
2020
).
https://doi.org/10.1038/s41598-020-64618-9
Google Scholar
Crossref
Search ADS
PubMed
 
15.
J.
Janssen
,
O.
Manca
, and
R.
Manca
,
Applied Diffusion Processes from Engineering to Finance
(
Wiley
,
2013
).
Google Scholar
Crossref
Search ADS
 
16.
J. P.
Bouchaud
 and
M.
Potters
,
Theory of Financial Risk and Derivative Pricing: From Statistical Physics to Risk Management
(
Cambridge U. P
.,
Cambridge
,
2003
).
Google Scholar
Crossref
Search ADS
 
17.
D. J.
Navarro
 and
I. G.
Fuss
, “
Fast and accurate calculations for first-passage times in Wiener diffusion models
,”
J. Math. Psychol.
53
(
4
),
222
–
230
(
2009
).
https://doi.org/10.1016/j.jmp.2009.02.003
Google Scholar
Crossref
Search ADS
 
18.
V.
Pieper
,
M.
Domin
, and
P.
Kurth
, “
Level crossing problems and drift reliability
,”
Math. Methods Oper. Res.
45
(
3
),
347
–
354
(
1997
).
https://doi.org/10.1007/BF01194784
Google Scholar
Crossref
Search ADS
 
19.
P. L.
Krapivsky
 and
S.
Redner
, “
Life and death in an expanding cage and at the edge of a receding cliff
,”
Am. J. Phys.
64
(
5
),
546
–
552
(
1996
).
https://doi.org/10.1119/1.18152
Google Scholar
Crossref
Search ADS
 
20.
S.
Redner
 and
P. L.
Krapivsky
, “
Capture of the lamb: Diffusing predators seeking a diffusing prey
,”
Am. J. Phys.
67
(
12
),
1277
–
1283
(
1999
).
https://doi.org/10.1119/1.19115
Google Scholar
Crossref
Search ADS
 
21.
H. S.
Carslaw
 and
J. C.
Jaeger
,
Condution of Heat in Solids
, 2nd ed. (
Oxford U. P
.,
London
,
1959
).
Google Scholar
 
22.
H. S.
Carslaw
 and
J. C.
Jaeger
, “
Some two-dimensional problems in conduction of heat with circular symmetry
,”
Proc. London Math. Soc.
s2-46
(
1
),
361
–
388
(
1940
).
https://doi.org/10.1112/plms/s2-46.1.361
Google Scholar
Crossref
Search ADS
 
23.
T.
Lagache
,
E.
Dauty
, and
D.
Holcman
, “
Quantitative analysis of virus and plasmid trafficking in cells
,”
Phys. Rev. E
79
(
1
),
011921
(
2009
).
https://doi.org/10.1103/PhysRevE.79.011921
Google Scholar
Crossref
Search ADS
 
24.
T.
Lagache
 and
D.
Holcman
, “
Quantifying intermittent transport in cell cytoplasm
,”
Phys. Rev. E
77
(
3
),
030901
(
2008
).
https://doi.org/10.1103/PhysRevE.77.030901
Google Scholar
Crossref
Search ADS
 
25.
R. E.
Blankenship
,
Molecular Mechanisms of Photosynthesis
, 2nd ed. (
John Wiley & Sons
,
2014
).
Google Scholar
 
26.
NIST Handbook of Mathematical Functions
, edited by
W. J.
Olver
,
D. W.
Lozier
,
R. F.
Boisvert
, and
C. W.
Clark
 (
Cambridge U. P.
,
2023
).
Google Scholar
 
27.
B.
De Bruyne
,
J.
Randon-Furling
, and
S.
Redner
, “
Optimization in first-passage resetting
,”
Phys. Rev. Lett.
125
(
5
),
050602
(
2020
).
https://doi.org/10.1103/PhysRevLett.125.050602
Google Scholar
Crossref
Search ADS
PubMed
 
28.
J.
Abate
 and
W.
Whitt
, “
A unified framework for numerically inverting Laplace transforms
,”
Informs J. Comput.
18
(
4
),
408
–
421
(
2006
).
https://doi.org/10.1287/ijoc.1050.0137
Google Scholar
Crossref
Search ADS
 
© 2024 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
Published open access through an agreement withJISC Collections