Trapping diffusive particles at surfaces is a key step in many systems in chemical and biological physics. Trapping often occurs via reactive patches on the surface and/or the particle. The theory of boundary homogenization has been used in many prior works to estimate the effective trapping rate for such a system in the case that either (i) the surface is patchy and the particle is uniformly reactive or (ii) the particle is patchy and the surface is uniformly reactive. In this paper, we estimate the trapping rate for the case that the surface and the particle are both patchy. In particular, the particle diffuses translationally and rotationally and reacts with the surface when a patch on the particle contacts a patch on the surface. We first formulate a stochastic model and derive a five-dimensional partial differential equation describing the reaction time. We then use matched asymptotic analysis to derive the effective trapping rate, assuming that the patches are roughly evenly distributed and occupy a small fraction of the surface and the particle. This trapping rate involves the electrostatic capacitance of a four-dimensional duocylinder, which we compute using a kinetic Monte Carlo algorithm. We further use Brownian local time theory to derive a simple heuristic estimate of the trapping rate and show that it is remarkably close to the asymptotic estimate. Finally, we develop a kinetic Monte Carlo algorithm to simulate the full stochastic system and then use these simulations to confirm the accuracy of our trapping rate estimates and homogenization theory.

Many systems in chemical and biological physics involve diffusive particles being trapped at surfaces. Trapping often occurs via localized “patches” (or “traps,” “binding sites,” “pores,” etc.) on the surface, meaning a particle reacts if it contacts a patch but reflects from the rest of the surface. Examples include proteins binding to receptors on a cell membrane;1 industrial processes, such as filtration2 and gas separation;3 reactions on porous catalyst support structures;4 diffusion current to collections of microelectrodes;5 and water transpiration through plant stomata.6,7

Mathematical models of trapping by such “patchy” surfaces typically use the diffusion equation to describe particle motion with absorbing boundary conditions at patches and reflecting boundary conditions on the rest of the surface (imposing absorbing boundary conditions is a common approximation, but it has recently come under scrutiny8–10). These mixed boundary conditions are often replaced by a single effective “trapping rate” κ > 0 in a technique called “boundary homogenization.”11,12 This trapping rate (also called the reactivity, leakage parameter, homogenization constant, reaction constant, or radiation constant) encapsulates the effective trapping properties of the patchy surface.

To illustrate boundary homogenization of a patchy surface [see Fig. 1(a)], let c(x, y, z, t) denote the concentration of diffusing particles at time t, where x ≥ 0 is the distance from the surface and (y,z)R2 is the location in the plane parallel to the surface. The concentration evolves according to the diffusion equation,
tc=Dtr2x2+2y2+2z2c,
where Dtr > 0 denotes the translational diffusivity of the particle. Trapping at the patchy surface is described by mixed boundary conditions,
c=0,x=0,(y,z) in a surface patch,xc=0,x=0,(y,z) not in a surface patch.
(1)
Boundary homogenization replaces (1) by the following partially absorbing (also called partially reactive, Robin, or third type) boundary condition involving the effective trapping rate κ1 > 0:
Dtrxc=κ1c,x=0.
In the case that the patches are roughly evenly distributed with common radius r1 > 0 and occupy a small fraction
f1N1πr124πR21
(2)
of the surface, the trapping rate is (see Refs. 1, 11, 13, and 14 or  Appendix A) given by
κ14Dtrπr1f1.
(3)
In (2), N1 > 0 is the average number of patches in an area of size 4πR2 for some length scale R > 0 [i.e., N1/(4πR2) is the density of patches]. Since the seminal derivation of the leading order result in (3), a large literature has been devoted to obtaining higher order corrections, which account for details, such as patch shape, patch arrangement and clustering, surface curvature, and patch diffusivity.11–21 
FIG. 1.

Schematic diagram of a diffusive particle trapping at a surface. (a) A perfectly reactive particle (red sphere) that is trapped by small patches on the surface (blue disks). (b) A patchy particle that is trapped when small patches on the particle (red disks) contact a perfectly reactive surface (blue wall). (c) A patchy particle that is trapped when small patches on the particle (red disks) contact small patches on the surface (blue disks).

FIG. 1.

Schematic diagram of a diffusive particle trapping at a surface. (a) A perfectly reactive particle (red sphere) that is trapped by small patches on the surface (blue disks). (b) A patchy particle that is trapped when small patches on the particle (red disks) contact a perfectly reactive surface (blue wall). (c) A patchy particle that is trapped when small patches on the particle (red disks) contact small patches on the surface (blue disks).

Close modal
Other systems involve so-called patchy particles in which trapping occurs via localized patches on the particles.22–26 To illustrate boundary homogenization of a patchy particle [see Fig. 1(b)], consider a concentration c(x, θ, φ, t) of diffusing spherical particles at time t, where x ≥ 0 is the distance from the surface and (θ, φ) ∈ [0, π] × [0, 2π) describes the particle orientation in spherical coordinates. The concentration evolves according to the diffusion equation,
tc=Dtr2x2c+Drot(sinθ)22φ2+(tanθ)1θ+2θ2c,
where Dtr > 0 and Drot > 0 denote the respective translational and rotational diffusivities. If trapping only requires that a patch on the particle contacts the surface (i.e., the surface is fully reactive), then the following mixed boundary conditions are imposed:
c=0,x=0,(θ,φ) in a particle patch,xc=0,x=0,(θ,φ) not in a particle patch.
(4)
If each particle is a sphere of radius R > 0 containing N2 ≫ 1 roughly evenly distributed patches of common radius r2R occupying a small fraction
f2N2πr224πR21
(5)
of the particle’s boundary, then (4) can be approximated by25 
Dtrxc=κ2c,x=0,
where the trapping rate is
κ2β4Dtrπr2f2,
(6)
where
βR2DrotDtr>0
compares the rotational diffusivity to the translational diffusivity. We note that if the diffusion of the particle satisfies the Stokes–Einstein relation, then β=3/40.87.
In this paper, we combine the two aforementioned scenarios and consider patchy surfaces trapping patchy particles [see Fig. 1(c)]. In particular, we suppose that the patchy particles are trapped by the patchy surface only when particle patches contact surface patches. We use matched asymptotic analysis to derive the following trapping rate:
κχf1f2κ1κ2,
(7)
where
χ=χ(βr1/r2)>0
is a dimensionless function of the ratio βr1/r2. In particular, χ depends on the electrostatic capacitance of a certain four-dimensional object called a duocylinder27 embedded in five-dimensional space. While we do not have an exact formula for χ, we employ a kinetic Monte Carlo algorithm to rapidly compute χ to very high accuracy.
In addition, we use the theory of Brownian local time28 to derive the following heuristic trapping rate:
κ0f1κ2+f2κ1.
(8)
Using the particular forms in (2), (3), (5), and (6), the heuristic formula in (8) is equivalent to
κ0=χ0f1f2κ1κ2,
where
χ0=βr1/r2+1βr1/r2/π.
We find that the heuristic trapping rate κ0 is a remarkably accurate approximation to the asymptotic trapping rate κ in (7). Indeed, the accuracy depends on the value of βr1/r2, but the relative error is never more than 16%,
κκ0κ0.16.

Finally, we develop a kinetic Monte Carlo algorithm to simulate the full stochastic system. This algorithm partitions the stochastic path of the particle into simpler steps involving radially symmetric geometries, which are amenable to analytical and numerical solution. Comparing these simulation results to the theory confirms the accuracy of boundary homogenization using the trapping rate in (7).

The rest of this paper is organized as follows. In Sec. II, we formulate the mathematical model precisely in terms of a stochastic process and an associated partial differential equation (PDE) in five spatial dimensions. In Sec. III, we apply matched asymptotic analysis to this PDE to derive the trapping rate (7) in the limit of small patches. Since this trapping rate involves the capacitance of a four-dimensional duocylinder, we employ a kinetic Monte Carlo method to rapidly compute this capacitance in Sec. IV. In Sec. V, we derive the heuristic trapping rate in (8). In Sec. VI, we compare the results of our analysis to stochastic simulations of the full problem. We conclude by discussing related work and highlighting future directions.

Consider a spherical particle with radius R > 0 diffusing in a slab of finite or infinite width 0 < L + 2R ≤ ∞. Suppose that the left wall of the slab contains infinitely many reactive patches with average density N1/(4πR2) > 0 and the particle contains N2 ≥ 1 reactive patches on its surface. Suppose that the patches on the wall are disks of radius r1 > 0 and the patches on the particle are spherical caps with polar angle r2/R > 0.

At time t ≥ 0, let X(t) ∈ [0, L] denote the distance between the particle and the left wall and let (Y(t),Z(t))R2 denote the position of the center of the particle in the plane parallel to the left wall. Let (Θ(t), Φ(t)) ∈ [0, π] × [0, 2π) describe the orientation of the particle in spherical coordinates. The state of the system at time t ≥ 0 is thus fully described by the five coordinates,
X(t),Y(t),Z(t),Θ(t),Φ(t)Ω̄,
where Ω̄ denotes the closure of
Ω(0,L)×R2×[0,π]×[0,2π).
Suppose that the particle diffuses with translational diffusivity Dtr > 0, which means that the coordinates (X(t), Y(t), Z(t)) satisfy the following stochastic differential equations (SDEs):
dX=2DtrdWX,X(0)=x0,dY=2DtrdWY,Y(0)=yR,dZ=2DtrdWZ,Z(0)=zR,
where WX, WY, and WZ are independent, standard Brownian motions. The particle reflects from the left wall of the slab and also from the right wall if L < ∞. Suppose that the particle has rotational diffusivity Drot > 0, which means its orientation (Θ(t), Φ(t)) satisfies the following SDEs:
dΘ=DrottanΘ(t)dt+2DrotdWΘ,Θ(0)=θ(0,π),dΦ=2DrotsinΘ(t)dWΦ,Φ(0)=φ[0,2π),
where WΘ and WΦ are independent, standard Brownian motions.
We are interested in the first time that one of the patches on the particle contacts one of the patches on the left wall, which we refer to as the reaction time or trapping time τ. To describe the patches on the left wall, let γ1(y0, z0, r) denote the disk of radius r > 0 centered at (y0,z0)R2,
γ1(y0,z0,r){(y,z):(yy0)2+(zz0)2<r2}.
If {(ym,zm)}m=1 denote the centers of the infinitely many patches on the left wall, then the set of all patches on the left wall is
Γ1m=1γ1(ym,zm,r1).
To describe the patches on the particle, let γ2(θ0, φ0, a) denote the spherical cap with polar angle a > 0 centered at (θ0, φ0),
γ2(θ0,φ0,a){(θ,φ):(θθ0)2+sin2(θ0)(φφ0)2<a2}.
If {(θn,φn)}n=1N denote the centers of the N2 ≥ 1 patches on the particle, then the set of all patches on the particle is
Γ2n=1N2γ2(θn,φn,r2/R).
The reaction time is then given by
τinft>0:X(t)=0,(Y(t),Z(t))Γ1,(Θ(t),Φ(t))Γ2.
(9)
In words, the reaction time τ in (9) is the first time that the following three events occur simultaneously: (i) the particle is touching the wall (X = 0), (ii) the particle’s planar coordinates are in a patch on the wall ((Y, Z) ∈ Γ1), and (iii) the particle’s angular orientation coordinates are in a patch on the particle ((Θ, Φ) ∈ Γ2).
The probability distribution of the reaction time τ in (9) is described by the survival probability conditioned on the initial state of the system,
S(x,y,z,θ,φ,t)Pτ>t|(X(0),Y(0),Z(0),Θ(0),Φ(0))=(x,y,z,θ,φ).
(10)
The survival probability satisfies the backward Fokker–Planck (backward Kolmogorov) equation,29 
tS=LS,(x,y,z,θ,φ)Ω,
(11)
where L=Ltr+Lrot and Ltr is the Laplacian acting on (x, y, z),
Ltr=Dtr2x2+2y2+2z2,
and Lrot is the Laplace–Beltrami operator acting on (θ, φ),
Lrot=Drot(sin(θ))22φ2+(tan(θ))1θ+2θ2.
The survival probability satisfies a unit initial condition,
S=1,t=0,
(12)
and either a reflecting condition at x = L if L < ∞,
xS=0ifx=L<,
(13)
or a unit far-field condition if L = ∞,
limxS=1ifL=.
(14)
Condition (13) encapsulates the assumption that the particle reflects from the right wall if L < ∞. Condition (14) means that a particle cannot react before any fixed time t > 0 if it started infinitely far from the left wall. The survival probability satisfies mixed boundary conditions at x = 0,
S=0,x=0,(y,z),(θ,φ)Γ1×Γ2,xS=0,x=0,(y,z),(θ,φ)Γ1×Γ2.
(15)
Boundary condition (15) reflects the assumption that the reaction only occurs if a patch on the particle contacts a patch on the left wall.
We now derive the effective trapping rate κ, assuming that the patches (i) are roughly evenly distributed and (ii) occupy a small surface area fraction of the wall (f1 ≪ 1) and the particle (f2 ≪ 1) (these two assumption imply that the patches are well-separated on each surface in the sense that the patch radius is much smaller than the typical distance between patches). Specifically, we derive the appropriate κ such that the solution S(x, y, z, θ, φ, t) to the five-dimensional PDE in (11)(15) is well approximated by the solution S̄(x,t) to the following one-dimensional diffusion equation:
tS̄=Dtr2x2S̄,x(0,L),
(16)
where S̄ satisfies (12)(14) and the following Robin boundary condition at x = 0:
DtrxS̄=κS̄,x=0.
(17)
Recall that the radius of the patches on the wall is r1 > 0 and the patches on the particle are spherical caps with polar angle r2/R > 0. Assume that
r1=εa1R,r2=εa2R,
(18)
where 0 < ɛ ≪ 1 and a1 > 0, a2 > 0 are order one constants, which allow the patches on the wall to differ in size from the patches on the particle (without loss of generality, we can take either a1 = 1 or a2 = 1, but we keep both parameters for pedagogical reasons). To derive the trapping rate κ, we apply matched asymptotic analysis19,25 to the five-dimensional PDE in (11)(15) in the limit ɛ → 0. In  Appendix A, we use this same asymptotic approach to derive (3).
Fix any time t > 0. As ɛ → 0, the surfaces of the particle and the left wall become perfectly reflecting, and thus,
limε0S=1.
(19)
We expect that S has a boundary layer in a neighborhood of each possible intersection of patches, so we introduce the outer expansion
S1+ε3S1+,
(20)
which is valid away from the patches. The leading order correction of ɛ3 in (20) can be derived by considering the case of zero rotational diffusion. In this case, reaction can only occur if a patch on the particle is initially oriented toward the wall, which has probability f1 = O(ɛ2) (assuming an initial uniform distribution of particle orientation). Given that a patch on the particle is oriented toward the wall (and with zero rotational diffusion), the problem is equivalent to a perfectly reactive particle and a patchy surface, which has a leading order correction of O(ɛ) (see  Appendix A). Combining these observations yields the leading order correction of O(ɛ3) in (20). We comment that the expansion in (20) assumes that t is a fixed time independent of ɛ [in particular, (20) is not valid if we take t → ∞ for fixed ɛ].
Plugging outer expansion (20) into (11) shows that S1 must satisfy
tS1=LS1,(x,y,z,θ,φ)Ω,xS1=0,x=0,(y,z,θ,φ)m=1(ym,zm)×n=1N2(θn,φn).
(21)
Observe that from the perspective of the outer solution, the patches have shrunk to points.
We now determine the behavior of S1 as the nth trap on the particle approaches the mth trap on the left wall,
(x,y,z,θ,φ)(0,ym,zm,θn,φn).
Toward this end, we introduce the following stretched coordinates for fixed m ∈ {1, 2, …} and n ∈ {1, …, N2}:
μ=ε1x/R,ν=ε1(yym)/R,η=ε1(zzm)/R,s=ε1(θθn)/β,p=ε1sin(θn)(φφn)/β,
where
βR2Drot/Dtr>0
compares the rotational diffusivity to the translational diffusivity. We then define the inner solution w as a function of the stretched coordinates,
w(μ,ν,η,s,p,t)SεRμ,ym+εRν,zm+εRη,θn+εβs,φn+εβsinθnp,t.
(22)
Using (11), the inner solution w satisfies
tw=Lw=ε2R2Dtr(μμ+νν+ηη)w+ε2Drotsin(θn)2β2pp+β2ssw+O(ε1)=ε2R2Dtr(μμ+νν+ηη+ss+pp)w+O(ε1)
(23)
since β is defined so that Drotβ−2 = R−2Dtr. Therefore, plugging the inner expansion
w=w0+O(ε)
into (23) implies that w0 is harmonic,
(μμ+νν+ηη+ss+pp)w0=0,μ>0.
(24)
Using (15), w0 satisfies the following mixed boundary conditions at μ = 0:
w0=0,μ=0,ν2+η2<b12,s2+p2<b22,μw0=0,μ=0,otherwise,
(25)
where b1 = a1 and b2 = a2/β.
It follows from electrostatics30 that w0 has the far-field behavior,
w0α1cρ3 as ρμ2+ν2+η2+s2+p2,
(26)
where α is a constant determined by matching with the outer solution and
cc0(a1,a2/β)>0
is a dimensionless constant depending on a1 and a2/β. In particular, c0(b1, b2) is the electrostatic capacitance of a so-called duocylinder27 embedded in R5,
D(0,ν,η,s,p)R5:ν2+η2<b12,s2+p2<b22.
(27)
We postpone the discussion of computing c0(a1, a2/β) until Sec. IV and proceed in terms of c = c0(a1, a2/β).
The matching condition is that the near-field behavior of the outer expansion as (x, y, z, θ, φ) → (0, ym, zm, θn, φn) must agree with the far-field behavior of the inner expansion as ρ → ∞. That is,
1+ε3S1w0as (x,y,z,θ,φ)(0,ym,zm,θn,φn),ρ.
Using (22) and (26), it follows that α = 1 and that S1 has the following singular behavior as (x, y, z, θ, φ) → (0, ym, zm, θn, φn):
cxR2+yymR2+zzmR2+θθnβ2+sin2(θn)φφnβ232.
(28)
Writing the singular behavior (28) in the distributional form for each n ∈ {1, …, N2} and m ∈ {1, 2, …}, the boundary condition in (21) at x = 0 becomes
xS1=Km=1n=1N2δ(yym)δ(zzm)δ(θθn)sin(θn)δ(φφn),
(29)
where K = 4π2cRβ2. The derivation of the distributional form in (29) of the singular behavior in (28) is given in  Appendix B.
We have assumed that the patches on the wall have average density N1/(4πR2), which means that
liml1(2l)2[l,l]2m=1δ(yym)δ(zzm)dydz=N14πR2
(30)
for some N1 ∈ (0, ∞). The integral in (30) is simply the number of patches in the square [−l, l]2, and thus, N1 is the average number of patches on the wall in an area equal to the surface area of the particle (note that N1 need not be an integer).
Define
S̄(x,t)liml14π(2l)2ΩlS(x,y,z,θ,φ,t)dΣ,
where
Ωl[l,l]2×[0,π]×[0,2π),dΣsinθdθdφdydz.
Furthermore, define the ratio
κ̄DtrxS̄(0,t)S̄(0,t)
(31)
so that it is a tautology that S̄ satisfies
DtrxS̄=κ̄S̄,x=0.
We now derive the effective trapping rate by determining the behavior of (31) as ɛ → 0.
The denominator in (31) approaches unity as ɛ → 0 due to (19). To determine the behavior of the numerator in (31), we interchange the derivative with the limit and integrals, recall the outer expansion in (20), and apply the boundary condition in (29),
κ̄Dtrliml14π(2l)2ΩlxS(0,y,z,θ,φ,t)dΣε3N1N2DtrRβ2c4as ε0.
(32)
Now, a simple scaling argument shows that
c0(b1,b2)=α3c0(b1/α,b2/α)for all α>0.
(33)
Therefore, recalling c = c0(a1, a2/β), we may rewrite (A11) as
κ̄ε3N1N2DtrRβ(a1a2)3/24c0r1βr2,r2r1βas ε0.
(34)
To express (34) in a more intuitive form, note first that the fractions of the wall and particle covered by patches are
fiNiε2ai24=Niri2R2,i{1,2}.
(35)
Furthermore, if (i) κ1 denotes the trapping rate in the case that the wall is patchy and the particle is perfectly reactive and (ii) κ2 denotes the trapping rate in the case that the particle is patchy and the wall is perfectly reactive,
κ14Dtrπr1f1,κ2β4Dtrπr2f2,
(36)
then (34) becomes
κ̄κχf1f2κ1κ2as ε0,
(37)
where χ is the following function of r1β/r2:
χ=χ(r1β/r2)=πc0(r1β/r2,1/r1β/r2).
(38)

The homogenized trapping rate κ in (37) involves the factor χ in (38), which depends on the electrostatic capacitance c0 of the four-dimensional duocylinder D in (27) embedded in R5. We calculate this capacitance by modifying the kinetic Monte Carlo method developed in Ref. 26 to calculate the capacitance of a different four-dimensional region embedded in R5. This previous method adapted a method devised by Bernoff et al.12 for related trapping problems in R3. These kinetic Monte Carlo methods modify the 1956 “walk-on-spheres” method by Muller.31 

We now summarize this calculation of c0. The details are given in  Appendix C. The leading order inner solution w0 satisfying (24) and (25) can be written in terms of a certain “splitting probability,” which is the probability that a Brownian particle diffusing in five-dimensional space eventually reaches the duocylinder D in (27). Since the capacitance c0 is determined by the far-field behavior of w0 [see (26)], we calculate c0 from this splitting probability.

We estimate this splitting probability by simulating the diffusive paths of M ≫ 1 independent Brownian particles and counting the fraction that reaches D before reaching some large distance ρ ∈ (0, ∞) from D. Simulating these stochastic paths can be done very efficiently by using a kinetic Monte Carlo algorithm. The only error in this calculation of the capacitance c0 is due to (i) the finite number of trials M and (ii) the finite outer radius ρ. Due to the computational efficiency of the kinetic Monte Carlo algorithm, M and ρ can both be taken very large, which ensures that the error in calculating c0 is very small.

We now present the results of calculating c0. It is immediate that the capacitance is symmetric,
c0(b1,b2)=c0(b2,b1)for all b1>0,b2>0,
(39)
and thus, by the scaling in (33), it is enough to compute
c0(1,b),b(0,1].
(40)
In Fig. 2, we plot the capacitance in (40) as a function of b ∈ (0, 1]. The black solid curve in Fig. 2 is from the kinetic Monte Carlo algorithm described above with outer radius ρ = 105 and M = 107 trials. As mentioned above, the only sources of error in this computation are the finite outer radius and the finite number of trials. As detailed in Ref. 26, the relative error stemming from the finite outer radius vanishes like ρ3 as ρ → ∞. Hence, setting ρ = 105 means that relative error in calculating the capacitance from the finite outer radius is on the order of 10−15.
FIG. 2.

Capacitance in (40) of duocylinder D in (27) with radii b1 = 1 and b2 = b ∈ (0, 1]. The black solid curve is the result of the kinetic Monte Carlo simulations described in Sec. IV, and the red dashed curve is the heuristic estimate in (44). The black circles are the heuristic estimate in (45), which is merely the heuristic estimate in (44) scaled so that it agrees with the Monte Carlo simulations when b = 1.

FIG. 2.

Capacitance in (40) of duocylinder D in (27) with radii b1 = 1 and b2 = b ∈ (0, 1]. The black solid curve is the result of the kinetic Monte Carlo simulations described in Sec. IV, and the red dashed curve is the heuristic estimate in (44). The black circles are the heuristic estimate in (45), which is merely the heuristic estimate in (44) scaled so that it agrees with the Monte Carlo simulations when b = 1.

Close modal

Using the confidence intervals from Ref. 26, we find that for each choice of b, the kinetic Monte Carlo estimate of the capacitance in Fig. 2 has a relative error of less than 1.1% with probability 0.95. For choices of b ≥ 0.5, the relative error is less than 0.13% with probability 0.95.

We now give a heuristic derivation of the effective trapping rate using the theory of Brownian local time.28 We then show that the resulting heuristic trapping rate, denoted as κ0, differs from the asymptotic trapping rate, κ in (37), by at most 16% in the small patch limit (ɛ → 0).

If we only homogenize the plane and use the trapping rate κ1 in (3), then the survival probability S1 = S1(x, θ, φ, t) satisfies
tS1=(Dtr2x2+Lrot)S1,DtrxS1=κ1S1,x=0,(θ,φ)Γ2,0,x=0,(θ,φ)Γ2
and (12)(14). In this case, S1 has the following probabilistic representation:28 
S1(x,θ,φ,t)=P(τ1>t|X(0)=x,Θ(0)=θ,Φ(0)=φ),
where
τ1inf{t>0:X,Θ,Φ(t)>E1},
where X,Θ,Φ(t) is the following local time:
X,Θ,Φ(t)limη01η0t1|X(s)|<η1(Θ(s),Φ(s))n(θn,φn)ds,
where 1A denotes the indicator function on an event A (i.e., 1A = 1 if A occurs and 1A = 0 otherwise) and E1 is an independent exponential random variable with rate κ1.
Alternatively, if we only homogenize the particle and use the trapping rate κ2 in (6), then the survival probability S2 = S2(x, y, z, t) satisfies
tS2=LtrS2,DtrxS2=κ2S2,x=0,(y,z)Γ1,0,x=0,(y,z)Γ1
and (12)(14). In this case, S2 has the following probabilistic representation:
S2(x,y,z,t)=P(τ2>t|X(0)=x,Y(0)=y,Z(0)=z),
where
τ2inf{t>0:X,Y,Z(t)>E2},
where X,Y,Z(t) is the following local time:
X,Y,Z(t)limη01η0t1|X(s)|<η1(Y(s),Z(s))m(ym,zm)ds,
and E2 is an independent exponential random variable with rate κ2.
To homogenize both the plane and the particle, we consider the following reaction time:
τ0min{τ1,τ2}=inf{t>0:X,Θ,Φ(t)>E1 or X,Y,Z(t)>E2}.
The survival probability of τ0 depends on the initial angular position (Θ(0), Φ(0)) and the initial planar position (Y(0), Z(0)) (since τ0 depends on X,Θ,Φ and X,Y,Z). However, averaging over the initial angular position (Θ(0), Φ(0)) yields
X,Θ,Φ(t)Θ(0),Φ(0)=f2(t),
(41)
where (t) is the local time of X at x = 0,
(t)=limη01η0t1|X(s)|<ηds.
Similarly, averaging over the initial planar position (Y(0), Z(0)) yields
X,Y,Z(t)Y(0),Z(0)=f1(t).
(42)
Therefore, we consider a reaction time τ0̄, which is identical to τ0 in (9) except that we replace X,Θ,Φ and X,Y,Z by their averaged versions in (41) and (42),
τ0̄=inf{t>0:f2(t)>E1 or f1(t)>E2}=inf{t>0:(t)>min{E1/f2,E2/f1}}.
Now, it follows from basic properties of exponential random variables that E1/f2 is exponentially distributed with rate f2κ1 (since E1 is exponentially distributed with rate κ1). Similarly, E2/f1 is exponentially distributed with rate f1κ2 (since E2 is exponentially distributed with rate κ2). Furthermore, the minimum of two exponential random variables with rates f1κ2 and f2κ1 is exponentially distributed with rate f1κ2 + f2κ1. Therefore, τ0̄ is equal in distribution to
τ0̄=distinf{t>0:(t)>E0},
where E0 is exponentially distributed with rate
κ0f1κ2+f2κ1.
(43)
Therefore, the survival probability of τ0̄,
S0(x,t)P(τ0̄>t|X(0)=x),
satisfies (12)(14), (16) with trapping rate κ0 in (43) at x = 0,
DtrxS0=κ0S0,x=0.
We now compare the heuristic trapping rate κ0 in (43) with the asymptotic trapping rate κ in (37). Using (35) and (36), we can write κ0 in the following form:
κ0=f1κ2+f2κ1=χ0f1f2κ1κ2
if we define
χ0=χ0(βr1/r2)=βr1/r2+1/βr1/r2.
Furthermore, we have that κ0κ if the following formula approximates the capacitance of the duocylinder:
c0(b1,b2)c̄(b1,b2)b1b2(b1+b2)π.
(44)
Since c̄ satisfies scaling (33) and the symmetry in (39), the approximation in (44) is equivalent to the approximation c0(1,b)c̄(1,b)=b(1+b)/π.
In Fig. 2, the black solid curve is c0(1, b) computed from kinetic Monte Carlo simulations (see Sec. IV) and the red dashed curve is c̄(1,b)=b(1+b)/π. Clearly, c0c̄, but the approximation is reasonably accurate. Indeed, the relative error |c0c̄|/c0 is at most 16%, and thus, the relative error between κ0 and κ is at most 16%,
κκ0κ0.16.
We can obtain an alternative approximation by simply multiplying c̄(1,b) by a constant ζ so that the resulting approximation agrees with c0(1, b) at b = 1. Specifically, the black circles in Fig. 2 are
ζc̄(1,b),where ζ=c0(1,1)/c̄(1,b)1.19.
(45)
The approximation in (45) is more accurate than c̄(1,b) for b ∈ (0.1, 1].

To validate the boundary homogenization theory, we now develop a kinetic Monte Carlo algorithm for simulating a diffusing, patchy particle (sphere) interacting with a patchy surface (plane). The algorithm adapts a kinetic Monte Carlo algorithm developed by Bernoff et al.12 for simulating a diffusing, fully-reactive, spherical particle interacting with a patchy plane.

Our algorithm requires simulating both translational and rotational diffusion of the spherical particle. For these simulations, the patches on the plane are placed on a periodic grid. The algorithm consists of two stages.

  • Stage I: Project from the bulk to the plane. The sphere starts away from the plane, where it is free to diffuse. The sphere is propagated to the location where it first touches the plane. Both this location and the time it takes for the sphere to diffuse to the plane are drawn from exact probability distributions. If the sphere and plane touch such that a patch on the sphere is touching a patch on the plane, the elapsed time of the simulation is recorded and the simulation ends. Otherwise, the algorithm proceeds to stage II.

  • Stage II: Project from the plane to the bulk. The goal of this stage is to propagate the sphere from the plane into the bulk the maximum distance it can be propagated without the possibility of a reaction occurring during the propagation. Stage II is split into two cases:

    • Case A: Reflect from the plane. If at the end of stage I, the location where the sphere touches the plane is not within a patch on the plane, then the shortest distance dA to the set of patches on the plane is calculated. The sphere is then propagated into the bulk to a hemisphere of radius dA. The time for the sphere to diffuse is drawn from an exact probability distribution and recorded, and the algorithm returns to stage I.

    • Case B: Reflect only from the sphere. If at the end of stage I, the location where the sphere touches the plane is within a patch on the plane but the sphere is not oriented such that a patch on the sphere is touching a patch on the plane, then the arc length dB on the sphere from the location of contact to the nearest patch on the sphere is calculated. The spherical orientation angles are then propagated to a distance of dB, and the time for it to do so is drawn from an exact probability distribution. The sphere location is propagated into the bulk according to an exact distribution for diffusion during that time. The algorithm then returns to stage I.

The numerical implementation of these two stages of the algorithm is described in  Appendix D.

We now calculate a numerically optimized trapping rate25  κopt from the stochastic realizations of the reaction times obtained from the algorithm described above. Specifically, given M ≫ 1 stochastic realizations of the reaction time, we order them t1 < t2 <⋯< tM and define the empirical survival probability,
Semp(tj)=1j12M,j=1,,M,
which counts the fraction of particles arriving after time tj.12 We then choose κopt so that the solution S̄ to the homogenized, one-dimensional problem in (16) and (17) approximates the empirical survival probability Semp.
More precisely, for any trapping rate κ, the solution S̄ to (12), (14), (16), and (17) is
S̄(x0,t;κ)=1erfcx04Dtrt+ex024Dtrterfcx2κt+x04Dtrt,
(46)
where erfc(u) denotes the complementary error function and erfcx(u) = exp(u2)erfc(u) denotes the scaled complementary error function. For any trapping rate κ, we use the Kolmogorov–Smirnov distance to calculate the error in distribution between S̄ and Semp,
E(κ)=maxj=1,,N|S(x0,tj;κ)Semp(tj)|.
(47)
The numerically optimized κopt is then the value of the trapping rate, which minimizes E.
We generate M = 106 stochastic realizations of the reaction times using the stochastic simulation algorithm described above for each of the following parameter combinations: r1 = r2 = ɛ ∈ {0.05, 0.10, …, 0.20}, N2 ∈ {11, 55, 99}, N1 = 4π, and Dtr = Drot = 1. Recalling (35), the absorbing fractions thus range from f1 ≈ 0.008 to f1 ≈ 0.126 and f2 ≈ 0.007 to f2 = 0.99. In the top panel of Fig. 3, we plot the numerically optimized κopt values and the theoretical κ values as a function of ɛ for the three values of N2. The theoretical values of κ are calculated from (35)(37) using the numerically-calculated value of c0 from Sec. IV. In the bottom panel of Fig. 3, we plot the relative error,
|κoptκ|κopt,
as a function of ɛ for the three values of N2. This plot shows excellent agreement between theoretical and numerically optimized trapping rate values.
FIG. 3.

Comparing the theoretical trapping rate values to the numerically optimized trapping rates. See the text for details.

FIG. 3.

Comparing the theoretical trapping rate values to the numerically optimized trapping rates. See the text for details.

Close modal

In the top panel of Fig. 4, we plot the Kolmogorov–Smirnov distance E(κ) in (47) (error in distribution) for these theoretical κ values as a function of ɛ for the three values of N2. This plot shows that the error in distribution is quite small for all these parameter values, which validates the boundary homogenization theory. To illustrate this excellent agreement, in the bottom panel of Fig. 4, we plot the empirical survival probabilities Semp and the theoretical survival probabilities S(x0, tκ) in (46) using the theoretical κ values for ɛ = 0.05 and the three values of N2.

FIG. 4.

Comparing the reaction time distribution from simulation to the boundary homogenization theory. See the text for details.

FIG. 4.

Comparing the reaction time distribution from simulation to the boundary homogenization theory. See the text for details.

Close modal
Our homogenization theory was in the limit ɛ → 0, and thus, assumed that the fraction of the surface covered in patches, f1, and the fraction of the particle covered in patches, f2, are both small,
f11,f21.
(48)
However, the kinetic Monte Carlo algorithm for simulating the full stochastic system allows us to investigate the accuracy of the homogenization theory away from regime (48).

In Fig. 5, we fix r1 = r2 = ɛ = 0.1 and let N1, N2 each take the values {11, 55, 99, 199, 299, 399} so that [recalling (35)] f1, f2 each range from 2.75% up to 99.75% (we take Dtr = Drot = 1). In the top panel of Fig. 5, the curves are the theoretical κ computed from (35)(37) and the circles are the numerically optimized values κopt. As expected, κ and κopt agree in regime (48), but the values of κ and κopt deviate away from this parameter regime. In the bottom panel of Fig. 5, we plot the Kolmogorov–Smirnov distance E(κ) in (47) (error in distribution) for the theoretical κ values.

FIG. 5.

Comparing the reaction time distribution from simulation to boundary homogenization theory away from the f1 ≪ 1, f2 ≪ 1 regime. See the text for details.

FIG. 5.

Comparing the reaction time distribution from simulation to boundary homogenization theory away from the f1 ≪ 1, f2 ≪ 1 regime. See the text for details.

Close modal
Interestingly, the bottom panel of Fig. 5 shows that the accuracy of the homogenization theory as measured by the error in distribution (as opposed to the difference between κ and κopt) is fairly high even away from regime (48). The error in distribution is largest in the regimes
f11,f21
(49)
and
f11,f21.
(50)
However, regimes (49) and (50) merely reduce to previously studied problems in homogenization. Indeed, the trapping rate for (49) reduces to the trapping rate for a patchy particle and a uniformly reactive surface [κ2 in (6)]. Similarly, the trapping rate for (50) reduces to the trapping rate for a patchy surface and a uniformly reactive particle [κ1 in (3)].

Summarizing, although our theory assumes that f1 and f2 are both small (and simulations confirm its high accuracy in this regime), the theory is fairly accurate even in the intermediate regimes of (i) f1 ≪ 1 and f2 ≪̸ 1, f2 ≉ 1, (ii) f2 ≪ 1 and f1 ≪̸ 1, f1 ≉ 1, and (iii) f1 ≪̸ 1, f1 ≉ 1 and f2 ≪̸ 1, f2 ≉ 1. Furthermore, the regimes of (49) and (50) reduce to previously studied problems.

In this paper, we studied the trapping of a diffusive particle by a surface, assuming that trapping requires reactive patches on the particle to contact reactive patches on the surface. We formulated a stochastic model of this trapping time and described its probability distribution via a PDE in five spatial dimensions. Assuming that the surface area fractions covered by patches are small, we applied the method of matched asymptotic analysis to this PDE to derive the effective trapping rate κ. This trapping rate κ depends on the electrostatic capacitance of a four-dimensional duocylinder embedded in five dimensions, and we employed a kinetic Monte Carlo simulation algorithm to rapidly and accurately compute this capacitance. Using the theory of Brownian local time,28 we derived a simple heuristic trapping rate κ0 and showed that it is always within 16% of κ in the small patch limit. We further developed a kinetic Monte Carlo algorithm to simulate the full, five-dimensional stochastic model of the trapping time, which verified the accuracy of the effective trapping rate κ.

This work is related to a long line of previous studies. Boundary homogenization has been used to study the reaction kinetics for patchy particles18–21,25 and patchy surfaces,11,12,15,16,20,32,33 where the patchy object (particle or surface) interacts with a uniformly reactive particle or surface. Related work on interactions of two patchy spherical particles includes Refs. 24, 34, and 35. In Ref. 34, Monte Carlo simulations were used to study the impacts of translational and rotational diffusion on interactions of two or more patchy particles. A computational approach to association and dissociation of patchy particles was developed in Ref. 24. The authors of Ref. 35 used lattice and lattice-adjacent models to study the interactions of pairs of patchy particles. We are not aware of prior work on the interactions between a patchy particle and a patchy surface. Mathematically, our use of the method of matched asymptotic analysis follows a similar method used in Refs. 19 and 26 and also in Refs. 21, 25, 3642. This approach is related to the strong localized perturbation analysis initially developed in Refs. 43 and 44.

Naturally, our model and mathematical analysis made a number of simplifying assumptions. For one, we assumed that the patchy surface (left wall) is flat (see Fig. 1). If the surface had some curvature, then our results require that the characteristic radius of curvature of the surface, Rc > 0, is much greater than the radius of the particle. That is, we require
Rc/R1
(51)
so that the surface is flat at the length scale of the particle. In addition, as in the analysis of a patchy particle and a uniformly reactive surface in Ref. 25, applying our results to a curved patchy surface requires that the particle “forgets” its orientation over timescales in which the wall curvature affects the particle. Since (i) the particle orientation relaxes to the uniform distribution at rate Drot and (ii) the surface curvature becomes relevant on the timescale Rc2/Dtr, we thus require
DrotRc2/Dtr=β2Rc2/R21.
(52)
Requirements (51) and (52) are equivalent for the typical case that the particle’s rotational diffusivity and translational diffusivity are comparable (i.e., β ≫̸ 1 and β ≪̸ 1).

We also assumed that the absorbing fractions were both small (i.e., f1 ≪ 1, f2 ≪ 1). Although simulations showed that the theory can be fairly accurate away from this parameter regime (see Sec. VI C), an interesting avenue for future work is to devise more accurate approximations when one or both of the absorbing fractions are not small. One possible approach is to interpolate between our results for the regime f1 ≪ 1, f2 ≪ 1 and the prior results for the regimes f1 ≈ 1, f2 ≪ 1 and f1 ≪ 1, f2 ≈ 1.

The authors acknowledge the support from the National Science Foundation (Grant Nos. NSF CAREER DMS-1944574 and DMS-1814832).

The authors have no conflicts to disclose.

Claire E. Plunkett: Conceptualization (equal); Formal analysis (equal); Software (lead); Writing – original draft (equal); Writing – review & editing (equal). Sean David Lawley: Conceptualization (equal); Formal analysis (equal); Funding acquisition (lead); Software (supporting); Writing – original draft (equal); Writing – review & editing (equal).

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

We now derive the trapping rate κ1 in (3) for a patchy surface and a perfectly reactive particle. The formula in (3) can be obtained via several different approaches,1,11,13,14 but we give here a simplified version of the approach used in Sec. III above.

The scenario of a patchy surface and a perfectly reactive particle fits into the framework of Sec. II by assuming that the particle is completely covered in patches, i.e., f2 = 1. The survival probability in (10) then satisfies (11)(15) with Γ2 = [0, π] × [0, 2π). That is, the angular orientation of the particle (θ, φ) is irrelevant, and the diffusion equation in (11) is only in time t and Cartesian coordinates (x, y, z).

Similar to Sec. III, we fix any time t > 0 and introduce an outer expansion, which is valid away from the patches,
S1+εS1+.
(A1)
Note that the ɛ term in (A1) replaces the ɛ3 term for the full patchy sphere and patchy particle problem in (20). Plugging (A1) into (11) and using that Γ2 = [0, π] × [0, 2π) show that S1 satisfies
tS1=LS1,(x,y,z)[0,)×R2,xS1=0,x=0,(y,z)m=1(ym,zm).
(A2)
To determine the behavior of S1 near the mth trap on the left wall (0, ym, zm), we introduce stretched coordinates for fixed m ∈ {1, 2, …},
μ=ε1x/R,ν=ε1(yym)/R,η=ε1(zzm)/R.
We then define the inner solution w as
w(μ,ν,η,t)S(εRμ,ym+εRν,zm+εRη,t).
(A3)
Using (11), the inner solution w satisfies
tw=ε2R2Dtr(μμ+νν+ηη)w+O(ε1).
(A4)
Therefore, plugging the inner expansion w = w0 + O(ɛ) into (A4) implies that w0 is harmonic,
(μμ+νν+ηη)w0=0,μ>0.
(A5)
Using (15), w0 satisfies the following boundary conditions:
w0=0,μ=0,ν2+η2<a12,μw0=0,μ=0,otherwise.
(A6)
The PDE in (A5) and (A6) is the so-called electrified disk problem, and the solution is known analytically.45 For our purposes, we only need the far-field behavior of w0, which is
w0α1a1(2/π)ρ as ρμ2+ν2+η2,
(A7)
where α is determined by matching with the outer solution. The matching condition is that the near-field behavior of the outer expansion as (x, y, z) → (0, ym, zm) must agree with the far-field behavior of the inner expansion as ρ → ∞. That is,
1+εS1w0as (x,y,z)(0,ym,zm),ρ.
Using (A3) and (A7), it follows that α = 1 and
S1a1(2/π)xR2+yymR2+zzmR2as (x,y,z)(0,ym,zm).
(A8)
Writing the singular behavior (A8) in the distributional form for each m ∈ {1, 2, …}, the boundary condition in (A2) becomes
xS1=4a1Rm=1δ(yym)δ(zzm),x=0.
(A9)
Define
S̄(x,t)liml1(2l)2[l,l]2S(x,y,z,t)dydz
and
κ̄DtrxS̄(0,t)S̄(0,t)
(A10)
so that it is a tautology that S̄ satisfies DtrxS̄=κ̄S̄ at x = 0. We now derive the effective trapping rate by determining the behavior of (A10) as ɛ → 0.
The denominator in (A10) approaches unity as ɛ → 0 due to (19). To determine the behavior of the numerator in (A10), we interchange the derivative with the limit and integral, recall the outer expansion in (A1), and apply (A9),
κ̄Dtrliml1(2l)2[l,l]2xS(0,y,z,t)dydzεN1DtrRa1πas ε0.
(A11)
Using (2) and (18) yields (3).
To derive the distributional form in (29) of the singular behavior in (28), suppose that a function f satisfies
tf=Lf,(x,y,z,θ,φ)Ω
and the following boundary condition at x = 0:
xf=Km=1n=1Nδ(yym)δ(zzm)δ(θθn)sin(θn)δ(φφn).
To obtain the singular behavior of f as (x, y, z, θ, φ) → (0, ym, zm, θn, φn), define the inner solution analogous to (22),
g(μ,ν,η,s,p,t)fεRμ,ym+εRν,zm+εRη,θn+εβs,φn+εβsinθnp,t,
and introduce the inner expansion gɛ−3g0 +⋯. By the same argument that yielded (24) and (25), we have that g0 is harmonic,
(μμ+νν+ηη+ss+pp)g0=0,μ>0,
(B1)
and g0 satisfies the following boundary condition at μ = 0:
μg0=ε4R4π2cRβ2δ(εRν)δ(εRη)δ(εβs)sin(θn)δεβpsin(θn)=4π2cδ(ν)δ(η)δ(s)δ(p),
(B2)
where we have used the identity δ(αx) = δ(x)/|α|. Ignoring arbitrary additive constants, the solution to (B1) and (B2) is26 
g0=c(μ2+ν2+η2+s2+p2)3/2.
Matching the far-field behavior of g with the near-field behavior of f shows that f indeed has the singular behavior in (28).
We now detail the method used in Sec. IV to calculate the electrostatic capacitance c0. The method relies on a probabilistic representation of the solution w0 of (24) and (25). Let Z(t)R5 be a standard five-dimensional Brownian motion. Define the first time that this process reaches the duocylinder D in (27),
τinf{t>0:Z(t)D}.
The leading order inner solution w0 satisfying (24) and (25) can be written as
w0(μ,ν,η,s,p)=1q(μ,ν,η,s,p),
where q is the probability that Z eventually reaches D, conditioned on the initial position of Z,
q(μ,ν,η,s,p)=P(τ<|Z(0)=(μ,ν,η,s,p)).
The function q must be harmonic for μ ≠ 0 and satisfy the boundary conditions at μ = 0,
q=1,μ=0,(ν,η,s,p)D,μq=0,μ=0,(ν,η,s,p)D.
Let q̄(ρ) denote the average of q over the surface of the five-dimensional ball of radius ρ=μ2+ν2+η2+s2+p2>0 centered at the origin. If μ = 0 and ρ > 0 is such that
ρ>ρ0(b1)2+(b2)2,
then μq = 0. Hence, integrating the harmonic PDE for q over the surface of the 5D ball of radius ρ > ρ0 yields the following ordinary differential equation for q̄(ρ):
(4/ρ)ρ+ρρq̄=0,ρ>ρ0.
Solving this equation and using the far-field behavior of w0 in (26) yield
q̄(ρ)=c0/ρ3,ρ>ρ0.
(C1)

Equation (C1) implies that we can calculate c0 by calculating the probability q̄(ρ) that the five-dimensional Brownian motion Z eventually reaches the duocylinder D defined in (27), conditioned that Z is initially uniformly distributed on a ball of radius ρ > ρ0. We can estimate q̄(ρ) by simulating M ≫ 1 realizations of Z and calculating the proportion of these M trials, which reach D before reaching some large outer radius ρρ0.

An efficient way to simulate these instances is via a kinetic Monte Carlo method. The kinetic Monte Carlo method breaks the simulation process into two stages of simpler diffusion processes that can be exactly simulated and then alternates between these two stages until reaching a break point. The only error in this method results from the finite outer radius ρ < ∞ and the finite number of trials M < ∞, both of which can be taken very large due to the computational efficiency of the method.

Simulations are initialized by placing the “particle” Z on the five-dimensional sphere of radius ρ centered at the origin according to a uniform distribution. The method then utilizes the following two stages developed by Bernoff et al.12 for a similar problem in R3.

  • Stage I: Project from bulk to plane: The particle is projected to the μ = 0 plane following the exact distribution given below. If the particle lands within D, then this event is recorded and the trial ends. If not, the algorithm proceeds to stage II.

  • Stage II: Project from plane to bulk: A distance d > 0 is calculated, which is less than or equal to the distance from the current particle location to D. The particle is then projected to a uniformly distributed point on the 5D sphere of radius d. If the particle reaches a distance larger than ρ, then this event is recorded and the trial ends. Otherwise, the algorithms returns to stage I.

The distribution in stage I is calculated by first sampling the random time that it takes the particle to reach μ = 0, which is12 
t*=14zerfc1(U)2,
where U is uniformly distributed on [0, 1]. Then, given Z at position (μ,ν,η,s,p)R5 at the beginning of stage I, the position at the end of stage I is
(0,ν,η,s,p)+2t*(0,ξ1,ξ2,ξ3,ξ4)R5,
where ξ1, ξ2, ξ3, ξ4 are independent standard normal random variables.
The goal of stage II is to propagate the particle as far as possible while ensuring that the particle cannot reach D during this propagation. Let (0, ν, η, s, p) be the position of the particle at the beginning of stage II. Then, given that the algorithm is in stage II, we have that ZD and, therefore,
d1=ν2+η2b1>0 and/or d2=s2+p2b2>0,
where b1 and b2/β are the two radii of the duocylinder in (27). If we define the distance d ≔ max{d1, d2}, then the five-dimensional sphere of radius d centered at (0,ν,η,s,p)D cannot intersect D. Therefore, stage II places the particle uniformly on the surface of this 5D sphere.

We now describe the numerical implementation of the kinetic Monte Carlo algorithm of Sec. VI used to simulate the full stochastic system.

The left wall has reactive patches of radius r1 > 0 with density N1/(4π) placed on a square lattice such that the centers are
(ȳi,z̄j)=4πiN1,4πjN1,i,jZ.
Simulations are initialized by placing the unit spherical particle such that the nearest point of the particle to the left wall is (x0, y0, z0) at t = 0, where we fix starting distance from the left wall x0 and choose (y0, z0) from a uniform distribution on [2πN1,2πN1]2. The sphere has N2 (with N2 odd) reactive spherical caps of polar angle r2 > 0 with centers (θ̄k,φ̄k) on a Fibonacci lattice,20 
θ̄k=arccos2kN21,φ̄k=4πk1+5,k=1,,N2.
The particle is given initial orientation (θ0, ϕ0) sampled from a uniform distribution on the unit sphere using
θ0=arccos(2ξ11),φ0=2πξ2,
where (ξ1, ξ2) is uniformly distributed on [0, 1]2.

1. Stage I: Project from the bulk to the plane

This stage follows the algorithm of Bernoff et al.12 for the translational diffusion of the sphere plus rotational diffusion. Given a particle located at (x0, y0, z0) at time t = 0, we sample the transit time to the particle’s first contact with the left wall by sampling a uniform random number ξ ∈ [0, 1] and calculating
t1=14Dtrx0erfc1(ξ)2.
This time is added to the cumulative elapsed time of the simulation. The spatial location of the particle at the arrival time t1 is (0, y1, z1), where y1 and z1 are independently drawn from normal distributions with respective means y0 and z0 and variances 2Dtrt1.
We also sample (θ1, φ1), the orientation of the sphere following the propagation of the sphere to the wall, given time t1. The authors of Ref. 46 gave series solutions f(θ, t) to the equation for rotational diffusion of a point in angular coordinates (θ, φ) on a sphere of radius 1, given that the point is initially at the pole of the sphere,
tf=DrotΔf,f(θ,0)=δ(θ).
We numerically integrate the solution f(θ, t) from Ref. 46, which yields the following cumulative distribution function:
F(θ,t)=2π0θf(θ,t)sinθdθ.
(D1)
From this, we can sample the change θd in θ, given t1. Due to rotational symmetry, we sample the change φd in φ from a uniform distribution on [0, 2π). We then update the orientation of the sphere to (θ1, φ1), given initial orientation (θ0, φ0) by rotating (θd, φd), implemented using rotation operators by interpreting (θ0, φ0) as a point on the Bloch sphere.47 
At this point, the algorithm checks to see if a patch on the sphere is touching a patch on the plane. That is, the algorithm checks if the position coordinates (y1, z1) are within a patch on the plane and if the orientation coordinates (θ1, φ1) are within a patch on the sphere. To do so, the signed distances dA and dB to the nearest patches on both the plane and the sphere, respectively, are calculated. The signed distance dA to the nearest patch on the plane is given by
dA=y1mod4πN12+z1mod4πN12r1.
The signed distance dB is the signed arc length between the point of contact of the sphere and the plane to the nearest patch on the sphere. This is calculated by taking the arc length48  ai to the center of the ith patch and subtracting the patch radius r2 from the minimum arc length min ai,
ai=tan1(sinθ̄isin(Δφi))2+(sinθ1cosθ̄icosθ1sinθ̄jcos(Δφi))2cosθ1cosθ̄i+sinθ1sinθ̄icos(Δφi) for Δφiφ̄iφ1,dB=minair2.
If dA < 0 and dB < 0, then the trial is complete and the final elapsed time is recorded. If dA > 0, the algorithm proceeds to case A of stage II. If dA < 0 and dB > 0, the algorithm proceeds to case B of stage II.

2. Stage II, case A: Reflect from the plane

Given that the particle is located at (0, y1, z1) at some time t > 0 at the start of stage II, we propagate the particle to the hemisphere of radius dA centered at (0, y1, z1), which ensures that the end condition could not have been met. We sample a time t2 from the well-known distribution for a three-dimensional diffusion process to exit a sphere of radius dA (see, for example, Ref. 12). The time t2 is added to the cumulative elapsed time.

The particle is placed at (x2, y2, z2), sampled uniformly on the hemisphere of radius dA centered at (0, y1, z1) by symmetry of Brownian motion. The orientation of the particle is updated as it was in stage I, where we sample from the distribution F(θ, t2) in (D1), given t2 for the change θd in θ, and sample the change φd in φ from a uniform distribution on [0, 2π) and update (θ1, φ1) with rotational change (θd, φd) to obtain new orientation (θ2, φ2).

At this point, the algorithm returns to stage I, updating the position and orientation labels (x2, y2, z2) and (θ2, φ2) to (x0, y0, z0) and (θ0, φ0).

3. Stage II, case B: Reflect only from the sphere

In the case that the particle touches the plane such that it is touching a patch on the plane at the end of stage I, we sample the time it takes (θ1, φ1) to rotate a distance of dB (the distance to the nearest patch on the particle).

To do so, we consider the survival probability SdB(θ,t) for rotating a distance of dB by time t, given that the initial spherical coordinates are (0, φ1) at t = 0. By rotational symmetry, we can ignore φ1. This survival probability satisfies the spherical diffusion equation,
tSdB=(tanθ)1θSdB+2θ2SdB,t>0,θ(0,dB),SdB=1,t=0,θSdB=0,θ=0,SdB=0,θ=dB.
We numerically solve this over a mesh of values of dB, 0 < θ < dB, and t using the MATLAB pdepe function.49 To sample an arrival time to the spherical cap of radius dB, we sample a uniform random number ξ ∈ [0, 1] and numerically solve
SdB(0,s)=ξ
for s and then compute the exit time t2 = s/Drot. The time t2 is added to the elapsed time, and the angular position is updated from (θ1, φ1) with rotational change (dB, φd), where φd is sampled from a uniform distribution on [0, 2π), to obtain new orientation (θ2, φ2), as described in stage I. The updated location (x2, y2, z2) of the particle after time t2 is sampled by generating three standard normal random numbers (ξ1, ξ2, ξ3) and computing
x2=2Dtrt2ξ1,y2=y1+2Dtrt2ξ2,z2=z1+2Dtrt2ξ3.
The algorithm then returns to stage I, updating the position and orientation labels (x2, y2, z2) and (θ2, φ2) to (x0, y0, z0) and (θ0, φ0).
1.
H. C.
Berg
and
E. M.
Purcell
, “
Physics of chemoreception
,”
Biophys. J.
20
(
2
),
193
219
(
1977
).
2.
A.
Saxena
,
B. P.
Tripathi
,
M.
Kumar
, and
V. K.
Shahi
, “
Membrane-based techniques for the separation and purification of proteins: An overview
,”
Adv. Colloid Interface Sci.
145
(
1-2
),
1
22
(
2009
).
3.
R. W.
Baker
and
B. T.
Low
, “
Gas separation membrane materials: A perspective
,”
Macromolecules
47
(
20
),
6999
7013
(
2014
).
4.
F.
Keil
, “
Diffusion and reaction in porous networks
,”
Catal. Today
53
(
2
),
245
258
(
1999
).
5.
B. R.
Scharifker
, “
Diffusion to ensembles of microelectrodes
,”
J. Electroanal. Chem. Interfacial Electrochem.
240
(
1-2
),
61
76
(
1988
).
6.
H. T.
Brown
and
F.
Escombe
, “
Static diffusion of gases and liquids in relation to the assimilation of carbon and translocation in plants
,”
Philos. Trans. R. Soc., B
193
(
185-193
),
223
291
(
1900
).
7.
A.
Wolf
,
W. R.
Anderegg
, and
S. W.
Pacala
, “
Optimal stomatal behavior with competition for water and risk of hydraulic impairment
,”
Proc. Natl. Acad. Sci. U. S. A.
113
(
46
),
E7222
E7230
(
2016
).
8.
F.
Piazza
and
S. D.
Traytak
, “
Diffusion-influenced reactions in a hollow nano-reactor with a circular hole
,”
Phys. Chem. Chem. Phys.
17
(
16
),
10417
10425
(
2015
).
9.
F.
Piazza
, “
The physics of boundary conditions in reaction-diffusion problems
,”
J. Chem. Phys.
157
,
234110
(
2022
).
10.
G.
Handy
and
S. D.
Lawley
, “
Revising Berg-Purcell for finite receptor kinetics
,”
Biophys. J.
120
(
11
),
2237
2248
(
2021
).
11.
A. M.
Berezhkovskii
,
Y. A.
Makhnovskii
,
M. I.
Monine
,
V. Y.
Zitserman
, and
S. Y.
Shvartsman
, “
Boundary homogenization for trapping by patchy surfaces
,”
J. Chem. Phys.
121
(
22
),
11390
11394
(
2004
).
12.
A. J.
Bernoff
,
A. E.
Lindsay
, and
D. D.
Schmidt
, “
Boundary homogenization and capture time distributions of semipermeable membranes with periodic patterns of reactive sites
,”
Multiscale Model. Simul.
16
(
3
),
1411
1447
(
2018
).
13.
D.
Shoup
and
A.
Szabo
, “
Role of diffusion in ligand binding to macromolecules and cell-bound receptors
,”
Biophys. J.
40
(
1
),
33
39
(
1982
).
14.
R.
Zwanzig
, “
Diffusion-controlled ligand binding to spheres partially covered by receptors: An effective medium treatment
,”
Proc. Natl. Acad. Sci. U. S. A.
87
(
15
),
5856
5857
(
1990
).
15.
A. M.
Berezhkovskii
,
M. I.
Monine
,
C. B.
Muratov
, and
S. Y.
Shvartsman
, “
Homogenization of boundary conditions for surfaces with regular arrays of traps
,”
J. Chem. Phys.
124
(
3
),
036103
(
2006
).
16.
C. B.
Muratov
and
S. Y.
Shvartsman
, “
Boundary homogenization for periodic arrays of absorbers
,”
Multiscale Model. Simul.
7
(
1
),
44
61
(
2008
).
17.
A. F.
Cheviakov
,
A. S.
Reimer
, and
M. J.
Ward
, “
Mathematical modeling and numerical computation of narrow escape problems
,”
Phys. Rev. E
85
(
2
),
021131
(
2012
).
18.
L.
Dagdug
,
M. V.
Vázquez
,
A. M.
Berezhkovskii
, and
V. Y.
Zitserman
, “
Boundary homogenization for a sphere with an absorbing cap of arbitrary size
,”
J. Chem. Phys.
145
(
21
),
214101
(
2016
).
19.
A. E.
Lindsay
,
A. J.
Bernoff
, and
M. J.
Ward
, “
First passage statistics for the capture of a Brownian particle by a structured spherical target with multiple surface traps
,”
Multiscale Model. Simul.
15
(
1
),
74
109
(
2017
).
20.
A. J.
Bernoff
and
A. E.
Lindsay
, “
Numerical approximation of diffusive capture rates by planar and spherical surfaces with absorbing pores
,”
SIAM J. Appl. Math.
78
(
1
),
266
290
(
2018
).
21.
S. D.
Lawley
and
C. E.
Miles
, “
How receptor surface diffusion and cell rotation increase association rates
,”
SIAM J. Appl. Math.
79
(
3
),
1124
1146
(
2019
).
22.
K.
Šolc
and
W. H.
Stockmayer
, “
Kinetics of diffusion-controlled reaction between chemically asymmetric molecules. I. General theory
,”
J. Chem. Phys.
54
(
7
),
2981
2988
(
1971
).
23.
K.
Šolc
and
W.
Stockmayer
, “
Kinetics of diffusion-controlled reaction between chemically asymmetric molecules. II. Approximate steady-state solution
,”
Int. J. Chem. Kinet.
5
(
5
),
733
752
(
1973
).
24.
H. C. R.
Klein
and
U. S.
Schwarz
, “
Studying protein assembly with reversible brownian dynamics of patchy particles
,”
J. Chem. Phys.
140
(
18
),
184112
(
2014
).
25.
S. D.
Lawley
, “
Boundary homogenization for trapping patchy particles
,”
Phys. Rev. E
100
(
3
),
032601
(
2019
).
26.
C. E.
Plunkett
and
S. D.
Lawley
, “
Bimolecular binding rates for pairs of spherical molecules with small binding sites
,”
SIAM J. Multiscale Model. Simul.
19
,
148
(
2021
).
27.
C.
McMullen
,
The Visual Guide to Extra Dimensions: Visualizing the Fourth Dimension, Higher-Dimensional Polytopes, and Curved Hypersurfaces
(
Custom Books
,
Lexington
,
2008
).
28.
D. S.
Grebenkov
, “
Paradigm shift in diffusion-mediated surface phenomena
,”
Phys. Rev. Lett.
125
(
7
),
078102
(
2020
).
29.
G. A.
Pavliotis
,
Stochastic Processes and Applications
(
Springer
,
2016
).
30.
J. D.
Jackson
,
Classical Electrodynamics
, 2nd ed.(
Wiley
,
New York
,
1975
).
31.
M. E.
Muller
, “
Some continuous Monte Carlo methods for the Dirichlet problem
,”
Ann. Math. Stat.
27
(
3
),
569
589
(
1956
).
32.
A. G.
Belyaev
,
G. A.
Chechkin
, and
R. R.
Gadyl’shin
, “
Effective membrane permeability: Estimates and low concentration asymptotics
,”
SIAM J. Appl. Math.
60
(
1
),
84
108
(
1999
).
33.
M.
Bruna
,
S. J.
Chapman
, and
G. Z.
Ramon
, “
The effective flux through a thin-film composite membrane
,”
Europhys. Lett.
110
(
4
),
40005
(
2015
).
34.
A. C.
Newton
,
J.
Groenewold
,
W. K.
Kegel
, and
P. G.
Bolhuis
, “
Rotational diffusion affects the dynamical self-assembly pathways of patchy particles
,”
Proc. Natl. Acad. Sci. U. S. A.
112
(
50
),
15308
15313
(
2014
).
35.
C. J.
Roberts
and
M. A.
Blanco
, “
Role of anisotropic interactions for proteins and patchy nanoparticles
,”
J. Phys. Chem. B
118
(
44
),
12599
12611
(
2014
).
36.
P. C.
Bressloff
and
S. D.
Lawley
, “
Escape from subcellular domains with randomly switching boundaries
,”
Multiscale Model. Simul.
13
(
4
),
1420
1445
(
2015
).
37.
P. C.
Bressloff
and
S. D.
Lawley
, “
Stochastically gated diffusion-limited reactions for a small target in a bounded domain
,”
Phys. Rev. E
92
(
6
),
062117
(
2015
).
38.
A. F.
Cheviakov
and
M. J.
Ward
, “
Optimizing the principal eigenvalue of the Laplacian in a sphere with interior traps
,”
Math. Comput. Modell.
53
(
7-8
),
1394
1409
(
2011
).
39.
A. F.
Cheviakov
,
M. J.
Ward
, and
R.
Straube
, “
An asymptotic analysis of the mean first passage time for narrow escape problems: Part II: The sphere
,”
Multiscale Model. Simul.
8
(
3
),
836
870
(
2010
).
40.
D.
Coombs
,
R.
Straube
, and
M.
Ward
, “
Diffusion on a sphere with localized traps: Mean first passage time, eigenvalue asymptotics, and Fekete points
,”
SIAM J. Appl. Math.
70
(
1
),
302
332
(
2009
).
41.
S. D.
Lawley
and
C. E.
Miles
, “
Diffusive search for diffusing targets with fluctuating diffusivity and gating
,”
J. Nonlinear Sci.
29
(
6
),
2955
2985
(
2019
).
42.
S. D.
Lawley
and
V.
Shankar
, “
Asymptotic and numerical analysis of a stochastic PDE model of volume transmission
,”
Multiscale Model. Simul.
18
(
2
),
887
915
(
2020
).
43.
M. J.
Ward
,
W. D.
Heshaw
, and
J. B.
Keller
, “
Summing logarithmic expansions for singularly perturbed eigenvalue problems
,”
SIAM J. Appl. Math.
53
(
3
),
799
828
(
1993
).
44.
M. J.
Ward
and
J. B.
Keller
, “
Strong localized perturbations of eigenvalue problems
,”
SIAM J. Appl. Math.
53
(
3
),
770
798
(
1993
).
45.
I. N.
Sneddon
,
Mixed Boundary Value Problems in Potential Theory
(
North-Holland Publishing Company
,
1966
).
46.
V.
Tulovsky
and
L.
Papiez
, “
Formula for the fundamental solution of the heat equation on the sphere
,”
Appl. Math. Lett.
14
(
7
),
881
884
(
2001
).
47.
M. A.
Nielsen
and
I. L.
Chuang
,
Quantum Computation and Quantum Information
, 10th anniversary ed. (
Cambridge University Press
,
Cambridge, UK
,
2010
).
48.
T.
Vincenty
, “
Direct and inverse solutions of geodesics on the ellipsoid with application of nested equations
,”
Surv. Rev.
23
(
176
),
88
93
(
1975
).
49.
The Mathworks, Inc.
, Natick, Massachusetts. MATLAB version 9.10.0.1739362 (R2021a) Update 5, 2021.