Trapping and escape in a turbid medium

Consider a beaker filled with a turbid medium—a fluid that contains many small reactive traps (Fig. 1). We assume that the trap concentration is sufficiently low that interactions between traps can be ignored. Suppose that a molecule diffuses in the fluid and is absorbed whenever the molecule touches the surface of any of the traps or the wall of the beaker. Our basic goal is to compute the probability that the molecule is absorbed by the traps or by the beaker wall. In the latter case, the diffusing molecule can be viewed as escaping from the beaker. This type of system arises naturally in the absorption of photons in a multiple scattering medium, where the extensive literature has discussed absorption and escape (see, e.g., Refs. 1–7).

For simplicity, we assume that the traps are identical spherical particles whose radii a > 0 are small compared to the beaker radius. To simplify matters, we view the traps as immobile. This is a reasonable assumption if the traps are much larger than the diffusing molecules, and the treatment of the general case is essentially identical as we will discuss. We assume that the number of traps is sufficiently large that fluctuations in their local density can be ignored. In this limit, it is generally possible to replace the discrete traps by an effective average trapping medium. The dynamics of the molecules can then be described by a reaction-diffusion equation, for which explicit solutions can be obtained by standard methods. The converse situation where the beaker contains a few traps or the trapping efficiency of each trap is distributed is challenging because the trapping probability depends in detail on the trap positions, and the effective-medium approach no longer applies.

Eventually all the molecules are absorbed by either the traps or by the beaker wall. What is the fraction E that reaches the beaker wall? What is the fraction T that gets absorbed by the traps? Our goal is to calculate these escape and trapping probabilities, E and T, respectively, for simple geometries. We will first study the case of a spherical beaker of radius R in three dimensions and then generalize to arbitrary spatial dimensions d and to beakers of other simple shapes, such as a parallelepiped and a cylinder.

We begin, in Sec. II, by writing the dynamical equation and governing the escape and trapping probabilities and develop an equivalent time-independent description for these probabilities. Although the molecule density evolves with time, it is possible to recast the infinite-time escape and trapping probabilities as a simpler time-independent problem. In Sec. III, we use this perspective to solve the idealized situation of a single trap of radius a at the center of a beaker of radius R. We also extend this solution to a non-concentric geometry in two dimensions.

In Sec. IV, we treat our primary example of the escape of a diffusing molecule from a spherical beaker that contains N $≫$ 1 traps. In principle, we have to solve the diffusion equation in the available volume—the region exterior to the traps and interior to the beaker—and then compute the time integrated fluxes to the traps and to the beaker walls. A full solution to this formidable problem is not feasible, and instead, we apply the powerful reaction-rate approach (RRA).^8–14 Here we replace the discrete traps by a spatially uniform effective trapping medium, leading to a problem that is readily soluble by standard methods.

Generically, one might anticipate that the exit probability should depend on two dimensionless parameters: N and a/R. We shall show, however, that the exit probability E(N, a/R) depends only on a single parameter, which is the product of N and a/R,

If one knew in advance that the behavior is determined by a single parameter, one might anticipate that it might be N(a/R)³, the ratio of the trap volume to the system volume, or perhaps N(a/R)², the ratio of the total surface of the traps to the total exit area. Thus the dependence on a single parameter, specifically on λ = Na/R, could be puzzling at first sight; we shall show, however, that this dependence follows naturally from the structure of the RRA.

In Sec. V, we obtain the explicit solution for the exit probability E(λ) for a spherical beaker in three dimensions by the RRA. The exit probability is given by the following compact formula:

The form of the exit probability is sensitive to the shape of the beaker, but the limiting behaviors

are universal, and they generally apply for non-pathological beaker shapes; only the amplitudes of these scaling laws are shape dependent. It bears emphasizing that the volume fraction occupied by the traps φ = N(a/R)³ is generally small in realistic systems. However, since λ = N^2/3φ^1/3, the parameter λ does not have to be small, even for φ$≪$ 1. Thus a dilute system of traps can be strongly absorbing. A similar situation arises in the case of small absorbing receptors on the surface of an otherwise reflecting sphere.¹⁵ This system can have an absorption coefficient that is nearly the same as that of a perfectly absorbing sphere for a small receptor density.

We seek the exit probability E(r₀) for a molecule to be absorbed at the surface of the beaker, given that the molecule starts at r₀. There is also the complementary trapping probability T(r₀) for a molecule to be absorbed at the surface of one of the traps. To determine the exit probability, we should solve the diffusion equation

for the molecular density ρ(r, t), with the initial condition ρ(r, t = 0) = δ(r − r₀), and with absorbing boundary conditions ρ = 0 on the surfaces of all traps and on the beaker wall. The exit probability E(r₀) equals the local flux −D$∇$ρ, integrated over the beaker surface and over all time. If we assume that the molecules are uniformly distributed throughout the beaker, the exit probability, averaged over the initial molecular positions, is

where the integrals extend over the allowed spatial region of the molecules.

Because the time does not appear in E(r₀), the problem can be recast into a time-independent form by integrating Eq. (4) over all time. By this approach, one finds that E(r₀) satisfies the backward Kolmogorov equation^10,11

where the subscript 0 indicates that the derivatives in the Laplacian refer to the initial coordinates of the molecule. This backward equation should be solved subject to the boundary conditions

The first condition merely states that a molecule starting on a trap surface cannot reach the beaker wall, while the second condition states that a molecule starting on the beaker wall necessarily exits. The exit probability thus satisfies the Laplace equation for the electrostatic potential in the same geometry, with the traps and the beaker wall substituted by conductors. We may now utilize well-known results for the corresponding electrostatic system^16,17 to determine the exit probability.

We first determine the exit probability E for a single spherical trap of radius a centered inside a spherical beaker of radius R. For ease of notation, we define E₋(r) as the probability that a molecule that starts at distance r from the origin will be absorbed by the inner sphere (the trap); thus E₋(r) corresponds to the trapping probability defined above. Similarly, let E₊(r) be the probability to first reach the outer sphere—the escape probability defined previously. As discussed above, E₊ and E₋ satisfy the Laplace equation

subject to the boundary conditions

By spherical symmetry, we only need to solve the radial part of the Laplace equation (7) subject to these boundary conditions, from which the exit probability is¹⁰ 

with λ = a/R. Since there is only a single “trap,” the quantity λ defined here is equivalent to that given in Eq. (1). The trapping probability is the complementary quantity: E₋(r) = 1 − E₊(r).

To compute the exit probability averaged over all molecules, we assume that they are uniformly distributed in the annulus a < r < R, so that the density in the range (r, r + dr) is proportional to r^d−1dr. The total exit probability is then

Using (8) and (9) gives, in three dimensions,

Notice that the fraction of molecules that exit from the beaker decreases from 1 to 1/2 as λ increases from 0 to 1 (Fig. 2). When λ → 0, the trap at the origin has negligible size so that nearly all molecules will escape the beaker. Conversely, as λ → 1, the physical region becomes an infinitesimally thin annulus, which can be approximated as a slab. In this case, the molecule is equally likely to escape via either boundary. In two dimensions, the solution to the Laplace equation (7) is again straightforward to obtain and the corresponding result for the exit probability is (Fig. 2)

This expression has the same limiting behaviors as in the case of d$≠$ 2.

We can also extend the above approach to a trap that is located non-symmetrically within a spherical beaker. In this case, the Laplace equation $∇2$E(r) = 0 does not have a compact closed-form solution in three dimensions.¹⁶ In two dimensions, however, the corresponding solution for E(x, y) is known.¹⁷ Suppose that the beaker is centered at the origin while the trap is centered at (c, 0) (with c < R − a to ensure that the trap is entirely inside the beaker). In the corresponding electrostatic problem, the potential between the conductors is identical to the potential induced by two equal and opposite point charges that are located at (c + d₁, 0) and (d₂, 0), with

We now let

be the distances between each charge and the diffusing molecule, with the molecule initially at (x, y) ∈ $R$⁠, with

the region between the conductors. The probability that a molecule that starts at (x, y) is absorbed by the beaker wall is

The exit probability is the spatial average of E(x, y) over the region $R$⁠,

This double integral does not seem to be expressible through known functions, however. Thus for the simplest setting of a single circular trap inside a circular beaker, the exit probability has a closed form only when these two objects are concentric.

The Smoluchowski reaction-rate approach (RRA) is an effective-medium approximation that replaces the complicated influences of each individual trap by their average influence (see, e.g., Refs. 8–14). In this approach, the diffusion equation (4) is replaced by the reaction-diffusion equation

where K is the reaction rate that accounts for the average absorption rate of diffusing molecules by the traps. For identical spherical traps, the reaction rate is^8–14 

where n = N/V (with N the number of traps and V the beaker volume) is the trap density. This expression for K is asymptotically exact when the volume fraction of the traps is small, so that their interactions can be neglected. The RRA can be generalized in a straightforward way to cases where the traps are mobile and where the size of the molecules cannot be ignored.

It is easy to compute the reaction rate via the electrostatic analogy (apparently first noticed in the classical paper;¹⁵ see Refs. 11 and 13 for reviews). This approach gives

where C is the electrical capacitance of the trap (assuming that it is a perfect conductor of the same shape). For example, a sphere of radius a has capacitance C = a, so that (14b) reduces to (14a); a disk of radius a has capacitance 2a/π, leading to K = 8Dan. The electrical capacitance of ellipsoids is also known. For instance, for oblate spheroids with semi-axes a = a ≤ b, the reaction rate is

In the following, we treat spherical traps and use (14a); arbitrary geometries can be obtained by straightforward generalization.

Within the RRA, we need to solve the reaction-diffusion equation with an absorbing boundary condition at the beaker wall

for a spatially uniform initial condition ρ(r, t = 0) = 1. Since the molecular flux to the wall is −D$∇$ρ, the total number of molecules that exit the beaker is

where the second integral is over the beaker wall and n is the unit normal to the wall. Dividing by the total initial number of molecules V gives the exit probability

We now compute the escape and trapping probabilities for (i) a spherical beaker in three directions, (ii) a spherical beaker in general dimensions, and finally (iii) other simple beaker shapes in three dimensions. Throughout, we consider the traps to be identical small spheres; however, our results can be extended to traps of other shapes if the electrostatic capacitance of this shape is known.

For spherical traps, the reaction rate is given by (14a). When the beaker is a ball of radius R, the surface integral in (16) is readily performed and the exit probability is

To compute this integral, we must solve the reaction-diffusion equation (13). This solution ( Appendix A) gives the concentration ρ as the Fourier series,

Substituting this expression in (17) and performing the integral, the exit probability is

The second equality has been obtained by standard residue calculus methods.

The limiting behaviors of Eq. (19) are instructive. As λ → 0, we obtain

while in the opposite λ → ∞ limit

The latter approximation is almost indistinguishable from the exact solution (19) for λ$≳$ 3 (see Fig. 3) because the next correction term in (20b) is of the order of $e−12λ$⁠. As we shall see, these two limiting behaviors of T ∼ λ for λ → 0 and E ∼ λ^−1/2 for λ → ∞ arise quite generally.

For general spatial dimensions, we must solve the d-dimensional reaction-diffusion (13). Using spherical symmetry, we assume that the initial condition is a normalized spherical shell of probability at r = r₀,

where Ω_d is the surface area of a d-dimensional unit sphere. We now Laplace transform (13) and re-express the radius in dimensionless units $x=rs/D$⁠, after which the reaction-diffusion equation (13) becomes

Here the tilde denotes the Laplace transform, and the prime denotes differentiation with respect to x. This equation should be solved subject to the boundary condition ρ(R, t) = 0.

From the solution to Eq. (21) (see  Appendix B), the probability E(r₀) that a diffusing molecule that starts at r₀ ultimately escapes to the beaker wall is given by

Here I_ν and K_ν are the modified Bessel functions of order ν, with ν = (2 − d)/2, and $μ=1+λ/s$⁠.

A more meaningful measure of the trapping efficiency of the system is the escape probability when the fluid initially contains a uniform density of diffusing molecules. Here we need to integrate (22) over all radii and also modify (22) for the initial condition of a uniform density in the entire beaker rather than a spherically symmetric density shell at radius r₀. For the latter, we merely need to multiply the result of integrating (22) over all radii by the factor d/R to account for the difference in the normalizations between a uniform density within a sphere and on a spherical surface.

For the spatial integration, we make use of the identity

as well as I_ν = I_−ν to write

In three dimensions, the above expression reduces to that given in (19). From (23), the exit probability has the same asymptotic behaviors (3) in all dimensions d > 2.

We studied the trapping characteristics of diffusing molecules that are absorbed whenever they either reach the boundary of a finite fluid-filled beaker or touch traps within the fluid. The solution for a single trap is simple for the example of a spherical trap at the center of a spherical beaker. It is impractical, however, to generalize this analytical approach for more than a single trap because it involves solving the diffusion equation subject to absorbing boundary conditions on the surface of all traps and the surface of the beaker. However, the problem greatly simplifies when there are many traps. In this limit, we can replace the traps by an effective and uniformly absorbing medium—the reaction rate approximation (RRA). We showed how to compute the fraction of molecules that escape the domain for the situation where the total number of traps is large. Remarkably, this escape probability depends only on a single parameter λ, the product of the total number of traps and the ratio of the characteristic size of the trap to the characteristic size of the domain. Our general result is that the exit probability from the beaker, E, exhibits universal asymptotic behaviors (3) in all dimensions greater than two.

While our model describes the absorption of diffusing molecules in a beaker that contains a reactive fluid, there are other potential scenarios. For example, the domain could be a cell and the diffusing molecules may exit the cell membrane or may be absorbed by traps within the cell and thereby perform some biologically useful function. Here, the exit may be treated as a loss and one would like to estimate the magnitude of the loss. In realistic systems, the parameter λ is usually large, so our results can be interpreted as the generic claim that this loss, viz., the fraction of molecules exiting through the cell membrane, universally scales as $A/λ$ when λ $≫$ 1. Only the amplitude A depends on the shape of the cell and on the shape of the traps.

Within the RRA, we determined the exit and trapping probabilities without needing to solve the underlying reaction-diffusion equation. In addition to simplifying the calculations, the RRA provides additional insights that would be very difficult to obtain by direct means. For example, for a cylindrical beaker of radius R and height H, the exit probability could, in principle, depend on three dimensionless parameters: N, a/R, and R/H. The RRA predicts that for very tall cylinders (R $≪$ H), the exit probability again actually depends on the single parameter λ = Na/H. It is worth noting that within the RRA formalism it is also straightforward to account for diffusing traps as well as non-zero radii for both the traps and the molecules. All that is needed is to amend the parameter λ from Na/R to

where the subscripts m and t refer to the molecules and the traps, respectively.

We thank Sergei Rudchenko, Nagendra K. Panduranga, and Kirill Korolev for helpful discussions. This work was partly supported by the National Science Foundation under Grant No. DMR-1608211.

To solve the reaction-diffusion equation (13) in three dimensions, we use a simplification that transforms the radial Laplacian that operates on the quantity u(r)/r to a one-dimensional Laplacian,

This identity suggests that we work with u = ρ r, after which the reaction-diffusion equation, with K given by (14a), becomes

It is useful to now introduce the dimensionless variables,

and then drop the overbars to simplify notation in the following. Equation (A2) becomes

We want to solve (A4) subject to the condition of a uniform initial density

an absorbing boundary condition at the beaker wall

and the condition

that ensures that ρ = u/r is well-defined in the origin. The exit probability (17) now becomes

We use the Fourier analysis (see, e.g., Ref. 18) to solve Eq. (A4) subject to (A5). Here it is useful to extend the domain of u(r, t) from [0, 1] to [−1, 1]. For convenience, we also extend u to be odd [i.e., u(−r) = −u(r)], so that the boundary condition (A5c) manifestly holds. Then the Fourier expansion of u(r, t) contains only sine terms,

From the initial condition (A5a), we obtain

Substituting (A7) into (A4) gives the differential equation for the amplitudes

Solving this equation, the concentration ρ(r, t) may be written as the infinite series in Eq. (18).

The Laplace transform of the reaction-diffusion equation (21) is

where $x=rs/D$ is the scaled coordinate, and x is in the range (0, X), with $X=Rs/D$⁠. This equation should be solved subject to the absorbing boundary condition at x = X. The left-hand side of (B1) is a Bessel differential equation that we solve separately for x < x₀ and x > x₀ and then patch together these two solutions by the standard joining conditions for Green’s functions. In each subdomain x < x₀ and x > x₀, the elemental solutions have the form

Here I_ν and K_ν are the modified Bessel functions of order ν, with ν = (2 − d)/2, and $μ=1+λ/s$⁠. The subscripts < and > denote the solution in the regions x < x₀ and x > x₀. The interior solution does not contain the function K_ν because K_ν diverges at the origin. Invoking the absorbing boundary condition ρ = 0 at x = X, as well as the continuity of Green’s function at r = r₀, the solution (B2) can be expressed in a form that is manifestly continuous and vanishes at x = X,

The constant A is determined by the joining condition

For differentiating Green’s function, we use

while for $ρ>′−ρ<′$⁠, we use the Wronskian relation

With these identities and performing some tedious but straightforward algebra, the amplitude of Green’s function is given by

Finally, we compute the flux −Dρ′ to the outer boundary, integrate over the surface of the sphere, and then take the s → 0 limit of the Laplace transform to obtain the probability E(r₀) that a diffusing molecule that starts at x₀ ultimately escapes to the beaker wall. This is Eq. (22).

Here we solve (13) and compute the exit probability for two additional cases.

Consider the parallelepiped of size L₁ × L₂ × L₃. The reaction-diffusion equation is

We seek a solution in the form of a Fourier sine series

which ensures that the absorbing boundary condition (15) is automatically satisfied on the walls x_j = 0, L_j. Substituting (C2) into (C1) and solving for F_n(t) give

where we use the shorthand notation

For the spatially uniform initial condition, we find that F_n(0) = 0 if at least one index n_j is even; when all n_j’s are odd, we have

After some straightforward calculations, the exit probability is

Here and below, $∑n$ runs over n₁, n₂, n₃ which are all positive and odd.

Let us consider two examples in more detail. For the cube with L₁ = L₂ = L₃ = L, we obtain

where we use the shorthand $n2=n12+n22+n32$ with λ = 4aN/πL.¹⁹ Using (C3) and the identity²⁰ 

we determine the trapping probability

The small λ expansion of this trapping probability is

where

To determine the large-λ behavior, we use symmetry and re-write (C3) as

Next, we replace the summation over n₁ by integration

from which, we can obtain the leading asymptotic behavior in the λ → ∞ limit. Computing the integral gives

In the λ → ∞ limit, we make the replacement $n22+n32+λ→λ$ and use the identity (C4) to find

As another example, consider the bar with dimensions L₁ = L₂ $≪$ L₃. We perform the summation over n₃ by using identity (C4). Then the exit probability becomes

where the summation runs over n₁, n₂ which are both positive and odd.

We obtain the asymptotic behaviors following similar steps as for the cube. The small-λ expansion of the trapping probability is given by (C6) with

Using (C10) and calculations similar to (C8) and (C9), we obtain

Generally for the L₁ × L₂ × L₃ parallelepiped, the exit probability decays (for N $≫$ L/a) as

From this general result, one can recover (C9) for the cube and (C11) for the bar.

If the cylinder height H greatly exceeds its radius R, the problem becomes two dimensional. In the dimensionless variables of (A3), the reaction-diffusion equation (13) becomes

The solution is the Bessel series

The absorbing boundary condition is satisfied when 0 < μ₁ < μ₂ < ⋯ are consecutive zeros of the Bessel function J₀(μ_n) = 0. For the uniform initial condition, and using the orthogonality condition $∫01drrJ0(μnr)J0(μmr)=0$ for n$≠$m, the coefficients A_n are

from which the exit probability is

Hence E ∝ λ⁻¹ for large λ. In contrast, for three-dimensional systems, the exit probability always scales as λ^−1/2 for λ$≫$ 1.