Full distribution of first exit times in the narrow escape problem

In the scenario of the narrow escape problem (NEP) a particle diffuses in a finite container and eventually leaves it through a small"escape window"in the otherwise impermeable boundary, once it arrives to this window and over-passes an entropic barrier at the entrance to it. This generic problem is mathematically identical to that of a diffusion-mediated reaction with a partially-reactive site on the container's boundary. Considerable knowledge is available on the dependence of the mean first-reaction time (FRT) on the pertinent parameters. We here go a distinct step further and derive the full FRT distribution for the NEP. We demonstrate that typical FRTs may be orders of magnitude shorter than the mean one, thus resulting in a strong defocusing of characteristic temporal scales. We unveil the geometry-control of the typical times, emphasising the role of the initial distance to the target as a decisive parameter. A crucial finding is the further FRT defocusing due to the barrier, necessitating repeated escape or reaction attempts interspersed with bulk excursions. These results add new perspectives and offer a broad comprehension of various features of the by-now classical NEP that are relevant for numerous biological and technological systems.


Introduction
The Narrow Escape Problem (NEP) describes the search by a diffusing particle for a small 'escape window' (EW) in the otherwise impermeable boundary of a finite domain (figure 1) [1,2]. The NEP represents a prototypical scenario, inter alia, in biophysics, biochemistry, molecular and cell biology, and in nanotechnology. Specifically, the particle could be an ion, a chemically active molecule or a protein confined in a biological cell, a cellular compartment or a microvesicle. The EW could be a veritable 'hole' in the boundary, for instance, a membrane pore, but it could also represent a reactive site right inside (or on) the boundary, such as a protein receptor waiting to be triggered by a specific diffusing molecule. Similarly, the NEP may pertain to a tracer molecule trying to leave from an interstice in the porous structure of an underground aquifer. When one is solely interested in how the particle reaches the target, the NEP corresponds to the first-passage time problem [3]. The situation becomes more complicated when some further action is needed after the EW is reached. For instance, the particle may have to react chemically with an imperfect, partially-reactive site, which will require the crossing of an activation energy barrier. To produce a successful reaction event the particle may repeatedly need to revisit the target after additional rounds of bulk diffusion. Similarly, when the particle has to leave the domain through a pore or a channel, an entropy barrier due to the geometric confinement at the entrance to the EW needs to be crossed [4]. In these scenarios, the relevant quantities are the first-reaction or first-exit times, which are often substantially longer than the first-passage time to the EW. For brevity, in what follows we will call the firstpassage times to the successful reaction event as the first-reaction time (FRT), for both the eventual passage through the EW constrained by an entropic barrier or for the reaction with a partially reactive site. This is a random variable dependent on a variety of physical parameters and here we unveil its non-trivial statistical properties.
The paradigmatic NEP has been extensively studied in terms of mean FRTs [2,[5][6][7]. While early analyses concentrated on idealised, spherically-symmetric domains enclosed by a perfect hard-wall with a perfect (that is, barrierless) EW [8][9][10][11][12][13][14], more recent work examined the NEP for domains with more complicated geometries and non-smooth boundaries, when, e.g. a perfect barrierless target is hidden in a cusp or screened by surface irregularities [15][16][17][18][19]. A striking example of the NEP with very severe escape restrictions is the dire strait problem [20]. Effects of short-ranged (contact) and long-ranged interactions with the confining boundary, which are quite ubiquitous and give rise to 'intermittent motion' characterised by alternating phases of bulk and surface diffusion, were shown to allow for an optimisation of the search times [4,[21][22][23][24][25][26]. It was also studied how molecular crowding in cellular environments impacts the search dynamics [13,27,28]. When the target is imperfect the contributions due to diffusive search and penetration through a barrier were shown to enter additively to the mean FRT [4,[29][30][31], akin to the celebrated Collins-Kimball relation in chemical kinetics [32]. Since these two rate-controlling factors disentangle the concept of search-versus barrier-controlled NEP was proposed [4,30].
Despite of this considerable body of works published on the NEP, information beyond the mean FRT is scarce, only two recent analytical and numerical works consider the full distributions of first-passage times to perfect targets and only in two-dimensional settings [33,34]. In fact, the NEP typically involves a broad spectrum of time scales giving rise to a rich and interesting structure of the probability density function (PDF) of FRTs with different time regimes. Knowledge of this full PDF is therefore a pressing case. As we show here the FRT may be strongly defocused over several orders of magnitude, effecting noticeable trajectory-to-trajectory fluctuations [35]. One therefore cannot expect that, in general, solely the first moment of the PDF can be sufficient to characterise its form exhaustively well. Moreover, it is quite common that the behaviour of the positive moments of a distribution is dominated by its long-time tail and, hence, stems from anomalously long and rarely observed trajectories [35][36][37]. Fluctuations between realisations are indeed common in diffusive processes and can be characterised by the amplitude fluctuations of time-averaged moments [38,39] and single-trajectory power spectra [40][41][42].
We here present the PDF of the FRTs for the NEP in the generic case of a spherical domain with an imperfect target and an arbitrary fixed starting point of the diffusing particle. Indeed, in cellular environments, this setting represents one of the most interesting applications of the NEP, as particles such as small signalling molecules or proteins are usually released at some fixed point inside the cell [43] and need to diffusively locate specific targets such as channel or receptor proteins embedded in the cell membrane [44,45]. Our analytical approach is based on a self-consistent closure scheme [4,46,47] that yields analytical results in excellent agreement with numerical solutions. We demonstrate that in these settings the PDF exhibits typically four distinct temporal regimes delimited by three relevant time scales, for which we also present explicit expressions. These characteristic times may differ from each other by several orders of magnitude. We fully characterise the PDF of FRTs by identifying the cumulative reaction depths (the fractions of FRT events corresponding to each temporal regime of the PDF) and analyse their dependence on the system parameters. Overall, our analysis provides a first complete and comprehensive framework for the understanding of multiple facets of the biologically and technologically relevant NEP, paving a way for a better theoretical understanding of the NEP and a meaningful interpretation of experimental data. Last but not least, it permits for a straightforward derivation of the FRT PDF for a more Figure 1. Spherical domain of radius R containing a target (a partially reactive site or an EW) depicted as a spherical cap (in magenta) of a polar angle e located at the North pole. The red dot denotes the starting point (r, θ) of the diffusing particle, and the zigzag line depicts a random diffusive path to the target. general problem with several diffusive particles starting their random motion from given locations. In appendices, we also discuss the forms of the PDF and the relevant time scales for other possible initial conditions.

The NEP
As depicted in figure 1 we consider a diffusing point-like particle with diffusion coefficient D inside a spherical domain of radius R. The boundary enclosing the domain is perfectly reflecting everywhere, except for the imperfect target-a spherical cap of angle e at the North pole. By symmetry, the behaviour is independent of the azimuthal angle, and we need to solve the diffusion equation for the survival probability S(r, θ; t) that a particle released at position r=(r, θ) at time 0 has not reacted with the partially-reactive site or has not escaped through the EW up to time t. The diffusion equation is completed by the initial condition S(r, θ; t=0)=1 and by the mixed boundary conditions of zero-current at the hard wall combined with the standard Collins-Kimball partially-reactive (or partially-reflecting) boundary condition imposed at the EW or the chemically-active site [32] (see also [48,49] for an overview) The latter Robin boundary condition signifies that the reaction is a two-stage process consisting of (a) the diffusive transport of the particle to the vicinity of the target and (b) the eventual imperfect reaction which happens with a finite probability. The proportionality factor κ here is the intrinsic reactivity, which is defined as κ=ωl, where l is the effective 'thickness' of the reaction zone in the vicinity of the target-the minimal reaction distance, while ω is the rate (frequency) describing the number of elementary reaction acts per unit of time within the volume of the reaction zone. Clearly, ω (and hence, κ) depends on the amplitude of the barrier against the reaction event: w = ¥ (k = ¥) corresponds to the case of a perfect barrierless reaction or an immediate escape upon the first encounter with the EW. In this case the FRT reduces to the first-passage time to the target. Conversely, ω=0 (κ=0) implies that the reaction/escape event is completely suppressed. In our calculations, we suppose that κ is a given parameter(for heterogeneous, space-dependent reactivity, see [50] and references therein). In section 4, we will elaborate on possible distinctions between its values for reaction and escape processes. Once S(r, θ; t) is obtained, the PDF H(r, θ; t) of the FRT is given by the negative derivative of the survival probability with respect to t. Further details on the NEP, our analytical derivation, and the comparison with numerical results are presented in appendices.

Results
Our principal analytical result is the explicit, approximate expression for the PDF H(r, θ; t) of the FRT, where f e e = P 1, cos n 0 ( ) ( ) are the Legendre polynomials, j n (z) are the spherical Bessel functions of the first kind, the prime denotes the derivative with respect to the argument, while α m are the positive solutions (organised in ascending order) of the transcendental equation Equations (2)-(5) completely determine the explicit form of H(r, θ; t). As we already mentioned above, this result is obtained by resorting to a self-consistent closure scheme, developed earlier for the calculation of the mean FRT in certain reaction-diffusion problems [4,46,47]. This approximation consists in replacing the actual mixed boundary condition (1) by an inhomogeneous Neumann condition and in the derivation of an appropriate closure relation, which ensures that the mixed boundary condition (1) holds on average. The adaptation of such a scheme to the calculation of the full PDF in NEP is discussed in detail in appendices A and B. We note that despite an approximate character of our approach, the obtained result turns out to be very accurate for arbitrary starting position (r, θ) (although not too close to the target), for a very wide range of ε, and for arbitrary value of κ. The accuracy was checked by comparing our approximation to a numerical solution of the original problem by a finite-elements method (FEM, see appendix E). We note, as well, that two important characteristic length-scales can be read off directly from expression (2) and are associated with the smallest and the second smallest solutions of (5), (that being, α 0 and α 1 ), which define the largest and the second largest decay times of the PDF, respectively. The analysis of all the solutions of (5), which define the full form of the PDF valid at arbitrary times, their relation to the characteristic time scales, the derivations of the asymptotic forms of H(r, θ; t), and its easy-to-implement approximate mathematical form are provided in the appendices. The latter also present a variety of other useful results based on H(r, θ; t), e.g. the surface-averaged and volume-averaged PDFs (see equations (B.14) and (B.15)).
In figure 2 we depict the PDF H(r, θ; t) defined by (2)-(5) for three different fixed initial positions (polar angles θ=0 (a), θ=π/2 (b), and θ=π (c), with fixed r/R=0.9), four values of the intrinsic reactivity κ, and for a particular choice e = 0.1. We observe that the PDF has a clearly defined structure, which depends on whether the starting point is close to the target or not. In the former case (as in panel (a)), the PDF consists of a hump-like region around the most probable FRT t mp , an extended plateau-like region after the crossover time t c (remarkably, within this region all values of the FRT are nearly equally probable), and, ultimately, exponential decay starting right after the mean FRT t mean . Note that up to t mean the particle experiences multiple collisions with the boundary of the container after unsuccessful reactions, thereby loosing the memory about its precise starting location. Due to this circumstance, the ultimate, long-time exponential decay (which is the most trivial part of the PDF), holds for all bounded domains, regardless of the precise initial condition (see also appendices). As evidenced in figure 2 the difference between t mp , t c and t mean may span orders of magnitude revealing a pronounced defocusing of the FRTs (see below). In contrast, when the starting point is far from the target (panels (b) and (c)), the hump-like region is much less apparent, and the extended plateau and the exponential decay are the major features of the PDF. Although the maximum of the PDF formally exists, the corresponding most probable time is not informative because the escape (or reaction) may happen at any moment over the plateau region with almost equal probabilities. We note that the overall shape of the curves presented in figure 2 is generic and gets only slightly modified upon parameter changes-see, e.g. figure A1 in appendices in which a similar plot for = r R 0.5 is presented. Note, as well, that the behaviour of the PDF in the case when the starting point is uniformly distributed over a spherical surface some distance r away from the origin, appears to be very similar to the one depicted in figure 2 in panels (b) and (c) (see figure C1, panel (b)). Conversely, in the case when the location of the starting point is random and uniformly distributed over the volume of the container, the PDF consists essentially of a plateau-like region and the ultimate exponential tail (see figure C1(a)).
The existence of the hump-like region is not an unexpected feature in the case of a fixed initial condition. For a particle initially placed some distance away from the target, it is impossible to reach the latter instantaneously such that the probability of having a very short FRT is close to zero. This defines the most probable reaction time t mp , which can be estimated as function of the initial distance σ to the target (in panel (a) of figure 2, s = -R r as θ=0) in the form (see appendices), , at which the particle first engages with the boundaries. The coloured bar-codes above each main panel indicate the cumulative reaction depths corresponding to the four values of k , in decreasing order from top to bottom. Each bar-code is split into ten regions of alternating brightness, representing the ten 10%-quantiles of the PDF. The numbers on top of the bar-codes indicate the values of Dt R 2 . For example, in panel (a) the first dark blue region for the case k = ¥ (perfect reaction or no entropic barrier at the EW) indicates that 10% of reaction events occur up to »´-Dt R 3 10 2 3 , which is close to the dimensionless most probable FRT Dt R mp 2 . The mean FRT t mean in this case is almost four orders of magnitude longer than the most probable FRT t mp , and over 70% of trajectories arrive to the target up to this time. Analogous estimates of the cumulative reaction depths for other initial conditions are presented in the text below and also in appendices. Compare figure A1 for the case of r/R = 0.5.
The most probable FTR corresponds to 'direct' trajectories [51] on which the particle moves fairly straight to the target, with immediate successful reaction. Consequently the value of t mp is 'geometry-controlled' [51,52] by the initial distance R−r. At such short time scales, the particle has not yet explored the boundaries of the domain, and thus the initial increase up to t mp and the subsequent power-law descent of the PDF are well-described by the Lévy-Smirnov-type law a -A t t exp 3 2 ( ) [3] (see appendix C.2 for more details). In our case of an imperfect target, the amplitude A also acquires the dependence on the reactivity κ and the target radius ρ. In particular, the descent from the peak value of the PDF obeys The independence on the domain size stems from the fact that the particle at such short times has no knowledge about the boundaries.
The hump-like region is delimited by the crossover time an important characteristic time-scale, at which the particle first engages with the domain boundaries and realises that it lives in a bounded domain. The plateau-like regime, which sets in at t c , is a remarkable feature for the NEP that has not been reported in literature on this important and generic problem. The existence of such a plateau for diffusion-reaction systems with an activation energy barrier at the target has only been recently discovered in other geometrical settings, such as, e.g. the diffusive search for a partially-reactive site placed on a hard-wall cylinder confined in a cylindrical container [47] and a spherical target at the centre of a larger spherical domain [52]. The fact that it also appears in typical NEP settings, as shown here, evidences that it is a Figure A1. PDF H(r, θ; t) of the FRT, determined from (2) to (5) and rescaled by R D 2 , as function Dt R 2 for e = = r R 0.5, 0.1, several values of dimensionless reactivity k k = R D (indicated in the plots), and θ=0 (a), θ=π/2 (b), θ=π (c). Self-consistent approximation is truncated to n=50 terms Vertical arrows indicate the mean FRT, t mean . Coloured bar-codes above each figure indicate the cumulative depths corresponding to four considered values of k , in decreasing order from top to bottom. Each bar-code is split into ten regions of alternating brightness, each representing ten 10%-quantiles of the distribution. The numbers on top of the bar-codes indicate the values of Dt R 2 . For example, in panel (a) the first dark blue region for the case k = ¥ (perfect reaction) indicates that 10% of reaction events occur till´-Dt R 3 10 2 1  . The dimensionless mean FRT in this case is around 10 1 , i.e. is almost two orders of magnitude bigger that Dt R mp 2 , and over 70% of all trajectories arrive to the target up to this time.
and several values of κ. Self-consistent approximation is truncated to n=100 terms Arrows indicate the mean FRT. Thin lines show the short-time asymptotic behaviours (C.11) and (C.14) for two plots, respectively. fundamental feature of systems in which the reaction event requires not only a diffusive search but also a barrier crossing.
The plateau typically persists for several decades in time, up to the mean FRT t mean (indicated by vertical arrows in figure 2), which is very close (but not exactly equal) to the decay-time ), associated with the smallest eigenvalue λ 0 of the Laplace operator (see appendices for more details). As a consequence, the FRT PDF features the exponential shoulder q~- Note that the right-tail of the PDF is thus fully characterised by the unique time scale t d . For a perfectly-reactive target, the asymptotic behaviour of t d in the narrow-escape limit (i.e. for e  0) was first obtained by Lord Rayleigh [8] more than a century ago within the context of the theory of sound, and then refined and adapted for the NEP by Singer et al [10], yielding (see also discussions in appendices and in [53][54][55]) For an imperfect target, however, this asymptotic result is no longer valid and one finds instead that (see appendices) which signifies that the decay time diverges much faster when e  0, as compared to the case of a perfect barrierless target. Lastly, we note that the fact that t mean and t d appear very close to each other is an apparent consequence of the existence of the plateau region with an exponential cut-off, suggesting that the major contribution to t mean (and most likely to all positive moments of the PDF) is dominated by the integration over this temporal regime.

Discussion
The results presented above reveal many interesting and novel features of the NEP that we now put into more general context. (i) In case of chemical reactions with a partially-reactive site,the intrinsic reactivity κ is independent of e. As a consequence, t d in equation (11) (as well as t mean , see appendices) diverges in the NEP limit as ẽ t 1 d 2 . The situation is much more delicate in the case of an entropy barrier ΔS. The impeding effect of the latter has been studied in various guises, including transport in narrow channels with corrugated walls (see, e.g. [56][57][58][59][60]) and also for the NEP (see, e.g. [61,62]), although for the latter problem no explicit results have been presented. In particular, for diffusion in channels represented as a periodic sequence of broad chambers and bottlenecks, the influential works [56,58] expressed this barrier through a profile h(x) of the channel cross-section,while in [60], for ratchet-like channels, this barrier was estimated as where β is the reciprocal temperature, h max and h min are the maximal and minimal channel apertures (a chamber versus a bottleneck). Therefore, one may expect that for the NEP the entropic barrier scales as b a e D~-S ln 1 1 ( ), where α>0 is a numerical constant characterising the precise form of the profile h(x). Assuming that the particle penetrates through the barrier solely due to a thermal activation,i.e. that k b -DS exp( ), one expects then that k ẽ a and thus vanishes in the NEP limit e  0. This means,in turn,that here ẽ i.e. t d diverges more strongly in the NEP limit than its counterpart in the case of a partially-reactive site. As a consequence, the overall effect of such an entropy barrier on the size of the plateau region and on the cumulative reaction depth on this temporal stage, is much more significant for the NEP than for diffusion-reaction systems. Importantly, both for the NEP and for reactions with a partially-reactive site, as e  0, the mean FRT is always controlled by the penetration through the barrier rather than by the stage of random search for the target whose contribution to t d diverges in the leading order only as e 1 . We also remark that, of course, the NEP limit has to be understood with an appropriate caution: taking e  0, one descends to molecular scales at which both the molecular structure of the container's boundaries, surface irregularities, presence of contaminants and other chemically active species come into play. At such scales, the solvent can also be structured close to the wall, affecting characteristic diffusion times and hence, the value of the diffusion coefficient D. Moreover, the size of a particle does matter here and one can no longer consider it as point-like. The entropic barrier is then a function of the ratio between the particle radius a and the radius ρ of the EW. Understanding theoretically the functional form of κ is so far beyond reach and one may address the problem only by numerical simulations or experiments. We are, however, unaware of any systematic approach to these issues, except for a recent paper [63], which studied the particle-size effect on the entropic barriers via numerical simulations. It was observed that the dependence of the inverse transition rate (and hence, of κ) of a Brownian particle in a model system consisting of a necklace of hard-wall spheres connected to each other by EWs, can be well-fitted by an empirical law of the form k r a 1 3 2 ( ) signifying that κ vanishes very rapidly when r  a . When the particle size a is comparable to that of the EW, the entropic barrier becomes very high resulting in a significant decrease of κ. Some other effects of the particle size and shape onto diffusive search, unrelated to our work, were discussed in [64,65].
(ii) The concept of the search-controlled versus barrier-controlled NEP proposed in [4,30] holds indeed for the mean FRT, see equation (C.6). In contrast, we observe that two other important characteristic time-scales-t mp and t c -are unaffected by the reactivity and solely depend on the diffusion coefficient and the geometry. This means that the two controlling factors represented by κ and D do not disentangle. Moreover, both affect the shape of the PDF and the corresponding cumulative reaction depths.
(iii) For the NEP with fixed starting point of the particle the characteristic time-scales t mp , t c and t mean differ by orders of magnitude meaning that the FRT are typically very defocused. This is a conceptually important point for the understanding of the NEP, but alone it does not provide a complete picture. Indeed, one should also analyse the cumulative reaction depths corresponding to each temporal regime. For example, for the settings used in figure 2 (panel (a)), the shortest relevant time-scale is t mp . For perfect, barrierless reactions the amount of reaction events appearing up to this quite short time ) is surprisingly non negligible, being of the order of 10%. Upon lowering the dimensionless reactivity k k = R D (i.e. upon increase of the activation barrier at the target or of the diffusion coefficient D, or upon lowering R), this amount drops significantly and the contribution of 'direct' trajectories becomes less important. This is quite expected because the latter now very rarely lead to the reaction event. In turn, the contribution of the whole hump-like region (for t stretching up to t c ) for barrierless reactions is more than 50%, meaning that more than a half of all reaction events take place up to the time-scale when the particle first realises that it moves in a confined domain. One infers next from figure 2 (panel (a)) that upon lowering k , this amount also drops rapidly. At the same time, the plateau-like region becomes increasingly more pronounced for smaller k and progressively more reaction events take place during this stage: we find 25%, 30%, 47% and 55% of all reaction events for k = ¥, 10, 1 and 0.1, respectively. Interestingly, the fraction of 'unsuccessful trajectories', which survive up to the largest characteristic time t mean , appears to be very weakly dependent on the value of k , and is rather universally about 30%. Note that this number is close to the fraction of outcomes of an exponentially distributed random variable above its mean, exp 1 0.37 ( )  . Therefore, even for this particular case when the starting point is very close to the target, one observes a rather diverse behaviour-there is no apparent unique time-scale and no overall dominant temporal regime: the region with the rise to the most probable value, the one with the subsequent descent to the plateau, the plateau-like region itself, and the eventual exponential decay may all become relevant for different values of the parameters and hence, their respective contributions depend on the case at hand. This illustrates our earlier statement that the knowledge of the full PDF of FRT is indispensable for the proper understanding of distinct facets of the NEP.
For a random starting point, when it is either uniformly distributed over the volume, or is located at a random point on a distance r away from the origin, the situation is somewhat different. Here, the hump-like region is smeared away and the PDF consists of just two regions: a plateau-like region and the ultimate exponential decay. As a consequence, in these cases t mean is a unique characteristic time-scale, as it was claimed previously (see, e.g. [7]). In appendices, we present explicit results for the general case of diffusion-mediated reactions with an activation energy barrier and imperfect EW with an entropy barrier. Concerning the cumulative reaction depths in case of a random starting point, we observe an analogous surprising behaviour to what was noted above: the cumulative reaction depth over the plateau region (up to t mean ) appears to be rather universally (i.e. independent of κ and other parameters) equal to ≈70%, while the remaining ≈30% of trajectories react within the final exponential stage.
(iv) Apart from the analysis of the cumulative depths, it is helpful to work out an explicit approximate formula for the PDF, which may be used to fit experimental or numerical data. In appendices, we outline a derivation of such a formula, which agrees very well for the whole range of variation of the FRT.
(v) One often deals with situations when several diffusing particles, starting at some fixed points, are present in a biophysical system and the desired random event is triggered by the one, which arrives first to a specific site on the surface of the container. Among many examples (see e.g. [66]), one may mention, for instance, calcium ions which activate calcium release in the endoplasmic reticulum at neuronal synapses. Our results for the fixed starting-point PDF H(r, θ; t) and the corresponding survival probability provide the desired solution to this problem too. As a matter of fact, if the particles are present at a nanomolar concentration, such that their dynamics are independent, the probability that neither of the particles arrived to the target site up to time t is just the product of the single-particle survival probabilities. The PDF of the time of the first reaction event then follows by a mere differentiation of this product.

Conclusion
With modern optical microscopy, production events of single proteins [43], their passage through the cell [67] as well as single-molecule binding events [68,69] can be monitored in living biological cells [70]. Typical processes in cells such as gene expression involve the production of proteins at a specific location, their diffusive search for their designated binding site, and the ultimate reaction with this partially-reactive site [71,72]. Given the often minute, nanomolar concentrations of specific signalling molecules in cells, the FRT PDF of a diffusing molecule with its target in the cell volume was shown to be strongly defocused and geometry-controlled [47,51,52]. These analyses demonstrate that the concept of mean FRTs (or reaction rates) is no longer adequate in such settings involving very low concentrations. Here one invariably needs knowledge of the full PDF to describe such systems faithfully.
Complementary to this situation when a particle is released inside a finite volume and needs to react with a target somewhere else in this volume, we here analysed the important case when a diffusing particle is released inside a bounded volume and needs to react with a target on the boundary. Our results for the FRT PDF including the passage of an energetic or entropic barrier leading to final reaction with the target or the crossing of the EW, demonstrate a pronounced defocusing of time scales. While the most likely FRT corresponds to a direct particle trajectory to the target and immediate reaction, indirect paths decorrelating the initial position by interaction with the boundary are further emphasised by unsuccessful reaction events, leading to further rounds of exploration of the bulk. This effect is shown to imply a pronounced defocusing, and thus a large reaction depth is reached long before the mean FRT. In particular, we demonstrate that the shapes of the PDF are distinctly different from the case when the particle needs to react with a target inside the volume of the bounded domain.
Our findings have immediate relevance to biological cells, in which specific molecules released in the cell need to locate and bind to receptors on the inside of the cell wall, or they need to pass the cell wall or compartment membrane through protein channels or membrane pores. Similar effects occur in inanimate systems, when chemicals in porous aquifers need to leave interstices through narrow holes in order to penetrate further. From a technological perspective our results will find applications in setups involving lipid bilayer micro-reaction containers linked by tubes [73].

Appendix A. Solution in Laplace domain
Our starting point is the diffusion equation for the survival probability S(r, θ; t) of a particle started at position (r, θ) inside a sphere of radius R, up to time t: subject to the initial condition S(r, θ; t=0)=1 and the mixed boundary conditions: where ∂ n =∂ r is the normal derivative directed outwards the ball, Δ the Laplace operator, D the diffusion coefficient, κ the intrinsic reactivity, and e the angular size of the EW (or target site) at the North pole ( figure 1).
where Q is an effective flux to be computed, and Θ is the Heaviside step function.
We where I ν (z) is the modified Bessel function of the first kind, and we fixed a particular normalisation (as shown below, the result does not depend on the normalisation). Substituting (A.7) into the modified boundary condition in (A.6), multiplying by q q P cos sin n ( ) and integrating over θ from 0 to π, one gets  As mentioned above, the solution does not depend on the normalisation of functions g n (r), as they always enter as a ratio of ¢ g n and g n . One can check that the PDF H(r, θ; t) is correctly normalised, i.e. q = = H r p , ; 0 1 ( ) . In fact, one has in the limit p→0:  In appendix E, we compare these approximate solutions to the numerical solution of the original problem obtained by a FEM. We show that the approximate solutions are remarkably accurate for a broad range of p, for both narrow (e = 0.1) and extended (e = 1) targets.
It is worth mentioning that the developed self-consistent approximation can also be applied to the exterior problem when a particle is released outside a ball of radius R and searches for a target on its boundary. The derivation follows the same lines, with two minor changes:

Appendix B. Solution in time domain
To get the PDF H(r, θ; t) in time domain, we need to perform the inverse Laplace transform of q H r p , ; ( ) from (A.21). As diffusion occurs in a bounded domain, the spectrum of the underlying Laplace operator is discrete and is determined by the poles of q H r p , ; ( ), considered as a function of a complex-valued parameter p. We note first that the zeros of the functions g 0 (R) and ¢ g R n ( ), standing in the denominator of q H r p , ; ( ), are not the poles, as the corresponding divergence is compensated by vanishing of the function η p at these points. So, we are left with the poles of the function η p from (A.17) which are determined by the following equation on Î  p : is the smallest strictly positive zero of ¢ j z 1 ( ) (i.e. the right-hand side is the first strictly positive eigenvalue of the Laplace operator for the full reflecting sphere). As expected, this bound is independent of the target size and determines the crossover time t c that we set as (as a qualitative border between two regimes, t c is determined up to an arbitrary numerical factor that we select here to be 1  As expected, when the target is small or its reactivity is small, the Laplacian spectrum is close to that of the fully reflecting sphere.

B.2. Perfectly reactive target
When k = ¥, the right-hand side of (B.4) is zero, and the asymptotic analysis for small e is different. Let us first consider the smallest pole which is expected to be small as e  0. As discussed in [4], the numerical coefficient of this asymptotic behaviour differs by around 10% from the exact coefficient π/3 [10]. With the help of these expressions, we plot figure 2 in the main text and figure A1 here for two different cases of the initial particle radius.

Appendix C. Asymptotic behaviour
In this section, we discuss the long-time and the short-time asymptotic behaviours of the PDF H(r, θ; t).

C.1. Long-time behaviour
As expected for diffusion in a bounded domain, the PDF exhibits an exponential decay at long times, with the rate determined by the pole with the smallest amplitude. The latter is given by the smallest positive solution of (B.4) whose asymptotic behaviour was discussed in appendices B.1 and B.2. In the narrow escape setting, this decay rate is close to the inverse of the volume-averaged mean FRT t mean . In turn, this mean, as well as higherorder moments, can be obtained from q H r p , ; ( ) in the limit  p 0. Using the basic properties of the functions g n (see appendix D), we get As a consequence, the Laplace-transformed PDF behaves as This result generalises our earlier formula for the volume-averaged mean FRT [4], which can be retrieved by averaging over the starting point uniformly distributed in the bulk: As mentioned earlier, the leading term in the latter expression is exact, while the numerical prefactor 32/(9π) in the subdominant term in this approximate relation exceeds by 10% the exact prefactor π/3 (for further details, see [4]).

C.2. Short-time behaviour
The short-time asymptotic behaviour is determined in the limit  ¥ p . We note that  Figure C1(right) illustrates the accuracy of this asymptotic behaviour.
The Heaviside function distinguishes two cases. When the starting point lies within the solid angle of the target ( i.e. q e < ), the distance to the target is just R−r, and this asymptotic formula holds. In turn, when the starting point lies outside this solid angle (i.e. q e > ), the distance to the target is larger than q -R r H r p , , ; ( )should decay faster, the contributions in (C. 18) and (C. 19) are thus cancelled, and one needs more refined asymptotic analysis. Without dwelling further on this analysis, we present below probabilistic arguments to get a reasonable approximation.
When the starting point is far from the small EW, small exit times are extremely unlikely, with the probability density vanishing exponentially fast as is the distance to the EW. To find the correct prefactor to this exponential tail, we note that (i) a small spherical cap, seen from far away, can approximated by a small sphere of the same radius r e = R 2 sin 2 ( ), centred at the North pole; the distance to the boundary of this sphere is s q r = + --R r rR 2 cos ; 2 2 and (ii) the reflecting boundary of the spherical confining domain can be removed because only almost straight trajectories to the target, not touching the confining boundary, do matter in the short-time limit. In other words, the short-  . Figure C2 reproduces the rescaled PDF from figure 2 to illustrate the quality of the short-time asymptotic formulae. For the case θ=0, the panel (a) confirms that the short-time asymptotic formula (C.19) accurately captures the left exponential tail of the PDF up to the most probable time, i.e. for  t t mp . However, this approximation fails at longer times due to the confinement. Note also that the estimation of the most probable time t mp from (C.19), suggested in [52] for a simpler geometric setting, is not accurate here, except for the perfectly reactive target. From figure C2(a), we only can say that < t R r D 6 mp 2 ( ) ( )for the perfectly reactive target (k = ¥), and < t R r D 2 mp 2 ( ) ( )for partially reactive targets (k < ¥). The panels (b) and (c) illustrate the accuracy of another short-time asymptotic formula, given by (C.20), which is applicable when the starting point is far away from the target (θ=π/2 and θ=π, respectively). Once again, the left exponential tail of H(r, θ; t) is captured reasonably well. In this case, there is no hump, and a plateau region emerges immediately after this tail. As the most probable time now belongs to this plateau region, it is very difficult to quantify. This is a radically new feature of the NEP, as compared to recent works [47,52].

C.3. Explicit closed-form approximation
The knowledge of both short-time and long-time asymptotic behaviours of the probability density motivated the search of closed-form approximations for the probability density H(r, θ; t) over the whole range of times. For instance, Isaacson and Newby developed a uniform asymptotic approximation for perfectly reactive small targets in bounded domains using a short-time correction calculated by a pseudopotential approximation [54]. In turn, a linear superposition of the short-time approximation in (C.20) with the long-time exponential decay was proposed in [52] for a partially reactive spherical target located in the centre of a larger confining sphere.
Here, we rationalise and improve the linear superposition by employing the simple idea of splitting trajectories toward a small target into two groups: (i) 'direct' trajectories that do not hit the boundary of the confining domain and thus do not 'know' about its presence, and (ii) 'indirect' trajectories that reach the boundary and thus explore the bounded confining domain until eventual reaction with the target. The statistics of direct trajectories is close to that of trajectories in the unbounded exterior of the target; in particular, if the starting point is far from a small target, the probability density of the first exit (reaction) times is close to s ) given by (C.20). In turn, the indirect trajectories are responsible for the long-time exponential decay of the probability density. If τ 1 and τ 2 denote the conditional first exit times for such direct and indirect trajectories respectively, then one can set τ=τ 1 with the probability q 1 of realising a direct trajectory, and τ=τ 2 with the probability =q q 1 2 1 of realising an indirect trajectory. If the related probability densities ρ 1 (t) and ρ 2 (t) were known, then the probability density of τ would be simply . C.22 The computation of the probability densities ρ 1 (t) and ρ 2 (t) is an even more complicated problem than the original mixed boundary problem for H(r, θ; t). So, we propose another, much simpler splitting into two groups by introducing a heuristic threshold T for the first exit time. In this approximation, all trajectories toward the target that are shorter than T are treated as 'direct', whereas all longer trajectories are treated as 'indirect'. According to this construction, the first exit times of direct trajectories, that can be described by the probability density s ), are conditioned to be shorter than T; in turn, the first exit times of indirect trajectories, that can be described by the probability density m m e t with μ=Dλ 0 being related to the smallest eigenvalue λ 0 of the Laplace operator in the confining domain, are conditioned to be longer than T. The related conditional probability densities are then given as The threshold T can either be fixed to a prescribed value (e.g. T=t c ), or determined by imposing an additional condition, i.e. the continuity of the derivative of ρ(t). Figure C3 illustrates the accuracy of this explicit approximation. Unfortunately, these relations are numerically unstable for ξ below n that prohibits their use even for moderate n.
In order to resolve this numerical problem, we propose the following trick. Rewriting (D.9) as To overcome numerical instabilities, we propose the following algorithm: (i) for any ξ, we select a large enough truncation parameter n max such that x n max | |  (in practice, we set n max =max{ 100, [10ξ]}); (ii) we approximate f n with n=n max according to its Taylor expansion  2  1  2  3  2  3 2  5   2  2  3 2  5 2  7  10  27  2  3 2  5 2  7 which is valid whenever x n | |  (even if x | | is large); (iii) applying (D.10) in a backward direction, we compute all f n with   n n 1 max . When the argument x | | is very large (as compared to n), one can use the asymptotic relation We checked that this algorithm provides very accurate results for real ξ. Qualitatively, the improved numerical stability of recurrent relations in the backward direction resembles the advantage of an implicit Euler scheme for solving differential equations over an explicit one. Note that the computation can be more subtle for some imaginary ξ, at which x coth diverges (i.e. at x p = k i with an integer k). However, we have not experienced such numerical difficulties for considered parameters.

Appendix E. Numerical verification
In order to assess the accuracy of the self-consistent approximation, we solve numerically the original boundary value problem for the modified Helmholtz equation  Once the solution is found, we also evaluate the volume average H p ( ) (by numerical integration over the triangular mesh), and the surface average H p r ( ) (by interpolating the solution to a line r = const and then integrating numerically over θ). We selected the maximal mesh size to be h max =0.01. In order to verify the accuracy of the FEM solution, we solved the problem twice, with two mesh sizes h max and 2h max , and checked that the resulting solutions were barely distinguishable. Figures E1 and E2 compare the volume-averaged and surface-averaged Laplace-transformed PDFs H p ( ) and H p r ( ) , obtained via our self-consistent approximation (lines) and by FEM (symbols). One can see that the self-consistent approximation is remarkably accurate over a broad range of p, for both narrow e = 0.1 and extended e = 1 target region. The accuracy is higher for less reactive targets.  Figure E2. Surface-averaged Laplace-transformed PDF H p r ( ) of the FRT through a spherical cap, with e = = = = R D r 1, 1, 0.5, 0.1 (left) and e = 1 (right), and several values of κ. Self-consistent approximation (lines), truncated to n=1000 terms, is compared to the FEM solution of the original problem (symbols), with h max = 0.01. Note that our theoretical predictions are almost indistinguishable from the FEM solution for several decades of variation of the dimensionless parameter R p D Figure E3 compares the Laplace-transformed PDF q H r p , ; ( ), obtained via our self-consistent approximation (lines) and by FEM (symbols), for the three starting points used in figures 2, C2 and C3. One can see that the self-consistent approximation is remarkably accurate over a broad range of p, for both narrow e = 0.1 and extended e = 1 target region. The accuracy is higher for less reactive targets. Minor deviations seen in panel (c) for the case k = ¥ are caused by the truncation to n=50 terms and can be corrected by increasing n (not shown). In turn, minor deviations in panel (e) cannot be corrected by increasing n and can thus be attributed to a limited accuracy of the self-consistent approximation for extended targets.