From single-particle stochastic kinetics to macroscopic reaction rates: fastest first-passage time of N random walkers

Physical kinetics studies the dynamics of molecular chemical reactions in terms of the time dependence of the reactant concentrations [1]. If only one reactant is present with concentration [A] (unimolecular reaction), the associated first-order rate equation is typically written in the form d[A]/dt = −k[A] with the reaction rate k. The resulting dynamics is exponential, [A] = [A]₀ exp(−kt), with initial concentration [A]₀ and half-life time t_1/2 = ln(2)/k. A physical derivation of chemical reaction rates in terms of the diffusivity of the molecular reactants was given in the seminal 1916 paper by von Smoluchowski [2] for immediate coagulation. For instance, the Smoluchowski rate of a particle with diffusion coefficient D to hit an immobile spherical target of radius a reads k_S = 4πDa[A]₀. A more general approach to calculate the reaction rate combining reaction and diffusion limitation was formulated by Collins and Kimball [3]. Different diffusion mechanisms can hereby lead to an effective renormalisation of the target size in these theories. Thus, when DNA-binding proteins search for their target site on a DNA chain, the combination of three-dimensional diffusion with one-dimensional diffusion along the DNA causes the target size a to be replaced by an 'effective sliding length'. In the 'facilitated diffusion picture' this quantity is based on the typical distance covered by the protein during sliding while it intermittently binds to the DNA [4–7]. Such predictions and their generalisations can by-now be measured routinely in single-molecule assays [8–14].

In chemical rate-based formulations of reaction kinetics such as those by Smoluchowski or Collins and Kimball it is tacitly assumed that the reactants are present at sufficiently high abundance, thus effecting smooth concentration levels. Moreover, in a 'well-stirred' reaction container spatial coordinates can be neglected. By-now the understanding is that such theories provide an adequate description of the chemical kinetics for most systems based on the interplay of molecular reaction and diffusion—upon some refinements and generalisations incorporating various physical factors missed in the original works [15–24]. Exceptions are provided by several particular reaction schemes in which, under some rather restrictive constraints, so-called fluctuation-induced behaviour emerges, see, e.g., [25–34] and references therein. One particular question concerns the possible concentration-dependence of the effective reaction rates [15, 35–38]. Indeed, the original Smoluchowski approach and many of its generalisations are only plausible for sufficiently low albeit finite concentrations. At higher concentrations diffusive transport as the rate-controlling factor becomes less important, and concurrently the diffusion coefficients themselves acquire a dependence on the concentrations of reactants and products.

Conversely, in a shift of general interest towards understanding the kinetic behaviour of reactions in diverse biochemical and biophysical systems low concentrations are increasingly considered to be a salient feature. Indeed, many intracellular processes of signalling, regulation, infection, immune reactions, and metabolism as well as transmitter release in neurons occur upon the arrival of one or few biomolecules to small specific regions [39, 40]. The well-stirredness assumption based on simple diffusion models leads to the conclusion that DNA-binding proteins rapidly homogenise in the cell. Experiments show, however, that even in relatively small bacterial cells such transcription factors are localised around their encoding gene and further inhomogeneity is caused by the nucleoid state [41]. Moreover, genes that are controlled by a specific transcription factor tend to locate next to the gene encoding this transcription factor [42, 43]. This is consistent with explicit models for intracellular gene regulation [44, 45]. Such gene regulatory signals often rely on nanomolar concentrations [46, 47], similar to molecular concentration levels of autoinducer molecules controlling the state of cell colonies in quorum sensing [48–52]. This causes strong fluctuations of regulation events [44, 53–55].

The nanomolar concentration range of relevant molecular species renders the concept of reaction rates rather ill-defined. Due to the lack of a sufficiently large number of molecules, reaction times become strongly defocused and one cannot describe the system in terms of a single time scale associated with a reaction, but rather in terms of the full distribution of random times to a reaction event [56–58, 68–71]. Indeed, the shortest relevant time scale of the associated first-passage of molecules to their reaction target in such situations is 'geometry-controlled' [57, 70] in terms of 'direct paths' from the molecules' initial position to its target. This shortest characteristic time scale is also the most probable in the reaction time distribution [56, 57]. The longest time scale, typically orders of magnitude longer than the most probable time scale, corresponds to 'indirect' trajectories that on their path hit the boundary of the confined volume and thus lose any signature of their original position [57, 70].

Most approaches determining reaction rates or full probability densities of reaction times rely on a single-particle scenario, yet, despite of the small concentrations we alluded to above, typically a given number of particles are searching for a common target in parallel, for instance, several transcription factors seeking to bind to a specific binding site on the cellular DNA. For particles searching for a common target in parallel two relevant questions can be asked: (i) how much faster is the search process when more than one searcher is present and (ii) which searcher comes first, for instance, when we think of two different species of transcription factors competing for the same binding site. We could also think of an entire colony of cells, for instance, the many thousands of cells in a bacteria colony [72, 73], competing with each other for which cell is able to react first to a common environmental challenge. On a different level such competitive scenarios come into play when sperm cells race towards the egg cell. Indeed, such 'particle number' effects on the reaction kinetics of few-molecule diffusion-limited reactions have recently been addressed [74–76]. Generally, there is a range of systems in molecular and cellular biology for which the number of species involved lies in some intermediate range—much more than a few, but still much less than a macroscopic number (i.e., of the order of Avogadro's number N_A ≈ 6 × 10²³) [77]. In particular, neuronal connections often occur on a dendritic spine, where calcium is present in big amounts and arrival of the fastest of the calcium ions can trigger a transduction [78]. In another instance, the post-synaptic current is generated when the first receptor of a neurotransmitter is activated, a process which is mediated by the release and transport of several thousands of neurotransmitters from the pre-synaptic terminal. In immunology, which hosts a variety of such examples, a gene mechanism responsible for the selection and expression of a specific membrane receptor on B-cells relies on many such receptors, which bind directly to recognise a molecular unit of a pathogen.

A mathematical model in which N Brownian particles start from the same position simultaneously and search for an immobile small target was first analysed by Weiss et al almost 40 years ago [79]. Various additional aspects of this problem were subsequently considered [80–91]. Specifically it was already shown by Weiss et al that the mean first arrival time of the fastest of the N Brownian searchers to the target is inversely proportional to ln N as N → ∞. This means that the change of the 'efficiency' of a reaction is quite modest as compared to the quite large investment in requiring a larger number of searchers. It was thus concluded in [79] that the speedup due to many searchers is a comparatively minor effect: namely, to reduce the reaction time scale in a noticeable way, one needs a very large amount of searchers. For instance, to have a reduction by just a factor of 10, more than 20 000 searchers need to be deployed, an expensive number for many biological scenarios. Indeed, it was recently argued in the context of the tens of millions of sperm cells competing for a single egg in higher mammals that when N such searchers are launched from the same initial location to find a given target, the decisive time scale is the arrival time of the first searcher, that is, the shortest time [74–76].

In this paper, we address the reaction dynamics between an immobile small target site and N searchers diffusing inside a bounded domain Ω of finite volume |Ω|. A major part of our analysis pertains to a very general geometry of the system such that a bounded domain may have an arbitrary shape, with the only constraint being that the boundary is smooth and thus having a finite area |∂Ω|. The target Γ with an area |Γ| and a characteristic extent R can be placed either on the boundary or in the bulk. For illustrative purposes, we will use some simple domains.

We take here a broader perspective beyond the mean shortest time associated with the arrival of the fastest among N Brownian searchers. Namely, we analyse the full distribution function of the time of the first reaction event. This comprises three different relevant aspects. First, as it is already known from the single-searcher scenario that the full distribution of reaction times spans several distinct characteristic time scales [56–58, 70] ranging from the above-mentioned most probable time, a crossover time scale from a hump-like region to a plateau-like regime in which all values of the first passage times are nearly equally probable, and ultimately, the mean first passage time to the reaction event, followed by an exponential decay. When N searchers operate in parallel, we establish the full probability density of the fastest first-passage time (fFPT).

Second, compared to previous works we include imperfect reactions characterised by a finite intrinsic reaction constant κ (with κ = ∞ corresponding to a perfect reaction). Namely, for a chemical reaction to be successful, it is not sufficient for the diffusing molecule just to arrive to the reaction centre, but a reaction activation barrier needs to be overcome [3, 18, 59–67]. On top of this many biomolecules present a specific binding area, further reducing the probability of immediate reaction on encounter. As a consequence repeated collisions with the target are required, leading to repeated excursions in the volume. The effected further defocusing of the reaction times, analysed in detail for the case N = 1 in [56–58], was partially studied in [84, 85] for searchers starting from the same location. Here we highlight some additional features.

The third and the most important difference to most previous works is that we consider the scenario in which the searching particles initially are placed at distinct random positions. The characteristic properties are then obtained by averaging over these fixed initial positions, which are supposed to be uniformly spread across the confining domain. Recall that for a target placed away from a boundary, in the thermodynamic limit when both the number of particles N and the volume |Ω| of the confining domain Ω tend to infinity while the concentration N/|Ω| is kept finite (neither very small nor too large), the Smoluchowski approach (see appendix A) provides an exact solution in the case of a perfect reaction [92–96]. As shown in [97, 98] this even holds for the case of imperfect reactions when the finite reactivity is modelled in terms of a stochastic Poisson gating process. In this sense it can be argued that our analysis provides a connection between single-molecule reaction kinetics and standard chemical kinetics based on effective reaction rates even in the case when targets are located on the boundary. We also stress the conceptual difference in dealing with the uniform initial distribution as compared to a fixed starting position for all searchers (see also [88] for a related discussion). In both cases, the randomness of reaction times follows from random realisations of the particle trajectories. However, an additional source of randomness comes into play due to the random initial placement of each particle within the confining domain. This creates an intriguing new aspect to the problem, and we analyse its impact on the reaction kinetics resorting to an analysis of the typical behaviour, the averaged behaviour, and quantify fluctuations around the averaged behaviour. A uniform initial condition leads to a drastically different behaviour as compared to the case of a fixed starting point considered previously in [81–85]. In particular, as first noticed in [79] and later explored in [86–88], in a finite domain the contribution to the effective rate due to a diffusive search for a target decreases in proportion to 1/N², while the contribution due to a penetration through a barrier against reaction vanishes as 1/N in the limit N → ∞. This means that (i) placing N searchers at distinct positions gives a substantial increase in the reaction efficiency, which can be orders of magnitude larger as compared to the case when all N searchers start from the same point, (ii) as compared to a standard Collins–Kimball relation, which is valid in the thermodynamic limit and in which both contributions have the same 1/N-dependence on the number of searchers, here the contribution due to a diffusive search acquires an additional power of a concentration of searchers, (which is a strong concentration effect), and hence, (iii) in the limit N ≫ 1 reactions become inevitably controlled by chemistry (kinetically-controlled reactions) rather than by diffusion (diffusion-controlled reactions). Overall, the analysis developed here provides a comprehensive insight into the binding kinetics of the extremes of first-passage phenomena in finite systems with multiple diffusing reactants competing for a reactive target.

The paper is organised as follows. In section 2 we discuss the implications of fixed versus uniform starting positions. Section 3 contains our main theoretical results including a brief overview of general properties for bounded domains (section 3.1), approximations for the mean first-reaction time (section 3.2) and the variance (section 3.3), the long-time behaviour of the volume-averaged probability density of the reaction time (section 3.4), the fluctuations between individual realisations (section 3.5), and the typical reaction time density (section 3.7). Section 4 is devoted to the discussion of these results and their implications to chemical physics and biological systems. Details of derivations and some additional analyses are presented in the appendices.

Before proceeding to the results we discuss the uniform initial condition analysed in this paper. In the conventional macroscopic description of chemical kinetics, the reaction rate is computed by averaging the probability density $\rho \left(t\vert {\mathbf{x}}_{0}\right)$ of the first-passage time (FPT) to the target from the starting point x₀,

$$\begin{equation}J\left(t\right)=\underset{{\Omega}}{\int }\mathrm{d}{\mathbf{x}}_{0}\enspace {c}_{0}\left({\mathbf{x}}_{0}\right)\enspace \rho \left(t\vert {\mathbf{x}}_{0}\right),\end{equation} \tag{ 1 }$$

where ${c}_{0}\left({\mathbf{x}}_{0}\right)$ is the initial concentration profile of particles. As we assume a uniform initial concentration in the bounded confining domain Ω with finite volume |Ω|, the reaction rate is simply proportional to the volume-averaged probability density, $J\left(t\right)\propto \overline{\rho \left(t\right)}$, with

$$\begin{equation}\overline{\rho \left(t\right)}=\frac{1}{\vert {\Omega}\vert }\underset{{\Omega}}{\int }\mathrm{d}{\mathbf{x}}_{0}\rho \left(t\vert {\mathbf{x}}_{0}\right).\end{equation} \tag{ 2 }$$

As a consequence, the reaction rate J(t) incorporates two intertwined sources of randomness: a random choice of the initial position x₀ and a random trajectory from x₀ to the target. Even though both affect the FPT distribution, their respective roles are in fact not that well understood, thus deserving a more specific investigation.

A volume-averaged description as entering equation (2) is justified whenever the number of diffusing particles is macroscopically large so that one can actually speak about concentrations. In many of the biologically relevant settings discussed above, however, the number of diffusing particles can be small or moderately large (say, a few tens or a few thousand). One may therefore question whether the above macroscopic approximation is still applicable. In particular, what is the role of stochastic fluctuations between different realisations of the starting points? This question becomes particularly relevant in the short-time limit. In fact, if there is a single searcher, the density $\rho \left(t\vert {\mathbf{x}}_{0}\right)$ decays exponentially fast at short times, as e^−c/t, where c is proportional to the squared distance to a target from x₀. In contrast, the volume average in equation (2) superimposes all $\rho \left(t\vert {\mathbf{x}}_{0}\right)$ including the contributions from those particles that started infinitely close to the target are dominant. As a result, $\overline{\rho \left(t\right)}$ exhibits not an exponential but a power-law behaviour as t → 0, as we detail below. One can conclude that $\overline{\rho \left(t\right)}$ is not representative of an actual behaviour here, as one could expect. But what happens for N particles, in particular, in the large N limit? In other words, we aim at uncovering the role of stochastic fluctuations related to random initial positions of the particles as we change N.

In what follows we investigate the gradual transition from the single-particle description, in which stochastic fluctuations are significant, to the macroscopic description. We consider N independent identical particles undergoing Brownian dynamics with diffusion coefficient D that search in parallel for a common target (figure 1). The FPT τ₁ for a single particle, started from x₁, is characterised by its survival probability $S\left(t\vert {\mathbf{x}}_{1}\right)={\mathbb{P}}_{{\mathbf{x}}_{1}}\left\{{\tau }_{1}{ >}t\right\}$, from which the probability density is obtained by differentiation,

$$\begin{equation}\rho \left(t\vert {\mathbf{x}}_{1}\right)=-\frac{\partial }{\partial t}S\left(t\vert {\mathbf{x}}_{1}\right).\end{equation} \tag{ 3 }$$

The independence of the diffusing particles immediately implies that the FPT ${\mathcal{T}}_{N}=\mathrm{min}\left\{{\tau }_{1},\dots ,{\tau }_{N}\right\}$ of the fastest particle among N, the so-called fFPT, follows from the N-particle survival probability

$$\begin{equation}{S}_{N}\left(t\vert {\mathbf{x}}_{1},\dots ,{\mathbf{x}}_{N}\right)=\mathbb{P}\left\{{\mathcal{T}}_{N}{ >}t\right\}=\prod _{n=1}^{N}S\left(t\vert {\mathbf{x}}_{n}\right).\end{equation} \tag{ 4 }$$

The knowledge of this survival probability for a single particle therefore fully determines the one for multiple non-interacting particles searching in parallel, and the problem of characterising the fFPT may look trivial at first thought.

Zoom In Zoom Out Reset image size

However, an exact explicit form of the survival probability S(t|x) is known for a very limited number of simple settings such as, for instance, the FPT to the endpoints of an interval or to the boundary of a sphere [99]. In turn, only the asymptotic behaviour or some approximate forms are known for most practically relevant cases such as the narrow escape problem [100–109] (see also the review [110]). Moreover, even if S(t|x) is known explicitly, finding the moments of the fFPT remains a challenging and quite involved problem. While most former studies of these extreme FPTs focused on the case when all particles start from the same fixed point (i.e., x₁ = ... = x_N) [74, 76, 79, 83–85] we here explore a different direction, namely, the case when the particles start from independent and uniformly distributed points. In other words, we aim at uncovering the role of stochastic fluctuations related to random initial positions of the particles. Qualitatively, as N increases, the particular random realisation of the starting points is expected to become irrelevant, and the macroscopic description should become increasingly accurate. Here, we investigate the transition from the single-particle setting to the macroscopic limit and address the practically important question of when and how such a description becomes applicable.

To develop the theoretical framework we consider a given realisation of starting points x₁, ..., x_N. Then the probability density function (PDF) of the fFPT to a target domain Γ follows from equation (4) in the form

$$\begin{align}\hfill {\rho }_{N}\left(t\vert {\mathbf{x}}_{1},\dots ,{\mathbf{x}}_{N}\right)& =-\frac{\partial }{\partial t}\left[S\left(t\vert {\mathbf{x}}_{1}\right){\times}\cdots {\times}S\left(t\vert {\mathbf{x}}_{N}\right)\right]=S\left(t\vert {\mathbf{x}}_{1}\right){\times}\cdots {\times}S\left(t\vert {\mathbf{x}}_{N}\right)\sum _{n=1}^{N}\frac{\rho \left(t\vert {\mathbf{x}}_{n}\right)}{S\left(t\vert {\mathbf{x}}_{n}\right)}.\hfill \end{align} \tag{ 5 }$$

If (x₁, ..., x_N) are independent and uniformly distributed points in the bounded confining domain Ω, the volume-averaged density reads

$$\begin{align}\hfill \overline{{\rho }_{N}\left(t\right)}& =\underset{{\Omega}}{\int }\frac{\mathrm{d}{\mathbf{x}}_{1}}{\vert {\Omega}\vert }\dots \underset{{\Omega}}{\int }\frac{\mathrm{d}{\mathbf{x}}_{N}}{\vert {\Omega}\vert }\rho \left(t\vert {\mathbf{x}}_{1},\dots ,{\mathbf{x}}_{N}\right)=-\frac{\mathrm{d}}{\mathrm{d}t}{\left[\overline{S\left(t\right)}\right]}^{N}=N{\left[\overline{S\left(t\right)}\right]}^{N-1}\enspace \overline{\rho \left(t\right)},\hfill \end{align} \tag{ 6 }$$

where

$$\begin{equation}\overline{S\left(t\right)}=\frac{1}{\vert {\Omega}\vert }\underset{{\Omega}}{\int }\mathrm{d}\mathbf{x}\enspace S\left(t\vert \mathbf{x}\right)\end{equation} \tag{ 7 }$$

is the volume-averaged survival probability.

3.1. Summary of general single-particle properties

For a bounded domain Ω of an arbitrary shape with a smooth boundary ∂Ω the survival probability admits in the most general case the spectral expansion [99, 111, 112]

$$\begin{equation}S\left(t\vert {\mathbf{x}}_{0}\right)=\sum _{n}{\mathrm{e}}^{-Dt{\lambda }_{n}}{u}_{n}\left({\mathbf{x}}_{0}\right)\underset{{\Omega}}{\int }\mathrm{d}\mathbf{x}\enspace {u}_{n}^{{\ast}}\left(\mathbf{x}\right),\end{equation} \tag{ 8 }$$

where the asterisk denotes the complex conjugate, λ_n are non-negative eigenvalues which are sorted in ascending order,

$$\begin{equation*}0{\leqslant}{\lambda }_{1}{\leqslant}{\lambda }_{2}{\leqslant}\dots {\leqslant}{\lambda }_{n}{\leqslant}\dots ,\end{equation*}$$

and the u_n(x) are L₂(Ω)-normalised eigenfunctions of the Laplace operator, −Δu_n = λ_nu_n, subject to appropriate boundary conditions. The eigenvalues and eigenfunctions encode all the information about the geometry of the domain, the location of the target and its size [113]. We consider three typical situations (in increasing order of generality):

• (a)  
  a fully reactive boundary described by the Dirichlet boundary condition (u_n)|_∂Ω = 0 (here, Γ = ∂Ω);
• (b)  
  a partially reactive boundary described by the Robin boundary condition ${\left(D{\partial }_{\mathbf{n}}{u}_{n}+\kappa {u}_{n}\right)}_{\vert \partial {\Omega}}=0$, where κ is the reactivity and ∂_n is the normal derivative oriented outwards the domain (here, Γ = ∂Ω); and
• (c)  
  a partially reactive target Γ on the otherwise reflecting boundary ∂Ω\Γ, described by the mixed boundary condition ${\left(D{\partial }_{\mathbf{n}}{u}_{n}+\kappa {u}_{n}\right)}_{\vert {\Gamma}}=0$ and ${\left({\partial }_{\mathbf{n}}{u}_{n}\right)}_{\vert \partial {\Omega}{\backslash}{\Gamma}}=0$. We emphasise that the last situation includes two distinct cases: a target located on the reflecting boundary and a target located in the bulk (figure 1). In the latter case, the boundary of the confining domain Ω has two disjoint components: the outer reflecting boundary and the inner reactive target. Even though these two locations of the target are usually distinguished in physical literature, their mathematical description is the same. Note that the target region Γ does not need to be connected, i.e., one can consider multiple targets.

The volume-averaged survival probability then reads

$$\begin{equation}\overline{S\left(t\right)}=\sum _{n}{c}_{n}\enspace {\mathrm{e}}^{-D{\lambda }_{n}t},\end{equation} \tag{ 9 }$$

where

$$\begin{equation}{c}_{n}=\frac{1}{\vert {\Omega}\vert }{\left\vert \underset{{\Omega}}{\int }\mathrm{d}\mathbf{x}\enspace {u}_{n}\left(\mathbf{x}\right)\right\vert }^{2}.\end{equation} \tag{ 10 }$$

Similar spectral expansions hold for the PDF and its volume average,

$$\begin{equation}\rho \left(t\vert {\mathbf{x}}_{0}\right)=\sum _{n}D{\lambda }_{n}\enspace {\mathrm{e}}^{-Dt{\lambda }_{n}}{u}_{n}\left({\mathbf{x}}_{0}\right)\underset{{\Omega}}{\int }\mathrm{d}\mathbf{x}\enspace {u}_{n}^{{\ast}}\left(\mathbf{x}\right)\end{equation} \tag{ 11 }$$

and

$$\begin{equation}\overline{\rho \left(t\right)}=\sum _{n}D{\lambda }_{n}{c}_{n}\enspace {\mathrm{e}}^{-D{\lambda }_{n}t}.\end{equation} \tag{ 12 }$$

From these expansions one can easily derive the volume-averaged moments of the FPT ${\mathcal{T}}_{1}$,

$$\begin{equation}\overline{\langle {\mathcal{T}}_{1}^{k}\rangle }=k!\sum _{n}\frac{{c}_{n}}{{\left(D{\lambda }_{n}\right)}^{k}}\quad \left(k=1,2,3,\dots \right).\end{equation} \tag{ 13 }$$

Throughout the paper we distinguish the volume average over random initial points (denoted by the overline) and the ensemble average over random trajectories (denoted by angular brackets).

At long times the survival probability and the PDF decay exponentially fast, with decay rate 1/(Dλ₁) determined by the smallest eigenvalue λ₁. The volume average does not affect this decay. In contrast, the short-time behaviour is drastically different for fixed-point and volume-averaged quantities. In fact, the volume-averaged survival probability behaves as

$$\begin{equation}\overline{S\left(t\right)}\simeq 1-\frac{2\vert {\Gamma}\vert }{\sqrt{\pi }\vert {\Omega}\vert }\enspace \sqrt{Dt}+O\left(t\right)\quad \left(\kappa =\infty \right),\end{equation} \tag{ 14a }$$

$$\begin{equation}\overline{S\left(t\right)}\simeq 1-\frac{\vert {\Gamma}\vert \kappa }{\vert {\Omega}\vert }\enspace t+O\left({t}^{3/2}\right)\quad \enspace \quad \left(\kappa {< }\infty \right),\end{equation} \tag{ 14b }$$

which follows from the short-time expansion of the heat content [114–117] (see also appendix B for the next-order correction). Qualitative arguments behind this asymptotic result are simple. In the diffusion-limited regime (κ = ∞) only those particles in a thin layer of width $\sqrt{Dt}$ near the reactive boundary Γ can reach this boundary and thus disappear within short time. The relative fraction of such particles is $\sqrt{Dt}\enspace \vert {\Gamma}\vert /\vert {\Omega}\vert$, where |Γ| is the surface area of Γ and |Ω| is the volume of the domain Ω. In contrast, in the reaction-limited regime, the decay of $\overline{S\left(t\right)}$ is limited by κ and is thus proportional to t. An immediate consequence of equation (14) is

$$\begin{equation}\overline{\rho \left(t\right)}\simeq \frac{\vert {\Gamma}\vert }{\vert {\Omega}\vert }{\times}\begin{cases}\frac{1}{\sqrt{\pi }}\sqrt{D/t}+O\left(1\right),\quad \hfill & \kappa =\infty ,\hfill \\ \kappa +O\left(\sqrt{t}\right),\quad \hfill & \kappa {< }\infty .\hfill \end{cases}\end{equation} \tag{ 15 }$$

In contrast, $S\left(t\vert {\mathbf{x}}_{0}\right)$ and $\rho \left(t\vert {\mathbf{x}}_{0}\right)$ exhibit a fast exponential decay as t → 0,

$$\begin{equation}1-S\left(t\vert {\mathbf{x}}_{0}\right)\simeq A{t}^{-\alpha }\enspace {\mathrm{e}}^{-\vert {\mathbf{x}}_{0}-\partial {\Omega}{\vert }^{2}/\left(4Dt\right)},\end{equation} \tag{ 16a }$$

$$\begin{equation}\rho \left(t\vert {\mathbf{x}}_{0}\right)\simeq \alpha A{t}^{-\alpha -1}\enspace {\mathrm{e}}^{-\vert {\mathbf{x}}_{0}-\partial {\Omega}{\vert }^{2}/\left(4Dt\right)},\end{equation} \tag{ 16b }$$

with the constants A and α, see [70, 76, 85, 118–120].

As shown in [88], the above short-time asymptotic behaviour of the survival probability for a single particle determines the leading behaviour of the moments of the fFPT and its volume-averaged probability density $\overline{{\rho }_{N}\left(t\right)}$ in the limit N → ∞. In particular, for a uniform distribution of the starting points, $\overline{{\rho }_{N}\left(t\right)}$ is close to the Weibull density

$$\begin{equation}\overline{{\rho }_{N}\left(t\right)}\simeq kBN{t}^{k-1}\enspace {\mathrm{e}}^{-BN{t}^{k}}\quad \left(N\to \infty \right),\end{equation} \tag{ 17 }$$

where k = 1/2 and $B=2\vert {\Gamma}\vert \sqrt{D}/\left[\sqrt{\pi }\vert {\Omega}\vert \right]$ for a perfectly reactive target, while k = 1 and B = |Γ|κ/|Ω| for a partially reactive target.

3.2. Volume-averaged mean fFPT

We now turn to the many-particle case and obtain a simple approximation for the volume-averaged mean fFPT

$$\begin{equation}\overline{\langle {\mathcal{T}}_{N}\rangle }=\underset{0}{\overset{\infty }{\int }}\enspace \mathrm{d}t\enspace t\enspace \overline{{\rho }_{N}\left(t\vert {\mathbf{x}}_{1},\dots ,{\mathbf{x}}_{N}\right)}=\underset{0}{\overset{\infty }{\int }}\enspace \mathrm{d}t\enspace {\left[\overline{S\left(t\right)}\right]}^{N},\end{equation} \tag{ 18 }$$

where we used equation (6). We first split this integral into two parts,

$$\begin{equation}\overline{\langle {\mathcal{T}}_{N}\rangle }=\underset{0}{\overset{T}{\int }}\enspace \mathrm{d}t\enspace {\left[\overline{S\left(t\right)}\right]}^{N}+\underset{T}{\overset{\infty }{\int }}\enspace \mathrm{d}t\enspace {\left[\overline{S\left(t\right)}\right]}^{N}.\end{equation} \tag{ 19 }$$

The first integral will be evaluated explicitly by using the short-time asymptotic formula (14). Under an appropriate choice of T, the contribution of the second term is greatly attenuated for large N and can be ignored. The 'threshold time' T can be chosen to minimise the error of such an approximation but we will see that the final result does not depend on T, as expected.

3.2.1. Perfect reactivity

For a perfectly reactive target (κ = ∞) and a finite domain we use equation (14a) to get

$$\begin{equation}\overline{\langle {\mathcal{T}}_{N}\rangle }\simeq \frac{\pi \vert {\Omega}{\vert }^{2}}{2D\vert {\Gamma}{\vert }^{2}}\frac{1-\left(1+z\left(N+1\right)\right){\left(1-z\right)}^{N+1}}{\left(N+1\right)\left(N+2\right)},\end{equation} \tag{ 20 }$$

where $z=\frac{2\vert {\Gamma}\vert \sqrt{DT}}{\sqrt{\pi }\vert {\Omega}\vert }$ is a number between 0 and 1 that is set by our choice of T. When

$$\begin{equation}N\vert \mathrm{ln}\left(1-z\right)\vert \gg 1\end{equation} \tag{ 21 }$$

the correction term (1 − z)^N+1 can be neglected, and one gets

$$\begin{equation}\overline{\langle {\mathcal{T}}_{N}\rangle }\simeq \frac{\pi \vert {\Omega}{\vert }^{2}}{2D\vert {\Gamma}{\vert }^{2}}\frac{1}{\left(N+1\right)\left(N+2\right)}\quad \left(N\gg 1\right).\end{equation} \tag{ 22 }$$

This approximate expression indeed does not depend on T, as it should. Figure 2(a) illustrates the remarkable accuracy of this approximation in the case of an interval (0, L) with absorbing endpoints (here, |Γ|/|Ω| = 2/L). This high accuracy results from the fact that the asymptotic relation (14) is exponentially accurate for the interval,

$$\begin{equation}\overline{S\left(t\right)}\simeq 1-\frac{4\sqrt{Dt}}{L\sqrt{\pi }}+O\left({\mathrm{e}}^{-{L}^{2}/\left(4Dt\right)}\right)\quad \left(t\to 0\right),\end{equation} \tag{ 23 }$$

see the exact representation (C6b). For this reason, we kept the form 1/[(N + 1)(N + 2)] in equation (22) because its replacement by the leading term 1/N² would yield the correction O(N⁻³).

Zoom In Zoom Out Reset image size

In general, however, the next-order term in the short-time behaviour of the survival probability is of order O(t) that deteriorates the exponential accuracy of equation (22). In appendix B, we compute the correction O(N⁻³) term to this formula:

$$\begin{equation}\overline{\langle {\mathcal{T}}_{N}\rangle }\simeq \frac{\pi \vert {\Omega}{\vert }^{2}}{2D\vert {\Gamma}{\vert }^{2}}\frac{1}{{N}^{2}}\left(1-\left(1-\frac{\pi \vert {\Omega}\vert }{2\vert {\Gamma}\vert \mathcal{R}}\right)\frac{3}{N}+O\left({N}^{-2}\right)\right),\end{equation} \tag{ 24 }$$

where $\mathcal{R}$ is the mean curvature radius of the target defined in equation (B2) (note that $\mathcal{R}$ can be negative, see below). The leading, 1/N²-term was first mentioned in the seminal paper [79] and later re-discovered in [86, 87]. Moreover, a transition from an intermediate 1/N scaling to the ultimate 1/N² behaviour was analysed in [86]. A mathematical derivation of the leading 1/N²-term for a rather general class of diffusion processes was recently reported in [88]. Figure 2(b) shows the relative error of our approximation for diffusion inside a ball of radius L toward its absorbing boundary (i.e., Γ = ∂Ω and $\mathcal{R}=L$). As expected, the accuracy is lower at small N but still remains excellent. Finally, figure 2(c) presents the relative error for the case of a perfectly absorbing spherical target of radius R = 0.1L surrounded by a reflecting sphere of radius L (here, the mean curvature radius is negative: $\mathcal{R}=-R$). In contrast to previous examples, the asymptotic regime is established at much larger N, in agreement with the condition (21). In fact, here z ∝ |Γ| ∝ (R/L)² is 100 times smaller than in the case of an absorbing sphere, i.e., one needs a hundred-fold increase of N to get a comparable accuracy.

How can one rationalise the unexpected 1/N² scaling? As the starting points of the searchers are spread in the domain, the closest particle to the target has more chances to reach it first. As discussed in appendix D, the average distance $\overline{\delta }$ between the closest particle and the target is of the order of 1/N (and not 1/N^1/d as one might expect from a standard estimate of the average inter-particle distance in d dimensions). Moreover, if N ≫ 1, there are many particles that start from a comparable distance. Even if the closest particle will diffuse far away from the target, there is a high chance that one of the other particle started at the distance of order 1/N will hit the target. In other words, this problem resembles the diffusion of a single particle on the interval $\left(0,\overline{\delta }\right)$ with absorption at 0 (on the boundary) and reflection on $\overline{\delta }$, for which the mean FPT is ${\overline{\delta }}^{2}/\left(2D\right)\propto 1/{N}^{2}$. We conclude that the scaling 1/N² arises from the fact that many particle can start close to the target. Moreover, the generic properties of the heat kernel ensure that this scaling still holds for non-uniform initial distributions which do not exclude particles from a close vicinity of the target (see appendix D). As we will discuss in section 4, the decay $\overline{\langle {\mathcal{T}}_{N}\rangle }\propto {N}^{-2}$ is much faster and drastically different from the scaling law $\langle {\mathcal{T}}_{N}\rangle \propto 1/\mathrm{ln}\enspace N$ obtained earlier in the case of a fixed starting point [79].

3.2.2. Partial reactivity

For a partially reactive target (κ < ∞), equation (14b) entails a different scaling for the volume-averaged mean fFPT with N,

$$\begin{equation}\overline{\langle {\mathcal{T}}_{N}\rangle }\simeq \frac{\vert {\Omega}\vert }{\kappa \vert {\Gamma}\vert }{N}^{-1}\left(1+{\left(\frac{\kappa \vert {\Omega}\vert }{D\vert {\Gamma}\vert }\right)}^{\frac{1}{2}}{N}^{-\frac{1}{2}}+O\left({N}^{-1}\right)\right),\end{equation} \tag{ 25 }$$

i.e. the contribution to the volume-averaged mean fFPT in the regime of a reaction control vanishes only as a first inverse power of N. The derivation of the correction term is presented in appendix B. The leading 1/N-term was recently obtained in [88].

Figure 3 illustrates the quality of this approximation for several values of κ for an interval of length L with a partially reactive endpoint (in which case Γ = {0} and thus |Ω|/|Γ| = L). The accuracy is lower than in the former case of absorbing boundaries, because the error of the approximation, O(N⁻¹), decreases slower than that of equation (24). In addition, we plot the relative error in the case of a spherical partially reactive target of radius R = 0.1L surrounded by a reflecting sphere of radius L. In this setting, the coefficient in front of the correction term, $\sqrt{\kappa \vert {\Omega}\vert /\left(D\vert {\Gamma}\vert \right)}\approx 5.8$, is not small so that this correction term turns out to deteriorate the quality of the approximation at small N, as can be seen from the dashed line. In turn, when N is large, the correction term improves the accuracy, as expected.

Zoom In Zoom Out Reset image size

Since the ratio N/|Ω| in (25) can be interpreted as a mean concentration [A] of diffusing particles, our result is consistent with standard chemical kinetics, in which the reaction rate scales linearly with [A]. Note, however, that we also found the next-order term, which behaves as $\sqrt{\left[A\right]}$. In contrast, expression (24) vanishes as 1/N² and hence, is proportional to a squared concentration of particles. This is very different from the Smoluchowski rate which is linear with the concentration, because it appears as a limiting form when both the volume and the number of particles grow to infinity, while their ratio is kept fixed and is sufficiently small (see appendix A). It means that, in contrast to a standard Collins–Kimball relation, in which both rate controlling factors enter with the same power of concentration, for bounded domains in the limit N → ∞ the reaction inevitably enters into a regime of kinetic control, while the controlling factor of diffusion becomes subdominant and defines only some finite corrections. In a way, this situation is similar to a transition to kinetic control upon reduction of the size of the escape window, which was predicted recently for the so-called narrow escape problem [58, 109].

3.3. Variance of the fFPT

As the fFPT ${\mathcal{T}}_{N}$ includes two sources of randomness (from Brownian trajectories and from random initial positions), there are at least two ways of characterising the fluctuations of the ${\mathcal{T}}_{N}$. In the first way, one may use the variance defined as

$$\begin{equation}\text{Var}\left\{{\mathcal{T}}_{N}\right\}=\overline{\langle {\mathcal{T}}_{N}^{2}\rangle }-{\left[\overline{\langle {\mathcal{T}}_{N}\rangle }\right]}^{2},\end{equation} \tag{ 26 }$$

where the second term is just the square of the mean volume-averaged fFPT. From a physical point of view, it is more convenient to consider the volume-averaged conditional variance

$$\begin{equation}\overline{\mathrm{V}\mathrm{a}\mathrm{r}\left\{{\mathcal{T}}_{N}\right\}}=\overline{\left(\langle {\mathcal{T}}_{N}^{2}\rangle -\langle {\mathcal{T}}_{N}{\rangle }^{2}\right)}.\end{equation} \tag{ 27 }$$

In other words, one first evaluates the variance of ${\mathcal{T}}_{N}$ for a fixed set of initial positions x₁, ..., x_N, and then averages this conditional variance over these positions distributed uniformly in Ω.

To find both variances, we first evaluate the second moment of ${\mathcal{T}}_{N}$ with respect to the probability density ρ_N(t|x₁, ..., x_N),

$$\begin{align}\hfill \langle {\mathcal{T}}_{N}^{2}\rangle & =\underset{0}{\overset{\infty }{\int }}\enspace \mathrm{d}t\enspace {t}^{2}\enspace {\rho }_{N}\left(t\vert {\mathbf{x}}_{1},\dots ,{\mathbf{x}}_{N}\right)\hfill \\ \hfill & =2\underset{0}{\overset{\infty }{\int }}\enspace \mathrm{d}t\enspace t\enspace {S}_{N}\left(t\vert {\mathbf{x}}_{1},\dots ,{\mathbf{x}}_{N}\right),\hfill \end{align} \tag{ 28 }$$

and then average it over the starting points to find

$$\begin{equation}\overline{\langle {\mathcal{T}}_{N}^{2}\rangle }=2\underset{0}{\overset{\infty }{\int }}\enspace \mathrm{d}t\enspace t\enspace {\left[\overline{S\left(t\right)}\right]}^{N}.\end{equation} \tag{ 29 }$$

To get the variance in (26), it is sufficient to subtract the squared mean $\overline{\langle {\mathcal{T}}_{N}\rangle }$ studied in section 3.2. In turn, for the conditional variance in (27), one needs to subtract the squared first moment $\langle {\mathcal{T}}_{N}\rangle$, averaged over the starting points,

$$\begin{align}\hfill \overline{{\langle {\mathcal{T}}_{N}\rangle }^{2}}& =\underset{0}{\overset{\infty }{\int }}\enspace \mathrm{d}{t}_{1}\underset{0}{\overset{\infty }{\int }}\enspace \mathrm{d}{t}_{2}\overline{S\left({t}_{1}\vert {\mathbf{x}}_{1},\dots ,{\mathbf{x}}_{N}\right)S\left({t}_{2}\vert {\mathbf{x}}_{1},\dots ,{\mathbf{x}}_{N}\right)}\hfill \\ \hfill & =\underset{0}{\overset{\infty }{\int }}\enspace \mathrm{d}{t}_{1}\underset{0}{\overset{\infty }{\int }}\enspace \mathrm{d}{t}_{2}{\left[\overline{S\left({t}_{1}\right)S\left({t}_{2}\right)}\right]}^{N}.\hfill \end{align} \tag{ 30 }$$

In order to further simplify this expression we use the spectral expansions (8) and (11) along with the orthonormality of Laplacian eigenfunctions to derive the identities

$$\begin{equation*}\overline{S\left({t}_{1}\right)S\left({t}_{2}\right)}=\sum _{n}{c}_{n}\enspace {\mathrm{e}}^{-D{\lambda }_{n}\left({t}_{1}+{t}_{2}\right)}\end{equation*}$$

$$\begin{equation}=\overline{S\left({t}_{1}+{t}_{2}\right)},\end{equation} \tag{ 31a }$$

$$\begin{equation*}\overline{\rho \left({t}_{1}\right)S\left({t}_{2}\right)}=\sum _{n}D{\lambda }_{n}{c}_{n}\enspace {\mathrm{e}}^{-D{\lambda }_{n}\left({t}_{1}+{t}_{2}\right)}\end{equation*}$$

$$\begin{equation}=\overline{\rho \left({t}_{1}+{t}_{2}\right)},\end{equation} \tag{ 31b }$$

$$\begin{equation*}\overline{\rho \left({t}_{1}\right)\rho \left({t}_{2}\right)}=\sum _{n}{D}^{2}{\lambda }_{n}^{2}{c}_{n}\enspace {\mathrm{e}}^{-D{\lambda }_{n}\left({t}_{1}+{t}_{2}\right)}\end{equation*}$$

$$\begin{equation}={\left(-{\partial }_{t}\overline{\rho \left(t\right)}\right)}_{t={t}_{1}+{t}_{2}}.\end{equation} \tag{ 31c }$$

In particular, the first identity implies

$$\begin{equation}\overline{{\langle {\mathcal{T}}_{1}\rangle }^{2}}=\sum _{n}\frac{{c}_{n}}{{\left(D{\lambda }_{n}\right)}^{2}},\end{equation} \tag{ 32 }$$

which is twice smaller than $\overline{\langle {\mathcal{T}}_{1}^{2}\rangle }$, compare equation (13).

Using the identity (31a) we complete the above computation,

$$\begin{equation}\overline{{\langle {\mathcal{T}}_{N}\rangle }^{2}}=\underset{0}{\overset{\infty }{\int }}\enspace \mathrm{d}{t}_{1}\underset{0}{\overset{\infty }{\int }}\mathrm{d}{t}_{2}{\left[\overline{S\left({t}_{1}+{t}_{2}\right)}\right]}^{N},\end{equation} \tag{ 33 }$$

from which the volume-averaged conditional variance of the fFPT follows in the form

$$\begin{equation}\overline{\mathrm{V}\mathrm{a}\mathrm{r}\left\{{\mathcal{T}}_{N}\right\}}=2\underset{0}{\overset{\infty }{\int }}\enspace \mathrm{d}t\enspace t\enspace {\left[\overline{S\left(t\right)}\right]}^{N}-\underset{0}{\overset{\infty }{\int }}\enspace \mathrm{d}{t}_{1}\underset{0}{\overset{\infty }{\int }}\enspace \mathrm{d}{t}_{2}{\left[\overline{S\left({t}_{1}+{t}_{2}\right)}\right]}^{N}.\end{equation} \tag{ 34 }$$

To proceed, we substitute again $\overline{S\left(t\right)}$ by its short-time approximation (14). For a perfectly absorbing boundary (κ = ∞) we get from equation (29),

$$\begin{align}\hfill \overline{\langle {\mathcal{T}}_{N}^{2}\rangle }& \simeq 2\underset{0}{\overset{T}{\int }}\enspace \mathrm{d}t\enspace t\enspace {\left(1-B\sqrt{t}\right)}^{N}=\frac{4}{{B}^{4}}\underset{1-B\sqrt{T}}{\overset{1}{\int }}\mathrm{d}x{\left(1-x\right)}^{3}{x}^{N}\hfill \\ \hfill & \simeq \frac{24}{{B}^{4}\left(N+1\right)\left(N+2\right)\left(N+3\right)\left(N+4\right)},\hfill \end{align} \tag{ 35 }$$

where T is the cut-off time and $B=2\vert {\Gamma}\vert \sqrt{D}/\left[\sqrt{\pi }\vert {\Omega}\vert \right]$ (note that exponentially small corrections due to ${\left(1-B\sqrt{T}\right)}^{N}$ are neglected in equation (35)). Substituting equations (35) and (24) into equation (26), we get

$$\begin{equation}\text{Var}\left\{{\mathcal{T}}_{N}\right\}\simeq \frac{20}{{B}^{4}}{N}^{-4}.\end{equation} \tag{ 36 }$$

Similarly, we get

$$\begin{align}\hfill \overline{{\langle {\mathcal{T}}_{N}\rangle }^{2}}& \simeq \underset{0}{\overset{T}{\int }}\mathrm{d}{t}_{1}\underset{0}{\overset{T}{\int }}\mathrm{d}{t}_{2}{\left(1-B\sqrt{{t}_{1}+{t}_{2}}\right)}^{N}\hfill \\ \hfill & =\frac{2\sqrt{2}}{{B}^{4}}\left(\underset{1-B\sqrt{T}}{\overset{1}{\int }}\mathrm{d}x{\left(1-x\right)}^{3}{x}^{N}+\underset{1-\sqrt{2}B\sqrt{T}}{\overset{1-B\sqrt{T}}{\int }}\mathrm{d}x\left(1-x\right)\left(2{B}^{2}T-{\left(1-x\right)}^{2}\right){x}^{N}\right)\hfill \\ \hfill & \simeq \frac{12\sqrt{2}}{{B}^{4}\left(N+1\right)\left(N+2\right)\left(N+3\right)\left(N+4\right)}.\hfill \end{align} \tag{ 37 }$$

Combining these results, we get the volume-averaged conditional variance of the fFPT in the leading in the limit N → ∞ order:

$$\begin{equation}\overline{\mathrm{V}\mathrm{a}\mathrm{r}\left\{{\mathcal{T}}_{N}\right\}}\simeq \frac{12\left(2-\sqrt{2}\right)}{{B}^{4}}\enspace {N}^{-4}.\end{equation} \tag{ 38 }$$

Dividing this variance by the squared volume-averaged mean value of the fFPT produces

$$\begin{equation}\frac{\overline{\mathrm{V}\mathrm{a}\mathrm{r}\left\{{\mathcal{T}}_{N}\right\}}}{{\left[\overline{\langle {\mathcal{T}}_{N}\rangle }\right]}^{2}}\simeq \frac{3\left(2-\sqrt{2}\right)}{4}\approx 0.44\quad \left(N\gg 1\right).\end{equation} \tag{ 39 }$$

Similarly, one can also obtain the squared coefficient of variation of ${\langle {\mathcal{T}}_{N}\rangle }_{{\mathbf{x}}_{1},\dots ,{\mathbf{x}}_{N}}$ with respect to the random variables x₁, ..., x_N,

$$\begin{equation}\frac{\overline{\langle {\mathcal{T}}_{N}{\rangle }^{2}}}{{\left[\overline{\langle {T}_{N}\rangle }\right]}^{2}}-1\simeq 3\sqrt{2}-1\approx 3.24\quad \left(N\gg 1\right).\end{equation} \tag{ 40 }$$

The following point is to be emphasised. In statistical analysis of the first-passage phenomena in bounded domains [68, 89, 90], the coefficient of variation of the PDF is a meaningful characteristic which probes its 'effective' broadness. Typically, in situations when this parameter is much less than unity, one deals with a narrow distribution and the actual behaviour is well-captured by its first moment—the mean FPT, which sets a unique time scale. Only in this case the FPTs are 'focused', i.e., concentrated around the mean value. Conversely, when the coefficient of variation is of order of unity or even exceeds it, the PDF, despite the fact that it possesses moments of arbitrary order, exhibits in some aspects a behaviour reminiscent of so-called broad distributions, like heavy-tailed distributions which do not possess all moments. In this case, fluctuations around the mean value are comparable to the mean value itself and hence, it is most likely that the values of FPTs observed in two realisations of the process will be disproportionately different. In the case at hand, we notice that the coefficient of variation of the volume-averaged fFPT, equation (39), is of order of unity, which implies that here the fluctuations of this quantity around its mean value are of order of this value itself. More striking, the fluctuations of the fFPT corresponding to randomly distributed initial positions exceed the mean value of the fFPT, see equation (40), which signifies that randomness in the initial positions of the searchers has a very pronounced effect on the spread of fFPTs. Therefore, we conclude that for both volume-averaged and 'bare' fFPTs no unique time scale exists and their actual behaviour cannot be fully characterised by their mean values.

3.4. Long-time behaviour

Now we consider the volume-averaged probability density for N particles. As usual, the long-time limit presents the simplest setting for bounded domains due to the exponential decay of the survival probability and the PDF, as well as their volume averages. In particular one gets from equation (9)

$$\begin{equation}\overline{{\rho }_{N}\left(t\right)}\simeq ND{\lambda }_{1}{c}_{1}^{N}\enspace \mathrm{exp}\left(-ND{\lambda }_{1}t\right)\quad \left(t\to \infty \right),\end{equation} \tag{ 41 }$$

that is, the characteristic decay time 1/(Dλ₁) for a single particle gets simply reduced by the factor N. This steeper exponential decay shifts the PDF to smaller times. The decay time is a natural timescale of the diffusive exploration of the whole domain. For a single particle, this timescale is usually close to the mean FPT. Surprisingly enough, this is not true in case when multiple particles are searching for perfect targets starting from distinct random locations. As evidenced by our result (24), here the mean fFPT scales as 1/N², i.e., it is much smaller than 1/(Dλ₁N). Therefore, the speedup of the reaction kinetics due to a deployment of N searchers starting at distinct random positions appears to be really striking, as compared to a much more modest logarithmic increase in the efficiency predicted in the case when all of them start from the same point.

It is also instructive to compare the exponential decay in (41) to the asymptotic Weibull density (17). For partially reactive targets, k = 1 and the functional form of (17) generally agrees with result (41), except for the coefficients. In contrast, for the case of perfectly reactive targets, k = 1/2, and (17) exhibits a stretched-exponential decay with time. The difference between these two asymptotic behaviours stems from the order of how the limits are taken: the Weibull density (17) was derived for a fixed t as N → ∞ [88], whereas expression (41) corresponds to the limit t → ∞ with fixed N.

3.5. Fluctuations between individual realisations

When the starting point x₀ is random, the survival probability $S\left(t\vert {\mathbf{x}}_{0}\right)$ and the PDF $\rho \left(t\vert {\mathbf{x}}_{0}\right)$ for each fixed t can be considered as random variables themselves. This circumstance has been emphasised in the insightful paper [122] studying a survival of an immobile target in presence of diffusive traps. Moments of the survival probability regarded as a random variable, and the difference between the typical and mean behaviour was analysed [123–125] for the problem of survival of a diffusive particle in the presence of immobilised traps. The volume-averaged PDF, $\overline{\rho \left(t\right)}$, is simply the mean of $\rho \left(t\vert {\mathbf{x}}_{0}\right)$. In order to characterise fluctuations it is instructive to compute the variance of $\rho \left(t\vert {\mathbf{x}}_{0}\right)$.

3.5.1. Single particle

We characterise fluctuations by the squared coefficient of variation

$$\begin{equation}{\gamma }_{1}\left(t\right)=\frac{\overline{{\rho }^{2}\left(t\right)}-{\left[\overline{\rho \left(t\right)}\right]}^{2}}{{\left[\overline{\rho \left(t\right)}\right]}^{2}},\end{equation} \tag{ 42 }$$

where the second moment can be evaluated using equation (31c). At long times, we get

$$\begin{equation}{\gamma }_{1}\left(t\right)\to \frac{1}{{c}_{1}}-1\quad \left(t\to \infty \right),\end{equation} \tag{ 43 }$$

where c₁ is defined by equation (10). For instance, for an interval (0, L) with absorbing endpoints, c₁ = 8/π² and thus γ(∞) = π²/8 − 1 ≈ 0.23 is of order of unity.

For short times, equation (15) implies $\overline{{\rho }^{2}\left(t\right)}\simeq \left(B/{2}^{7/2}\right){t}^{-3/2}$ with $B=2\sqrt{D}\vert {\Gamma}\vert /\left[\sqrt{\pi }\vert {\Omega}\vert \right]$, and thus

$$\begin{equation}{\gamma }_{1}\left(t\right)\simeq \frac{\sqrt{\pi }\vert {\Omega}\vert }{16\vert {\Gamma}\vert \sqrt{2Dt}}+O\left(1\right)\quad \left(t\to 0\right).\end{equation} \tag{ 44 }$$

In other words, the coefficient of variation diverges in this limit, meaning that the volume-averaged density is no longer representative. This is confirmed by figure 5(a) see below.

3.5.2. Multiple particles

This illustrative computation can be extended to the case of N particles. We get

$$\begin{align}\hfill \overline{{\rho }_{N}^{2}\left(t\right)}& =N\overline{{\rho }^{2}\left(t\right)}{\left(\overline{{S}^{2}\left(t\right)}\right)}^{N-1}+N\left(N-1\right){\left(\overline{\rho \left(t\right)S\left(t\right)}\right)}^{2}{\left(\overline{{S}^{2}\left(t\right)}\right)}^{N-2}.\hfill \end{align} \tag{ 45 }$$

Using the identities (31), we find that the squared coefficient of variation obeys

$$\begin{align}\hfill 1+{\gamma }_{N}\left(t\right)& ={\left(\frac{\overline{{S}^{2}\left(t\right)}}{{\left[\overline{S\left(t\right)}\right]}^{2}}\right)}^{N-2}\left\{\frac{1}{N}\frac{\overline{{\rho }^{2}\left(t\right)}}{{\left[\overline{\rho \left(t\right)}\right]}^{2}}\frac{\overline{{S}^{2}\left(t\right)}}{{\left[\overline{S\left(t\right)}\right]}^{2}}+\left(1-1/N\right){\left(\frac{\overline{\rho \left(t\right)S\left(t\right)}}{\overline{\rho \left(t\right)}\cdot \overline{S\left(t\right)}}\right)}^{2}\right\}.\hfill \end{align} \tag{ 46 }$$

At long times, the ratios entering this expression approach 1/c₁ so that we get

$$\begin{equation}{\gamma }_{N}\left(t\right)\simeq {c}_{1}^{-N}-1\quad \left(t\to \infty \right).\end{equation} \tag{ 47 }$$

The Hölder inequality, applied to the first (positive) eigenfunction u₁, implies that

$$\begin{equation}{c}_{1}{\leqslant}1.\end{equation} \tag{ 48 }$$

As a consequence, γ_N(t) exponentially diverges with N for large t. In section 3.6 below we explain this paradoxical blow-up of fluctuations.

At short times, we find

$$\begin{equation}{\gamma }_{N}\left(t\right)\simeq \frac{\sqrt{\pi }\vert {\Omega}\vert }{16\vert {\Gamma}\vert \sqrt{2Dt}}{N}^{-1}-\left(1+1/N\right)\frac{1}{2}\quad \left(t\to 0\right)\end{equation} \tag{ 49 }$$

(note that one can also get correction terms O(1)/N from the first expression). As the number of particles increases, there are two effects. On one hand, the first (divergent) term is progressively attenuated. In other words, for any fixed t, taking N large enough we diminish γ_N(t) and hence, get narrower distribution of the random variable ρ = ρ_N(t|x₁, ..., x_N) (here we forget that ρ_N is itself the probability density and consider it as a random variable due to the randomness of x₁, ..., x_N). On the other hand, the second term increases and rapidly reaches a constant contribution 1/2 to the squared coefficient of variation. Even though the second term comes with the negative sign, it cannot render the coefficient of variation negative: indeed, the short-time asymptotic relation (49) is only valid for t small enough such that the first term in (49) provides the dominant contribution.

Figure 4 shows an excellent agreement between the exact expression for γ_N(t), its short-time asymptotic behaviour and Monte Carlo simulations results for an interval (0, L). Note that we replaced the O(1) term in equation (44) by −1 because for an interval, the short-time asymptotic formulas (15) are exponentially accurate, and the only correction −1 comes from the definition of the coefficient of variation.

Zoom In Zoom Out Reset image size

Figure 5 illustrates the PDFs ρ_N(t|x₁, ..., x_N) for 100 random combinations of the starting points x₁, ..., x_N chosen uniformly and independently on the interval (0, L). Fluctuations of the PDFs ρ_N(t|x₁, ..., x_N) around their volume-averaged mean $\overline{{\rho }_{N}\left(t\right)}$ are present and significant for all N ranging from 1 to 1000. The broadness of these fluctuations can be even better seen on figure 6 which illustrates 10 000 random realisations of the starting points. This is a rather counter-intuitive observation as a sort of convergence to the volume-averaged PDF is expected as N increases. Note that the asymptotic Weibull density (17) accurately describes the volume-averaged mean $\overline{{\rho }_{N}\left(t\right)}$ already for N = 100.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We emphasise that the PDF is rapidly shifted towards smaller times as N increases. Intuitively, one could expect that this shift is controlled by the long-time exponential decay and its timescale, T₁/N, see section 3.4. However, the appropriate time scale here is the volume-averaged mean fFPT, $\overline{\langle {\mathcal{T}}_{N}\rangle }$, which decays much faster, as 1/N², see the vertical dashed lines. We stress again that the decay time here is not related to the mean fastest FPT.

The general shape of the PDFs ρ_N(t|x₁, ..., x_N) is significantly different from the volume-average $\overline{{\rho }_{N}\left(t\right)}$ in figure 5. Thus for a single realisation the starting point is a finite distance away from the target, which effects the exponential cutoff to very short times and the emergence of the peak (the most probable FPT). At longer times we see the exponential shoulder with time scale T₁/N. In figure 5(a) individual realisations feature a slight bend at intermediate times after the most probable FPT and the exponential shoulder, however, this feature is getting lost for increasing N. In this one-dimensional setting the pronounced plateau in the FPT density uncovered in [57] is not present, a further separation of the most probable time and the longest time scale T₁/N would build up for decreasing reactivity. The volume-averaged PDF $\overline{{\rho }_{N}\left(t\right)}$, in strong contrast, includes particles starting arbitrarily closely to the target, and here we do not see the initial exponential suppression. Instead $\overline{{\rho }_{N}\left(t\right)}$ decays monotonically, and the only remaining relevant scale remaining is $\overline{\langle {\mathcal{T}}_{N}\rangle }$. Nevertheless the FPT density remains broad, as quantified by the 10%-quantiles in figure 5 and our results for the coefficient of variation.

Note that the small q-quantiles can be approximated by using again the short-time asymptotic formula (14). In fact, writing

$$\begin{equation}1-q=\overline{{S}_{N}\left(t\right)}={\left[\overline{S\left(t\right)}\right]}^{N}\simeq {\left(1-B\sqrt{t}\right)}^{N},\end{equation} \tag{ 50 }$$

which is valid for small t and thus small q, one gets easily

$$\begin{equation}{t}_{q}\simeq {\left(\frac{1-{\left(1-q\right)}^{1/N}}{B}\right)}^{2}\simeq \frac{{q}^{2}}{{B}^{2}{N}^{2}}\quad \left(N\gg 1\right),\end{equation} \tag{ 51 }$$

where $B=2\sqrt{D}\vert {\Gamma}\vert /\left[\sqrt{\pi }\vert {\Omega}\vert \right]$. As can be seen in figure 5 the relative locations of the quantiles hardly change for N = 10, 100, to 1000 but are just shifted to shorter times, as expected. This approximate relation is indeed quite accurate. Here, the scaling is again as 1/N², which is a direct consequence of the leading term $O\left(\sqrt{t}\right)$ in the short-time expansion of the volume-averaged survival probability.

3.6. Paradoxical divergence and its rationalisation

The constant c₁ from equation (10) determines the behaviour of the coefficient of variation in the long-time limit. This constant is defined to be independent of the size of the domain, i.e., any dilation of the domain does not change c₁. For instance, one has c₁ = 8/π² ≈ 0.81 for an interval with absorbing endpoints (of any length) and c₁ = 6/π² ≈ 0.61 for an absorbing sphere (of any radius). According to inequality (48), this coefficient cannot exceed 1, and it is actually equal to 1 only when the eigenfunction u₁ is constant, in other words, for a reflecting boundary without any reaction. As a consequence, the coefficient of variation γ_N at long times exhibits an exponential divergence with N. This observation appears paradoxical. We here resolve this counter-intuitive behaviour.

First, we recall that assuming the limit N → ∞ while keeping the volume |Ω| fixed corresponds to a diverging concentration, which is evidently unphysical. Moreover, in typical applications, the diffusing particles search for a small target contained within a much larger confining domain with reflecting boundary. In this setting, the ground eigenfunction u₁(x) is close to a constant so that c₁ is close to 1. In other words, if the target is fixed but the confining domain grows, c₁ tends to 1, and the double limit N → ∞ and |Ω| → ∞ with a fixed concentration N/|Ω| leads to a finite value of γ_N and thus mends this paradoxical divergence. To clarify these issues, it is therefore important to investigate its behaviour in the small target limit. We first analyse an exactly solvable setting and then briefly discuss the general case.

We consider the case of a spherical target of radius R surrounded by a larger reflecting sphere of radius L. In appendix C.4, we provide the exact formula (C22) for the coefficients c_n, which depend on the solutions α_n of equation (C20). For small R the smallest solution α₀ of this equation is expected to be small. Denoting = R/L ≪ 1, we consider separately the cases of infinite and finite reactivity.

• (a)  
  For infinite reactivity (κ = ∞) one substitutes α₀ ≃ a₁ + a₂² + O(³) into equation (C20) and expands it into a Taylor series to determine the expansion coefficients a_i, from which we get

  $$\begin{equation}{c}_{1}=1-\frac{108}{175}{\left(R/L\right)}^{2}+O\left({{\epsilon}}^{3}\right).\end{equation} \tag{ 52 }$$

  As a consequence, the coefficient of variation at long times can be approximated as

  $$\begin{equation}{\gamma }_{N}\left(\infty \right)\simeq \mathrm{exp}\left(\pi \frac{144}{175}\left[A\right]{R}^{2}L\right)-1,\end{equation} \tag{ 53 }$$

  where [A] = N/|Ω| is the concentration of the diffusing particles.
• (b)  
  In case of a finite reactivity (κ < ∞) more terms are needed. Substituting α₀ ≃ a₁ + a₂² + a₃³ + a₄⁴ + O(⁵) into equation (C20), we find

  $$\begin{equation}{c}_{1}=1-\frac{108}{175}\frac{{\kappa }^{2}{R}^{4}}{{D}^{2}{L}^{2}}+O\left({{\epsilon}}^{5}\right),\end{equation} \tag{ 54 }$$

  and thus

  $$\begin{equation}{\gamma }_{N}\left(\infty \right)\simeq \mathrm{exp}\left(\pi \frac{144}{175}\frac{{\kappa }^{2}{R}^{4}L\left[A\right]}{{D}^{2}}\right)-1.\end{equation} \tag{ 55 }$$

  To ensure the validity of this asymptotic analysis, the first correction term in equation (54) has to be small, which is realised when the inequality κR² ≪ DL holds. This inequality is evidently valid in the low concentration limit when the coefficient of variation vanishes, or when diffusion is very fast. Moreover, it holds in case of a sufficiently low reactivity κ; in this case, evidently, before a reaction eventually takes place there are multiple encounters with a target interspersed with bulk excursions such that the initial spatial heterogeneity due to a random distribution particles gets speared away, effecting γ_N(∞) to be small. In contrast, a larger domain favours a higher degree of spatial heterogeneity and thus increases γ_N. Note also that, in principle, one can relax this restrictive relation between the parameters and use equations (C22) and (47) directly. This will permit us to compute c₁ and hence, γ_N(∞) for any values of the parameters.

In the above setting, the target was a small ball located in the bulk and surrounded by a large reflecting sphere. If the target is located on the outer sphere (e.g., a small circular hole), one has to consider mixed Dirichlet–Neumann boundary condition on the sphere, for which there is no exact solution for the survival probability. In this case, one can either resort to the approximate solution in [58], or rely on the asymptotic analysis in the narrow escape limit. In the latter case, Ward and Keller showed for a perfectly reactive target that

$$\begin{equation}{c}_{1}=1-{E}_{1}{\epsilon}+o\left({\epsilon}\right),\end{equation} \tag{ 56 }$$

where the coefficient E₁ was expressed as an integral of the first-order correction to the ground eigenfunction [121]. This coefficient can also be expressed in terms of the surface Neumann Green function. When the confining domain Ω is a ball, the explicit form of this Green function was given in [105]. As a consequence, one can access the coefficient E₁ and thus the asymptotic behaviour of c₁ in the narrow escape limit. Moreover, this computation is valid for multiple small targets located on the sphere. Importantly, the asymptotic relation (56) remains valid even for nonspherical three-dimensional domains, even though the computation of E₁ is much harder. We emphasise the distinct scaling of 1 − c₁ with : it is linear for a target on the boundary, and quadratic for a target in the bulk, see equations (52) and (56). In two dimensions, the approach of c₁ to 1 is much slower:

$$\begin{equation}{c}_{1}=1-\frac{{E}_{1}}{\mathrm{ln}\left(2/{\epsilon}\right)}+o\left(1/\mathrm{ln}\left(2/{\epsilon}\right)\right),\end{equation} \tag{ 57 }$$

where 2L is the length of the target region, see [104, 121] for details. An extension of these asymptotic results to partially reactive targets located on the boundary presents an interesting perspective.

3.7. Typical PDF

Following reference [122] we aim at determining the typical PDF of the fastest FPT

$$\begin{equation}{\rho }_{N,\mathrm{t}\mathrm{y}\mathrm{p}}\left(t\right)\equiv \frac{1}{\tau }\mathrm{exp}\left(\overline{\mathrm{ln}\left(\tau {\rho }_{N}\left(t\vert {\mathbf{x}}_{1},\dots ,{\mathbf{x}}_{N}\right)\right)}\right),\end{equation} \tag{ 58 }$$

where τ is a timescale to render the expression in the logarithm dimensionless. For this purpose, we need to compute

$$\begin{equation}\overline{\mathrm{ln}\left(\tau {\rho }_{N}\left(t\vert {\mathbf{x}}_{1},\dots ,{\mathbf{x}}_{N}\right)\right)}=U\left(t\right)+\sum _{n=1}^{N}\overline{\mathrm{ln}\left(S\left(t\vert {\mathbf{x}}_{n}\right)\right)},\end{equation} \tag{ 59 }$$

where

$$\begin{equation}U\left(t\right)\equiv \overline{\mathrm{ln}\left(\tau \rho \left(t\vert {\mathbf{x}}_{1}\right)/S\left(t\vert {\mathbf{x}}_{1}\right)+\cdots +\tau \rho \left(t\vert {\mathbf{x}}_{N}\right)/S\left(t\vert {\mathbf{x}}_{N}\right)\right)}.\end{equation} \tag{ 60 }$$

As x₁, ..., x_N are independent uniformly distributed starting points the sum in equation (59) contains N identical terms. In turn, the first term U(t) is more sophisticated as the logarithm couples different particles. Using the following representation for the logarithm

$$\begin{equation}\mathrm{ln}\left(p\right)=\underset{0}{\overset{\infty }{\int }}\frac{\mathrm{d}z}{z}\left({\mathrm{e}}^{-z}-{\mathrm{e}}^{-pz}\right),\end{equation} \tag{ 61 }$$

we get

$$\begin{equation}U\left(t\right)=\underset{0}{\overset{\infty }{\int }}\frac{\mathrm{d}z}{z}\left({\mathrm{e}}^{-z}-{\left[\eta \left(z;t\right)\right]}^{N}\right),\end{equation} \tag{ 62 }$$

where

$$\begin{equation}\eta \left(z;t\right)=\frac{1}{\vert {\Omega}\vert }\underset{{\Omega}}{\int }\mathrm{d}\mathbf{x}\enspace \mathrm{exp}\left(-z\tau \rho \left(t\vert \mathbf{x}\right)/S\left(t\vert \mathbf{x}\right)\right).\end{equation} \tag{ 63 }$$

3.7.1. Long-time behaviour

The long-time behaviour is easy to determine. As ρ(t|x) = −∂S(t|x)/∂t the ratio ρ(t|x)/S(t|x) is a positive function that is monotonously increasing from 0 at t = 0 to Dλ₁ as t → ∞. As a consequence equation (62) implies

$$\begin{equation}U\left(t\right)\to \mathrm{ln}\left(\tau ND{\lambda }_{1}\right)\quad \left(t\to \infty \right).\end{equation} \tag{ 64 }$$

Moreover, we have

$$\begin{equation}\overline{\mathrm{ln}\left(S\left(t\right)\right)}\simeq -Dt{\lambda }_{1}+{b}_{1}\quad \left(t\to \infty \right),\end{equation} \tag{ 65 }$$

where

$$\begin{align}\hfill {b}_{1}& =\mathrm{ln}\left(\vert {\Omega}{\vert }^{-1/2}\underset{{\Omega}}{\int }\mathrm{d}\mathbf{x}\enspace {u}_{1}\left(\mathbf{x}\right)\right)+\frac{1}{\vert {\Omega}\vert }\underset{{\Omega}}{\int }\mathrm{d}\mathbf{x}\enspace \mathrm{ln}\left(\vert {\Omega}{\vert }^{1/2}{u}_{1}\left(\mathbf{x}\right)\right),\hfill \end{align} \tag{ 66 }$$

and we recall that u₁(x) can be defined to be positive (factors |Ω|^±1/2 were included to get dimensionless quantities in logarithms). In this limit, one also has $\rho \left(t\vert {\mathbf{x}}_{0}\right)/S\left(t\vert {\mathbf{x}}_{0}\right)\simeq D{\lambda }_{1}$ so that

$$\begin{equation}{\rho }_{N,\mathrm{t}\mathrm{y}\mathrm{p}}\left(t\right)\simeq {\mathrm{e}}^{{b}_{1}N}ND{\lambda }_{1}\enspace {\mathrm{e}}^{-Dt{\lambda }_{1}N}\quad \left(t\to \infty \right).\end{equation} \tag{ 67 }$$

In contrast, the analysis of the short-time behaviour of the typical PDF is much more difficult and remains an open problem for future research.

We studied the fFPT distribution for diffusion of N Brownian particles in a finite domain Ω with smooth boundary ∂Ω. The particles react with a surface region Γ of perfect or partial reactivity contained in ∂Ω. We obtained the fFPT moments and probability density in the case of a uniform initial condition. In the previously considered case of a common, fixed initial position, the mean fFPT was shown to exhibit the very slow, 1/ln N scaling [74, 76, 79, 83–85]. In strong contrast, when the N particles are released random-uniformly in the domain Ω, the mean fFPT shows the much faster decay 1/N² for a perfectly reactive target and 1/N for a target with finite reactivity. In the former case, this scaling is also different from that of the decay time, T₁/N. We provided a rationale of this significantly altered behaviour in the dependence on the specific initial condition by scaling arguments and supported it by direct numerical computations.

While the scenario of a fixed initial condition bears relevance for the example of the released sperm cells considered in [74–76], spread-out initial conditions are relevant in other biological systems. Thus, in biofilms autoinducer molecules are released on the cell surfaces all over the colony [72, 73]. In intracellular gene regulation, so-called global regulatory proteins are distributed approximately evenly throughout the cell volume [43]. Other regulatory proteins may abound around their encoding gene [41], and genes that are controlled by a specific transcription factor tend to locate next to the gene encoding this transcription factor [42, 43]. These examples demonstrate the need for a theory that considers spread-out initial conditions as considered here.

In many biological systems the number of searching entities such as sperm cells or regulatory proteins are produced at the expense of chemical and energy resources. At the same time a large number of searchers guarantees a high degree of redundancy, which is of particular importance for the case of sperm cells in singular, crucial events such as reproduction. In contrast, for processes that are constantly running off, such as gene regulatory processes, the biochemical resources would quickly be drained if excessively many molecules had to be produced. In this sense, and given the above scenarios when the molecules abound in the vicinity of their target, the scaling 1/N and 1/N² provides a new perspective on why relatively few molecules may already be sufficient to effect comparatively speedy regulation.

While in our discussion we contrasted a fixed initial position with the scenario of random-uniform initial positions, figure 5 indicates that the value of the mean fFPT provides the correct scale even for individual random realisations of starting points. In this sense, one may roughly distinguish two classes of initial conditions: (i) either the particles are allowed to start near their target (distributed in the entire domain or concentrated in a subdomain around the target), (ii) or the initial particle position is characterised by a minimal distance to the target. The former class corresponds to our analysis here (see also [88]), the latter class is described by the results in previous works [74, 76, 83–85].

In summary, we believe that our results for the statistic of the fFPT for spread initial particle positions provide a fresh perspective to the field of many-particle search for an immobile target in a finite domain and its applications. It should be of interest to further develop this approach both from a mathematical point of view and with regard to concrete applications.

DSG thanks MJ Ward for fruitful discussions and acknowledges a partial financial support from the Alexander von Humboldt Foundation through a Bessel Research Award. RM acknowledges funding from the German Science Foundation (DFG, Grant ME 1535/7-1) as well as support from the Foundation for Polish Science (Fundacja na rzecz Nauki Polskiej, FNP) within an Alexander von Humboldt Polish Honorary Research Scholarship. GO is grateful to B Meerson for helpful discussions. The authors also thank the unknown referee for bringing our attention to reference [88] that was published after the submission of our manuscript and provided complementary views on the fFPT. We acknowledge the support of German Research Foundation (DFG) and Open Access Publication Fund of Potsdam University.

In this appendix we recall the behaviour in the conventional thermodynamic limit when both N and |Ω| tend to infinity with the concentration [A] = N/|Ω| being fixed. In this limit, one gets

$$\begin{align*}\hfill \overline{{S}_{N}}& ={\left[\overline{S\left(t\right)}\right]}^{N}\simeq {\left(1-\frac{1}{\vert {\Omega}\vert }\underset{{\Omega}}{\int }\mathrm{d}{\mathbf{x}}_{0}\left(1-S\left(t\vert {\mathbf{x}}_{0}\right)\right)\right)}^{N}\hfill \\ \hfill & \simeq \mathrm{exp}\left(-\left[A\right]\underset{{\Omega}}{\int }\mathrm{d}{\mathbf{x}}_{0}\left(1-S\left(t\vert {\mathbf{x}}_{0}\right)\right)\right)=\mathrm{exp}\left(-\left[A\right]F\left(t\right)\right),\hfill \end{align*}$$

where

$$\begin{equation}F\left(t\right)=\underset{{\Omega}}{\int }\mathrm{d}{\mathbf{x}}_{0}\left(1-S\left(t\vert {\mathbf{x}}_{0}\right)\right).\end{equation} \tag{ A1 }$$

For the survival probability of a spherical target of radius R one gets (see, e.g., [120])

$$\begin{align}\hfill F\left(t\right)& =4\pi R\left[\frac{Dt}{1+D/\left(\kappa R\right)}+\frac{2\sqrt{Dt}\enspace R}{\sqrt{\pi }{\left(1+D/\left(\kappa R\right)\right)}^{2}}-\frac{RD/\kappa }{{\left(1+D/\left(\kappa R\right)\right)}^{3}}\left(1-\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}\mathrm{x}\left(\sqrt{Dt}\left(1/R+\kappa /D\right)\right)\right)\right].\hfill \end{align} \tag{ A2 }$$

In the limit t → ∞ the dominant behaviour of F(t) is provided by the first term, which yields the Collins–Kimball relation

$$\begin{equation}\frac{t}{F\left(t\right)}=\frac{1}{4\pi DR}+\frac{1}{K},\end{equation} \tag{ A3 }$$

where K = 4πR²κ is the intrinsic reaction constant.

In the limit of an infinitely large intrinsic reactivity (κ → ∞), one retrieves from equation (A2) the classical Smoluchowski result

$$\begin{equation}F\left(t\right)=4\pi RDt\left(1+\frac{2R}{\sqrt{\pi Dt}}\right).\end{equation} \tag{ A4 }$$

Our computation of the volume-averaged mean fFPT relied on the leading term of the short-time expansion of the survival probability. We can improve this computation by accounting for the next-order correction.

B.1. Perfect reactivity

For the perfect reactivity, one has [114, 115]

$$\begin{equation}S\left(t\right)\simeq 1-B\sqrt{t}+Ct+O\left({t}^{3/2}\right)\quad \left(t\to 0\right),\end{equation} \tag{ B1 }$$

where $B=\left(2/\sqrt{\pi }\right)\sqrt{D}\frac{\vert {\Gamma}\vert }{\vert {\Omega}\vert }$ and

$$\begin{equation}C=\frac{D\vert {\Gamma}\vert }{\vert {\Omega}\vert \mathcal{R}},\quad \frac{1}{\mathcal{R}}=\frac{1}{\vert {\Gamma}\vert }\underset{{\Gamma}}{\int }\mathrm{d}\mathbf{s}\enspace H\left(\mathbf{s}\right),\end{equation} \tag{ B2 }$$

where H(s) is the mean curvature of the boundary at the point s. Substituting this expression into the first integral in equation (19), we get

$$\begin{equation*}\overline{\langle {\mathcal{T}}_{N}\rangle }\simeq \underset{0}{\overset{T}{\int }}\enspace \mathrm{d}t\enspace {\left(1-B\sqrt{t}+Ct\right)}^{N}=\underset{{z}_{0}}{\overset{1}{\int }}\enspace \mathrm{d}z\enspace f\left(z\right)\enspace {z}^{N},\end{equation*}$$

where we changed the integration variable as $z=1-B\sqrt{t}+Ct$, with ${z}_{0}=1-B\sqrt{T}+CT$ and

$$\begin{equation*}f\left(z\right)=\frac{B-\sqrt{{B}^{2}-4{B}_{1}\left(1-z\right)}}{C\sqrt{{C}^{2}-4{B}_{1}\left(1-z\right)}}.\end{equation*}$$

For large N, the main contribution comes from the vicinity of z = 1 so that one can expand the function f(z) into a Taylor series around this point

$$\begin{equation*}f\left(z\right)={f}^{\prime }\left(1\right)\left(z-1\right)+\frac{1}{2}{f}^{{\prime\prime}}\left(1\right){\left(z-1\right)}^{2}+O\left({\left(z-1\right)}^{3}\right),\end{equation*}$$

and we used that f(1) = 0. Integrating term by term, one gets

$$\begin{equation*}\overline{\langle {\mathcal{T}}_{N}\rangle }\simeq \frac{-2{f}^{\prime }\left(1\right)}{\left(N+1\right)\left(N+2\right)}+\frac{2{f}^{{\prime\prime}}\left(1\right)}{\left(N+1\right)\left(N+2\right)\left(N+3\right)}+\dots ,\end{equation*}$$

with f'(1) = −1/B² and f''(1) = 6C/B⁴. As a consequence, we get

$$\begin{equation}\overline{\langle {\mathcal{T}}_{N}\rangle }\simeq \frac{2/{B}^{2}}{{N}^{2}}+\frac{12C/{B}^{4}-6/{B}^{2}}{{N}^{3}}+O\left({N}^{-4}\right),\end{equation} \tag{ B3 }$$

which takes the form of equation (24).

B.2. Partial reactivity

For partial reactivity, one has [116]

$$\begin{equation}S\left(t\right)\simeq 1-Bt+C{t}^{3/2}+O\left({t}^{2}\right)\quad \left(t\to 0\right),\end{equation} \tag{ B4 }$$

where

$$\begin{equation}B=\frac{\kappa \vert {\Gamma}\vert }{\vert {\Omega}\vert },\quad C=\frac{4{\left(\kappa /D\right)}^{2}\vert {\Gamma}\vert }{3\sqrt{\pi }\vert {\Omega}\vert }{D}^{3/2}.\end{equation} \tag{ B5 }$$

Substituting this expression into the first integral in equation (19), we get

$$\begin{equation*}\overline{\langle {\mathcal{T}}_{N}\rangle }\simeq \underset{0}{\overset{T}{\int }}\enspace \mathrm{d}t\enspace {\left(1-Bt+C{t}^{3/2}\right)}^{N}=\underset{{z}_{0}}{\overset{1}{\int }}\enspace \mathrm{d}z\enspace f\left(z\right)\enspace {z}^{N},\end{equation*}$$

where we changed the integration variable as z = 1 − Bt + Ct^3/2, with z₀ = 1 − BT + CT^3/2, $f\left(z\right)=1/\left(B-3C\sqrt{t\left(z\right)}/2\right)$, and t(z) is the inverse of the above function z(t). Even though f(z) can be written explicitly in terms of the root of the cubic polynomial, it is not needed as we only use its expansion around z = 1, which reads

$$\begin{equation*}f\left(z\right)=\frac{1}{B}+\frac{3C}{2{B}^{5/2}}{\left(1-z\right)}^{1/2}+O\left(1-z\right).\end{equation*}$$

Integrating term by term, one gets

$$\begin{equation}\overline{\langle {\mathcal{T}}_{N}\rangle }\simeq \frac{1}{B}{N}^{-1}+\frac{3\sqrt{\pi }C}{4{B}^{5/2}}{N}^{-3/2}+O\left({N}^{-2}\right),\end{equation} \tag{ B6 }$$

which takes the form of equation (25).

In this Appendix, we summarise the well-known explicit solutions for two simple settings that we used for illustrations: an interval and a sphere.

C.1. Interval with absorbing endpoints

For the interval (0, L) with absorbing endpoints at 0 and L, the survival probability admits two equivalent explicit expansions

$$\begin{equation}S\left(t\vert {x}_{0}\right)=1-\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}\left(\frac{{x}_{0}}{\sqrt{4Dt}}\right)+\sum _{k=1}^{\infty }{\left(-1\right)}^{k}\left[\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}\left(\frac{kL-{x}_{0}}{\sqrt{4Dt}}\right)-\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}\left(\frac{kL+{x}_{0}}{\sqrt{4Dt}}\right)\right]\end{equation} \tag{ C1a }$$

$$\begin{equation}=2\sum _{k=1}^{\infty }\mathrm{sin}\left(\pi k{x}_{0}/L\right)\frac{1-{\left(-1\right)}^{k}}{\pi k}{\mathrm{e}}^{-{\pi }^{2}{k}^{2}Dt/{L}^{2}}.\end{equation} \tag{ C1b }$$

The first relation is numerically more convenient at small t, whereas the second at large t. We also have

$$\begin{equation}\rho \left(t\vert {x}_{0}\right)=\frac{1}{\sqrt{4\pi D{t}^{3}}}\left[{x}_{0}\enspace {\mathrm{e}}^{-{x}_{0}^{2}/\left(4Dt\right)}+\sum _{k=1}^{\infty }{\left(-1\right)}^{k+1}\left(\left(kL-{x}_{0}\right){\mathrm{e}}^{-{\left(kL-{x}_{0}\right)}^{2}/\left(4Dt\right)}-\left(kL+{x}_{0}\right){\mathrm{e}}^{-{\left(kL+{x}_{0}\right)}^{2}/\left(4Dt\right)}\right)\right]\end{equation} \tag{ C2a }$$

$$\begin{equation}=\frac{2D}{{L}^{2}}\sum _{k=1}^{\infty }\mathrm{sin}\left(\pi k{x}_{0}/L\right)\left(1-{\left(-1\right)}^{k}\right)\left(\pi k\right){\mathrm{e}}^{-{\pi }^{2}{k}^{2}Dt/{L}^{2}}.\end{equation} \tag{ C2b }$$

Note also that the Laplace-transformed probability density reads

$$\begin{equation}\tilde {\rho }\left(p\vert {x}_{0}\right)=\frac{\mathrm{sinh}\left({x}_{0}\sqrt{p/D}\right)+\mathrm{sinh}\left(\left(L-{x}_{0}\right)\sqrt{p/D}\right)}{\mathrm{sinh}\left(L\sqrt{p/D}\right)},\end{equation} \tag{ C3 }$$

from which one can easily compute positive-order moments.

The volume-averaged probability density is then

$$\begin{equation}\overline{\rho \left(t\right)}=\frac{2\sqrt{D}}{L\sqrt{\pi t}}\left(1+2\sum _{k=1}^{\infty }{\left(-1\right)}^{k}\enspace {\mathrm{e}}^{-{\left(kL\right)}^{2}/\left(4Dt\right)}\right)\end{equation} \tag{ C4a }$$

$$\begin{equation}=\frac{8D}{{L}^{2}}\sum _{k=1}^{\infty }{\mathrm{e}}^{-{\pi }^{2}{\left(2k-1\right)}^{2}Dt/{L}^{2}},\end{equation} \tag{ C4b }$$

and its Laplace transform of this density reads

$$\begin{equation}\mathcal{L}\left\{\overline{\rho \left(t\right)}\right\}=\frac{2}{z}\left(\frac{2\enspace \mathrm{cosh}\enspace z}{\mathrm{sinh}\enspace z}-\frac{\mathrm{cosh}\left(z/2\right)}{\mathrm{sinh}\left(z/2\right)}\right),\end{equation} \tag{ C5 }$$

where $z=L\sqrt{p/D}$. One also has

$$\begin{equation}\overline{S\left(t\right)}=1-\frac{4\sqrt{Dt}}{L\sqrt{\pi }}\left(1+2\sum _{k=1}^{\infty }{\left(-1\right)}^{k}\left({\mathrm{e}}^{-{\left(kL\right)}^{2}/\left(4Dt\right)}-\frac{\sqrt{\pi }kL}{\sqrt{4Dt}}\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}\left(kL/\sqrt{4Dt}\right)\right)\right)\end{equation} \tag{ C6a }$$

$$\begin{equation}=\frac{8}{{\pi }^{2}}\sum _{k=1}^{\infty }\frac{{\mathrm{e}}^{-{\pi }^{2}{\left(2k-1\right)}^{2}Dt/{L}^{2}}}{{\left(2k-1\right)}^{2}}.\end{equation} \tag{ C6b }$$

For an interval (0, L), one has λ_n = π²n²/L², ${u}_{n}\left(x\right)=\sqrt{2/L}\enspace \mathrm{sin}\left(\pi nx/L\right)$, and ${c}_{n}=2{\left(1-{\left(-1\right)}^{n}\right)}^{2}/\left({\pi }^{2}{n}^{2}\right)$.

We also easily get

$$\begin{equation}\overline{{\delta }^{2}\left({\mathbf{x}}_{0}\right)}=\frac{{L}^{2}}{12}\end{equation} \tag{ C7 }$$

and

$$\begin{equation}{b}_{1}=\frac{1}{L}\underset{0}{\overset{L}{\int }}\enspace \mathrm{d}x\enspace \mathrm{ln}\left(\sqrt{2}\mathrm{sin}\left(\pi x/L\right)\right)=-\frac{1}{2}\mathrm{ln}\enspace 2.\end{equation} \tag{ C8 }$$

C.2. Interval with a partially reactive endpoint

We also consider the interval (0, L) with partially reactive endpoint at L and reflecting endpoint at 0, for which the Laplacian eigenvalues and eigenfunctions are (see, e.g., [126])

$$\begin{equation}{\lambda }_{n}={\alpha }_{n}^{2}/{L}^{2},\quad {u}_{n}\left(x\right)={\beta }_{n}\enspace \mathrm{cos}\left({\alpha }_{n}x/L\right),\end{equation} \tag{ C9 }$$

with

$$\begin{equation}{\beta }_{n}={\left(\frac{{h}^{2}+{\alpha }_{n}^{2}}{{h}^{2}+h+{\alpha }_{n}^{2}}\right)}^{1/2},\end{equation} \tag{ C10 }$$

and α_n are the positive zeros of the equation α_n sin α_n = h cos α_n, and h = κL/D. For each n = 1, 2, ..., α_n belongs to an interval [π(n − 1), π(n − 1/2)) that facilitates their numerical computation. Note that this problem is equivalent to diffusion on the interval (−L, L) whose both endpoints are partially reactive. The spectral expansion of the survival probability is [126]

$$\begin{align}\hfill S\left(t\vert {x}_{0}\right)& =2\sum _{n=1}^{\infty }\frac{{\left(-1\right)}^{n-1}h\sqrt{{h}^{2}+{\alpha }_{n}^{2}}}{{\alpha }_{n}\left({h}^{2}+h+{\alpha }_{n}^{2}\right)}{\mathrm{e}}^{-Dt{\alpha }_{n}^{2}/{L}^{2}}\enspace \mathrm{cos}\left({\alpha }_{n}{x}_{0}/L\right).\hfill \end{align} \tag{ C11 }$$

Its volume average reads

$$\begin{equation}\overline{S\left(t\right)}=2{h}^{2}\sum _{n=1}^{\infty }\frac{{\mathrm{e}}^{-Dt{\alpha }_{n}^{2}/{L}^{2}}}{{\alpha }_{n}^{2}\left({h}^{2}+h+{\alpha }_{n}^{2}\right)}.\end{equation} \tag{ C12 }$$

The probability density ρ(t|x₀) and its volume average $\overline{\rho \left(t\right)}$ follow by taking the time derivative.

C.3. Absorbing sphere

We also consider the first-passage problem to the absorbing boundary of a ball of radius L, for which the survival probability is also explicitly known,

$$\begin{equation}S\left(t\vert {\mathbf{x}}_{0}\right)=1-\frac{R}{\vert {\mathbf{x}}_{0}\vert }\sum _{k=1}^{\infty }\left[\mathrm{e}\mathrm{r}\mathrm{f}\left(\frac{\left(2k-1\right)L+\vert {\mathbf{x}}_{0}\vert }{\sqrt{4Dt}}\right)-\mathrm{e}\mathrm{r}\mathrm{f}\left(\frac{\left(2k-1\right)L-\vert {\mathbf{x}}_{0}\vert }{\sqrt{4Dt}}\right)\right]\end{equation} \tag{ C13a }$$

$$\begin{equation}=2\sum _{n=1}^{\infty }{\left(-1\right)}^{n+1}\frac{\mathrm{sin}\left(\pi n\vert {\mathbf{x}}_{0}\vert /L\right)}{\pi n\vert {\mathbf{x}}_{0}\vert /L}{\mathrm{e}}^{-{\pi }^{2}{n}^{2}Dt/{L}^{2}}.\end{equation} \tag{ C13b }$$

As a consequence, the PDF reads

$$\begin{align}\hfill \rho \left(t\vert {\mathbf{x}}_{0}\right)& =\frac{L}{\vert {\mathbf{x}}_{0}\vert \sqrt{4\pi D{t}^{3}}}\sum _{k=1}^{\infty }\left(\left(\left(2k-1\right)L-\vert {\mathbf{x}}_{0}\vert \right){\mathrm{e}}^{-{\left(\left(2k-1\right)L-\vert {\mathbf{x}}_{0}\vert \right)}^{2}/\left(4Dt\right)}\right.\hfill \\ \hfill & \quad \left.-\left(\left(2k-1\right)L+\vert {\mathbf{x}}_{0}\vert \right){\mathrm{e}}^{-{\left(\left(2k-1\right)L+\vert {\mathbf{x}}_{0}\vert \right)}^{2}/\left(4Dt\right)}\right)\hfill \\ \hfill & =\frac{2D}{{L}^{2}}\sum _{n=1}^{\infty }{\left(-1\right)}^{n+1}\frac{\mathrm{sin}\left(\pi n\vert {\mathbf{x}}_{0}\vert /L\right)}{\pi n\vert {\mathbf{x}}_{0}\vert /L}{\left(\pi n\right)}^{2}\enspace {\mathrm{e}}^{-{\pi }^{2}{n}^{2}Dt/{L}^{2}}.\hfill \end{align} \tag{ C14 }$$

We also get

$$\begin{equation}\overline{S\left(t\right)}=1+\frac{3Dt}{{L}^{2}}-\frac{6\sqrt{Dt}}{L\sqrt{\pi }}\left\{1+2\sum _{k=1}^{\infty }\left({\mathrm{e}}^{-{k}^{2}{L}^{2}/\left(Dt\right)}-\sqrt{\pi }\left(kL/\sqrt{Dt}\right)\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}\left(kL/\sqrt{Dt}\right)\right)\right\}\end{equation} \tag{ C15a }$$

$$\begin{equation}=\frac{6}{{\pi }^{2}}\sum _{n=1}^{\infty }\frac{{\mathrm{e}}^{-{\pi }^{2}{n}^{2}Dt/{L}^{2}}}{{n}^{2}}\end{equation} \tag{ C15b }$$

and

$$\begin{equation}\overline{\rho \left(t\right)}=-\frac{3D}{{L}^{2}}+\frac{3\sqrt{D}}{L\sqrt{\pi t}}+\frac{6\sqrt{D}}{L\sqrt{\pi t}}\sum _{k=1}^{\infty }{\mathrm{e}}^{-{k}^{2}{L}^{2}/\left(Dt\right)}\end{equation} \tag{ C16a }$$

$$\begin{equation}=\frac{6D}{{L}^{2}}\sum _{n=1}^{\infty }{\mathrm{e}}^{-{\pi }^{2}{n}^{2}Dt/{L}^{2}}.\end{equation} \tag{ C16b }$$

We also easily get

$$\begin{equation}\overline{{\delta }^{2}\left({\mathbf{x}}_{0}\right)}=\frac{4\pi }{4\pi {L}^{3}/3}\underset{0}{\overset{L}{\int }}\mathrm{d}r\enspace {r}^{2}\enspace {\left(L-r\right)}^{2}=\frac{{L}^{2}}{10}\end{equation} \tag{ C17 }$$

and

$$\begin{align}\hfill {b}_{1}& =\mathrm{ln}\left(\frac{1}{\sqrt{4\pi {L}^{3}/3}}\underset{{\Omega}}{\int }\mathrm{d}\mathbf{x}{\beta }_{1}\frac{\mathrm{sin}\left(\pi \vert \mathbf{x}\vert /L\right)}{\pi \vert \mathbf{x}\vert /L}\right)+\frac{1}{4\pi {L}^{3}/3}\underset{{\Omega}}{\int }\mathrm{d}\mathbf{x}\enspace \mathrm{ln}\left(\sqrt{4\pi {L}^{3}/3}{\beta }_{1}\frac{\mathrm{sin}\left(\pi \vert \mathbf{x}\vert /L\right)}{\pi \vert \mathbf{x}\vert /L}\right)\hfill \\ \hfill & =-\frac{3\zeta \left(3\right)}{2{\pi }^{2}}-2\enspace \mathrm{ln}\enspace \pi ,\hfill \end{align} \tag{ C18 }$$

where ζ(z) is the Riemann zeta function with ζ(3) ≈ 1.2021. Here we used the following expression for the first eigenfunction: ${u}_{1}\left(\mathbf{x}\right)={\beta }_{1}\frac{\mathrm{sin}\left(\pi \vert \mathbf{x}\vert /L\right)}{\pi \vert \mathbf{x}\vert /L}$, where ${\beta }_{1}=\sqrt{\pi /\left(2{L}^{3}\right)}$ is the normalisation factor. Note also that c₁ = 6/π².

C.4. Spherical target surrounded by a reflecting sphere

We also consider a more elaborate situation of a spherical target of radius R, which is surrounded by a larger reflecting concentric sphere of radius L. The related FPT distribution for a single particle was studied in [57]. In particular, the exact spectral expansion for the survival probability was provided:

$$\begin{equation}S\left(t\vert {\mathbf{x}}_{0}\right)=\sum _{n=1}^{\infty }{u}_{n}\left(\vert {\mathbf{x}}_{0}\vert \right)\enspace {\mathrm{e}}^{-Dt{\lambda }_{n}},\end{equation} \tag{ C19 }$$

where

$$\begin{align*}\hfill {\lambda }_{n}& ={\alpha }_{n}^{2}/{\left(L-R\right)}^{2},\hfill \\ \hfill {u}_{n}\left(r\right)& =\frac{{b}_{n}}{{\lambda }_{n}{\left(L-R\right)}^{2}}\frac{L{\alpha }_{n}\enspace \mathrm{cos}\left({\alpha }_{n}\frac{L-r}{L-R}\right)-\left(L-R\right)\mathrm{sin}\left({\alpha }_{n}\frac{L-r}{L-R}\right)}{r{\alpha }_{n}},\hfill \\ \hfill {b}_{n}& =\frac{2{R}^{2}{\alpha }_{n}^{2}}{\left(RL+\mu \left({L}^{2}+{R}^{2}\right)\right){\alpha }_{n}\enspace \mathrm{sin}\enspace {\alpha }_{n}+R\left(\mu L{\alpha }_{n}^{2}-R\right)\mathrm{cos}\enspace {\alpha }_{n}},\hfill \end{align*}$$

μ = D/(κ(L − R)), and α_n (with n = 1, 2, .... are positive solutions of the equation

$$\begin{equation}\mathrm{tan}\enspace {\alpha }_{n}=\frac{{\alpha }_{n}\left(RL+\mu {\left(L-R\right)}^{2}\right)}{R\left(L-R\right)+\mu {\left(L-R\right)}^{2}+\mu LR{\alpha }_{n}^{2}}.\end{equation} \tag{ C20 }$$

One easily gets $H\left(t\vert {\mathbf{x}}_{0}\right)$ by taking the time derivative.

The volume-averaged survival probability is obtained by integration over x₀,

$$\begin{align}\hfill \overline{S\left(t\right)}& =\frac{1}{\frac{4\pi }{3}\left({L}^{3}-{R}^{3}\right)}\sum _{n=1}^{\infty }{\mathrm{e}}^{-Dt{\lambda }_{n}}4\pi \underset{R}{\overset{L}{\int }}\mathrm{d}r\enspace {r}^{2}\enspace {u}_{n}\left(r\right)\hfill \\ \hfill & =\sum _{n=1}^{\infty }{c}_{n}\enspace {\mathrm{e}}^{-D{\lambda }_{n}t},\hfill \end{align} \tag{ C21 }$$

with

$$\begin{equation}{c}_{n}=\frac{3{\left(L-R\right)}^{3}}{{L}^{3}-{R}^{3}}\frac{{b}_{n}\left[\left(1+{\alpha }_{n}^{2}\frac{RL}{{\left(L-R\right)}^{2}}\right)\mathrm{sin}\enspace {\alpha }_{n}-{\alpha }_{n}\enspace \mathrm{cos}\enspace {\alpha }_{n}\right]}{{\alpha }_{n}^{5}}.\end{equation} \tag{ C22 }$$

In this appendix, we generalise and rationalise the scaling relation $\overline{\langle {\mathcal{T}}_{N}\rangle }\propto 1/{N}^{2}$ in terms of the closest particle to the target. Let us start with an example of the interval (0, L) with the absorbing endpoint 0 and reflecting endpoint L. For N independent particles randomly placed on the interval at positions x₁, ..., x_N, the distance to the target at 0 is δ = min{x₁, ..., x_N}. This is a random variable characterised by the simple law: $\mathbb{P}\left\{\delta { >}x\right\}={\left(\mathbb{P}\left\{{x}_{1}{ >}x\right\}\right)}^{N}={\left(1-x/L\right)}^{N}$, where we assumed the uniform distribution for x_i. As a consequence, the mean shortest distance is

$$\begin{equation}\overline{\delta }=\underset{0}{\overset{L}{\int }}\enspace \mathrm{d}x\enspace x\enspace \left(-{\partial }_{x}\mathbb{P}\left\{\delta { >}x\right\}\right)=\underset{0}{\overset{L}{\int }}\enspace \mathrm{d}x\enspace \mathbb{P}\left\{\delta { >}x\right\}=\frac{L}{N+1}.\end{equation} \tag{ D1 }$$

Similarly, the second, the third, etc closest particle are located on average at the distance 2L/(N + 1), 3L/(N + 1), etc. This is expected: if one had to place N particles at equal distances between any two neighbours (and from the endpoints), one would get their positions at kL/(N + 1), with k = 1, 2, ..., N. If N is large, there are many particles that start at the distance of the order L/N from the target, and there is high chance that one of these particles reaches the target within the time of the order of ${\overline{\delta }}^{2}/D$, in agreement with equation (22), up to a numerical prefactor.

The above argument can be easily extended to arbitrary bounded domains with a target having a smooth boundary Γ. Here, we define again the shortest distance to the target as δ = min{|x₁ − Γ|, ..., |x_N − Γ|}. This random variable is again characterised by the law $\mathbb{P}\left\{\delta { >}x\right\}={\left(\mathbb{P}\left\{\vert {\mathbf{x}}_{1}-{\Gamma}\vert { >}x\right\}\right)}^{N}$ due to the independence of the starting points. Even though the probability law for the distance from a single particle is now complicated, its asymptotic behaviour remains simple. Indeed, for small x, $\mathbb{P}\left\{\vert {\mathbf{x}}_{1}-{\Gamma}\vert {< }x\right\}$ is the probability of locating x₁ within a thin boundary layer of width x near the target Γ. For smooth Γ and small x, the volume of this layer can be estimated as x|Γ|, where |Γ| is the surface area of Γ. As a consequence, $\mathbb{P}\left\{\vert {\mathbf{x}}_{1}-{\Gamma}\vert {< }x\right\}\simeq x\vert {\Gamma}\vert /\vert {\Omega}\vert$ as x → 0. The mean shortest distance is then

$$\begin{align*}\hfill \overline{\delta }& =\underset{0}{\overset{L}{\int }}\enspace \mathrm{d}x\enspace x\enspace \left(-{\partial }_{x}\mathbb{P}\left\{\delta { >}x\right\}\right)=\underset{0}{\overset{L}{\int }}\enspace \mathrm{d}x\enspace \mathbb{P}\left\{\delta { >}x\right\}\hfill \\ \hfill & =\underset{0}{\overset{L}{\int }}\enspace \mathrm{d}x\enspace {\left(1-\mathbb{P}\left\{\vert {\mathbf{x}}_{1}-{\Gamma}\vert {< }x\right\}\right)}^{N},\hfill \end{align*}$$

where L is the maximal distance from the target: L = max_x∈Ω|x − Γ|. For large N, the main contribution comes from x close 0, for which $\mathbb{P}\left\{\vert {\mathbf{x}}_{1}-{\Gamma}\vert {< }x\right\}$ is small. Substituting the above approximation in this region, one gets

$$\begin{equation}\overline{\delta }\simeq \underset{0}{\overset{{L}^{\prime }}{\int }}\enspace \mathrm{d}x\enspace {\left(1-x\vert {\Gamma}\vert /\vert {\Omega}\vert \right)}^{N}\simeq \frac{\vert {\Omega}\vert }{\vert {\Gamma}\vert \left(N+1\right)},\end{equation} \tag{ D2 }$$

where the integral was truncated at L' = |Ω|/|Γ| to ensure that the integrand function remains positive. Similar computation can be performed to estimate that the mean distance to the target from the next-to-the-closest points is of the same order. If $\overline{\delta }$ is small as compared to the mean curvature of the target, diffusion in the lateral direction does not matter, and the problem is again reduced to diffusion of a single particle on the interval $\left(0,\overline{\delta }\right)$, for which the mean FPT is again ${\overline{\delta }}^{2}/\left(2D\right)$. This rationalises the scaling relation $\overline{\langle {\mathcal{T}}_{N}\rangle }\propto 1/{N}^{2}$ that we derived in equation (22).

An important consequence of this reasoning is that the scaling 1/N² holds far beyond the considered case of the uniform distribution of starting points. First, one can realise that only the behaviour of the distribution near the boundary of the target matters in the limit of large N; indeed, the integrand function ${\left(1-\mathbb{P}\left\{\vert {\mathbf{x}}_{1}-{\Gamma}\vert {< }x\right\}\right)}^{N}$ becomes exponentially small when $\mathbb{P}\left\{\vert {\mathbf{x}}_{1}-{\Gamma}\vert {< }x\right\}$ is not small. Second, even in the vicinity of the target, the distribution does not need to be uniform. For instance, if $\mathbb{P}\left\{\vert {\mathbf{x}}_{1}-{\Gamma}\vert {< }x\right\}\simeq {\left(x/A\right)}^{\alpha }$ with some scale A > 0 and any exponent α > 0, the integral in equation (D2) yields $\overline{\delta }\simeq A/\left(\alpha N+1\right)$, i.e., the scaling 1/N is just affected by a prefactor. This result confirms the generality of the scaling relation (22). This conclusion can also be obtained on a more rigorous basis by considering the heat kernel short-time asymptotic behaviour. In fact, the asymptotic relation (14) for the volume-averaged survival probability holds in a much more general setting, in particular, after the average with a smooth nonuniform distribution of the starting point [114, 116, 117]. We do not explore this direction in the paper.

Only if the distribution of the starting point excludes points near the target, the volume-averaged mean fFPT would scale differently. A standard example is the Dirac distribution that fixes the starting point, which was studied previously and yielded the $\overline{\langle {\mathcal{T}}_{N}\rangle }\propto 1/\mathrm{ln}\enspace N$ behaviour [79] (see also section 1 for more references). A similar behaviour is expected for any distribution that excludes points from some δ₀-vicinity of the target, i.e., if there exists δ₀ > 0 such that $\mathbb{P}\left\{\vert {\mathbf{x}}_{1}-{\Gamma}\vert {< }{\delta }_{0}\right\}=0$. Moreover, this strict exclusion condition can be relaxed. For instance, if $\mathbb{P}\left\{\vert {\mathbf{x}}_{1}-{\Gamma}\vert {< }x\right\}\simeq {\mathrm{e}}^{-A/x}$ as x → 0 (with some scale A > 0), the integral in equation (D2) yields again the 1/ln N behaviour for large N. A more systematic analysis of this problem and finding the sharp exclusion condition, present interesting problems for future research.