Analytical expressions for the first passage time distribution and hit distribution in two and three dimensions

The distribution of the time elapsed before a random variable reaches a threshold value for the ﬁrst time, called the ﬁrst 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 ﬁnance. Additionally, the central equations often appear in physics and engineering with only trivial variations, making the solution techniques widely applicable


I. INTRODUCTION
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][5][6][7] rupture, 8 timing precision in intracellular events, 9 random search, [10][11][12] disease spreading, 13 radiotherapy planning, 14 economics and finance, 15,16 psychology 17 as well as reliability theory. 18A 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 wðtÞ.Analytical expressions are readily available in textbooks for wðtÞ in one-dimensional (1D) diffusion processes 1,2 and widely used. 19,20However, 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 wðtÞ, and its angular decomposition in higher dimensions is, therefore, in high demand.The calculation of wðtÞ in higher dimensions is not extremely difficult and can be performed using techniques well known in heat conduction theory. 21,22However, 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. 25he remainder of this paper is organized as follows.In Sec.II, the standard diffusion equation is solved to establish wð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. 22In 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.

II. FIRST PASSAGE TIME DISTRIBUTION
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 d ðx À x 0 Þ, where d is the spatial dimension and d d denotes the d-dimensional Dirac function.The subsequent evolution of Gðx; x 0 ; tÞ is governed by the diffusion equation, @ t Gðx; x 0 ; tÞ þ @ x Á Jðx; x 0 ; tÞ ¼ 0; (1) 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Þ ¼ Ð 1 0 dt f ðtÞe Àst for a general function f(t), to the diffusion equation.After the transform, Eq. (1) becomes 1 Here, G is the Laplace transform of G.Note that the righthand 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 wðr; tÞ, whose Laplace transform is, as we denote, wðr; sÞ.Physically, wðr; tÞdt 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 þ dt.As aforementioned, the target is a d-dimensional sphere of radius a < r, as shown in Fig. 1.To get wðtÞ, one may solve Eq. (3) subjected to the following boundary conditions: Gðx 2 @X; x 0 ; sÞ ¼ 0; (5)   where @X 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 wð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 wðr; sÞ ¼ ð @X J Á ds; J ¼ ÀD@ x Gðx; x 0 ; sÞ: Here, ds 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 w 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 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 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 ! 1, we find B ¼ -1.Requiring that G vanishes for Now the flux density is Noting that the inward normal of the target surface is to the negative x-axis, we obtain To get wðr; tÞ, we need to do an inverse Laplace transform, wðr; tÞ ¼ Ð cþi1 cÀi1 ðds=2piÞe st wðr; sÞ, where c is some positive 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 spacetime is discretized as usual: the space is a lattice, and time is counted in steps.
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 which is the celebrated result for 1D.One may note that here wðr; tÞ is just the product of the average diffusion speed ðr À aÞ=t with the probability density Gðr; a; tÞ.

A. Results for 2D
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 h; sin hÞ, where R is the length of x and h is the polar angle (i.e., the angle with the x-axis) lying between Àp and p. 3) now takes on the following form: where c ¼ ffiffiffiffiffiffiffiffi s=D p and PðR; h; sÞ ¼ D Gðx; x 0 ; sÞ.Here, we have also used that d 2 ðx À x 0 Þ ¼ ð1=rÞdðR À rÞdðhÞ.This follows from the facts that valid for arbitrary function f ðxÞ, and that ð1=RÞdðR À rÞ ¼ ð1=rÞdðR À rÞ.
To make progress, we convert Eq. ( 13) into an ordinary differential equation involving R only.To this end, we note that dðhÞ is an even function and, hence, can be expanded as a cosine series, Now, as dictated by the linearity of Eq. ( 13), PðR; h; sÞ can also be expanded as a cosine series, 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 where J n ðR; sÞ ¼ À@ R P n ðR; sÞ.Integrating this over the entire target yields wðr; sÞ ¼ À ð 2p 0 dh aJða; h; sÞ ¼ a@ R P 0 ðR; sÞj R¼a : ( 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: 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. 26Since the inhomogeneous term vanishes everywhere except for R ¼ r, we may write down the solutions separately for R > r and R < r.Namely, Here, A, B, C, and D are coefficients, while I n and K n 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 (2) In the limit R ! 1; I n ðcRÞ diverges, whereas K n ðcRÞ vanishes.By requiring PðR; h; 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, Here, P 0 n ðR; sÞ ¼ @ R P n ðR; sÞ and 0 þ denotes the positive infinitesimal.Explicitly, this gives where the prime indicates the derivative to R. (4) At the absorbing boundary at R ¼ a, P n ða; sÞ ¼ 0, which yields These conditions allow us to fully determine the coefficients.
Using the identity that and after some algebra, we obtain With this, we arrive at Now an inverse Laplace transform to this expression gives wðr; tÞ.This is done by the technique of contour integration in Appendix A, where we find Here, J 0 and Y 0 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 and similarly for J 0 r ffiffiffiffiffiffiffiffi ffi x=D p and then using the property that J 0 ðzÞY 0 0 ðzÞ À J 0 0 ðzÞY 0 ðzÞ ¼ ð2=pzÞ, we obtain wðr $ a; tÞ % 2ðr À aÞ 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.

B. Results for 3D
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 h cos /; sin h sin /; cos hÞ, where h 2 ½0; p and / 2 ½0; 2pÞ 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; h; sÞ ¼ D Gðx; x 0 ; sÞ, which satisfies PðR; h; sÞ where c ¼ ffiffiffiffiffiffiffiffi s=D p as before, and we have used Þð1=sin hÞ@ h sin h@ h and This relation follows from the fact that valid for arbitrary function f ðxÞ.Here, the factor 2p arises from the integration over / due to rotational symmetry.Since h 2 ½0; p, there is a one-to-one correspondence between q ¼ cos h and h.It turns out to be more convenient to work with q instead of h directly.We shall hereafter write pðR; q; sÞ ¼ PðR; h; sÞ.Now, the source term can be written as where we have used dðq À 1Þ ¼ dðhÞ= sin h.Note that 26 where p l ðqÞ are Legendre polynomials. 26In accord with the linearity of Eq. (30), we expand pðR; q; sÞ ¼ 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 wðr; sÞ ¼ 2pa 2 ð 1 À1 dq @ R pðR; q; sÞj R¼a ¼ 2pa 2 @ R P 0 ðR; sÞj R¼a : Here, we have used the following identity: 26 1 2 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), In obtaining this equation, we have also used 26 1 sin h @ h sin h@ h p l ðqÞ ð Þ¼À lðl þ 1Þp l ðqÞ: Multiplying Eq. ( 35) by R 2 , we obtain Note that this equation is also encountered when solving the Schr€ odinger's equation for a particle moving in a central potential in quantum mechanics.In particular, for l ¼ 0 This equation can readily be solved.Let P 0 ¼ u=R.Then, @ R ðR 2 @ R P 0 Þ ¼ R@ 2 R u, and the equation becomes 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 Ð 1 0 wðr; tÞdt, 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

III. HIT DISTRIBUTION
While we have so far been focused on the FPT wð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 wðr; tÞ ¼ ð 2p 0 dhwðr; h; tÞ; where dtdhwðr; h; tÞ gives the probability of the circle being hit where the polar angle lies in ½h; h þ dh for the first time by the particle in the time interval ½t; t þ dt.The Laplace transform of wðr; h; tÞ, denoted by wðr; h; sÞ, corresponds to the flux through the arc spanning from h to h þ dh normalized by dh, i.e., ÀaJða; h; sÞ with J being the radial flux density given by Eq. ( 16).Using Eq. ( 26 (48) where P 0 l ðR; sÞ ¼ @ R P l ðR; sÞ, and p l ðcos hÞ are Legendre polynomials defined before in Sec.II B. In Sec.II B, we have shown that P l satisfies the following equation: 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Þ ¼ Ai l ðcRÞ þ Bk l ðcRÞ; while for R < r P l ðR; sÞ ¼ Ci l ðcRÞ þ Dk l ðcRÞ: Here, i l and k l are the lth order modified spherical Bessel functions of the first and the second kind, respectively. 26The coefficients A, B, C, and D are determined by the boundary conditions, P l ðR !1; sÞ ¼ 0; P l ða; sÞ ¼ 0; and the discontinuity at r where P 0 l ¼ @ R P l , and 0 þ denotes the positive infinitesimal as before.Using the Wronskian i 0 l ðxÞk l ðxÞ À i l ðxÞk 0 l ðxÞ ¼ p 2x 2 ; we find Hence, using the Wronskian again, It follows that 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 , denoted by w l ðr; tÞ, has the same form as Eq. ( 28) except that the Y 0 and J 0 therein are replaced by spherical Bessel functions y l and j l , respectively.See that p l ð1Þ ¼ 1 regardless of l.
In Figs. 2 and 3, we display the hit distribution wðr; h; tÞ and the cumulative hit distribution Ð h 0 dh 0 wðr; h 0 ; 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=2l þ 1Þd=dxð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 h ¼ 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 h ¼ 0. Rather, wðr; h; tÞ vanishes at the pole h ¼ 0; this is so because the area of the ribbon is proportional to sin ðhÞ [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 h ¼ p.This feature persists at all times.So, unlike the 2D case, wðr; h; tÞ for 3D does not peak at h ¼ 0. Instead, as the time elapses, the position of the peak shifts from small h toward the equator at h ¼ p=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.

IV. CONCLUSIONS
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,27PPENDIX A: INVERSE LAPLACE TRANSFORM (ILT) OF EQUATIONS ( 11) AND (27)   We first do the ILT of Eq. ( 11).This is done as follows.With l ¼ r À a and c being some positive number, Here, C is a contour shown in Fig. 4, where we have chosen Arg ðsÞ 2 ðÀp; p, so that the branch cut of ffiffi s p lies on the negative real axis.It is composed of six pieces, C 1 ; …; C 6 , • C 1 : the line s 2 ½c À iR; c þ iR for some large value of R. Thanks to Jordan's lemma, 26 the integral vanishes on C 2 and C 6 in the limit R ! 1.One may also show that the integral vanishes on C 4 for !0. Note that the integrand has no poles inside C, and, hence, the integral (A2) vanishes according to Cauchy's theorem, yielding ð   : Then, the integral (A8) equals ÀI 0 ð1Þ, where I 0 ðlÞ ¼ dI=dl.
To evaluate the ILT of Eq. ( 27), we follow the same procedures.Now, we look at the following integral on the same contour Again, the integrals on C 2 and C 6 as well as C 4 vanish, for K 0 ðzÞ $ e Àz = ffiffi z p for large argument z and K 0 ðzÞ $ Àln z for small argument.On C 3 , we get (A12) which, upon using K 0 ðikÞ ¼ À p=2ÞðY 0 ðkÞ þ iJ 0 ðkÞÞ with J 0 and Y 0 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).

APPENDIX B: NUMERICAL SIMULATIONS
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 s be the time interval for a step.One may define ¼ ðp=2dÞð1=sÞ as the attempt rate for the walker to hop to an adjacent site.Note that 2d is the coordinate number for the lattice.For s ( 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 ¼ 5b 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.
We are aware that the scripts for the simulations are likely to be of interest to college education.
where T ¼ tD=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:

• C 2 :
the circular arc of radius R from the top of C 1 to just above the negative real axis, plus the top straight segment joining the arc to C 1 .•C 3 : the line just above the negative real axis between ½ÀR; À for some small .• C 4 : the circular arc of radius about the origin.•C 5 : the line just below the negative real axis between ½À; ÀR. • C 6 : the circular arc of radius R from just below the negative real axis to the bottom of C 1 , plus the bottom straight segment joining the arc to C 1 .

Fig. 4 . 1 À1
Fig. 4. The contour in the complex-s plane used for doing the ILT.
n ðR;sÞ ¼K n ðcRÞK n ðcrÞ The Laplace transform of wðr; h; tÞ, which corresponds to the flux through the ribbon normalized by the angular span, reads wðr; h; sÞ ¼ 2pa 2 sin h @ R PðR; h; sÞj R¼a 0 dh wðr; h; tÞ;where dtdhwðr; h; tÞ now gives the probability of the sphere being hit within a ribbon spanning from h to h þ dh for the first time by the particle in the time interval ½t; t þ dt. ffi