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.

## I. INTRODUCTION

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 filtration^{2} 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 scrutiny^{8–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.

*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)\u2208R2$ is the location in the plane parallel to the surface. The concentration evolves according to the diffusion equation,

*D*

_{tr}> 0 denotes the translational diffusivity of the particle. Trapping at the patchy surface is described by mixed boundary conditions,

*κ*

_{1}> 0:

*r*

_{1}> 0 and occupy a small fraction

*N*

_{1}> 0 is the average number of patches in an area of size 4

*πR*

^{2}for some length scale

*R*> 0 [i.e.,

*N*

_{1}/(4

*πR*

^{2}) 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}

^{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,

*D*

_{tr}> 0 and

*D*

_{rot}> 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:

*R*> 0 containing

*N*

_{2}≫ 1 roughly evenly distributed patches of common radius

*r*

_{2}≪

*R*occupying a small fraction

^{25}

*βr*

_{1}/

*r*

_{2}. In particular,

*χ*depends on the electrostatic capacitance of a certain four-dimensional object called a duocylinder

^{27}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.

^{28}to derive the following heuristic trapping rate:

*κ*

_{0}is a remarkably accurate approximation to the asymptotic trapping rate

*κ*in (7). Indeed, the accuracy depends on the value of

*βr*

_{1}/

*r*

_{2}, but the relative error is never more than 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.

## II. MATHEMATICAL MODEL

### A. Stochastic formulation

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

*t*≥ 0, let

*X*(

*t*) ∈ [0,

*L*] denote the distance between the particle and the left wall and let $(Y(t),Z(t))\u2208R2$ 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,

*D*

_{tr}> 0, which means that the coordinates (

*X*(

*t*),

*Y*(

*t*),

*Z*(

*t*)) satisfy the following stochastic differential equations (SDEs):

*W*

_{X},

*W*

_{Y}, and

*W*

_{Z}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

*D*

_{rot}> 0, which means its orientation (Θ(

*t*), Φ(

*t*)) satisfies the following SDEs:

*W*

_{Θ}and

*W*

_{Φ}are independent, standard Brownian motions.

*τ*. To describe the patches on the left wall, let

*γ*

_{1}(

*y*

_{0},

*z*

_{0},

*r*) denote the disk of radius

*r*> 0 centered at $(y0,z0)\u2208R2$,

*γ*

_{2}(

*θ*

_{0},

*φ*

_{0},

*a*) denote the spherical cap with polar angle

*a*> 0 centered at (

*θ*

_{0},

*φ*

_{0}),

*N*

_{2}≥ 1 patches on the particle, then the set of all patches on the particle is

*τ*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}).

### B. Survival probability

*τ*in (9) is described by the survival probability conditioned on the initial state of the system,

^{29}

*x*,

*y*,

*z*),

*θ*,

*φ*),

*x*=

*L*if

*L*< ∞,

*L*= ∞,

*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,

## III. MATCHED ASYMPTOTIC ANALYSIS

*κ*, assuming that the patches (i) are roughly evenly distributed and (ii) occupy a small surface area fraction of the wall (

*f*

_{1}≪ 1) and the particle (

*f*

_{2}≪ 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\u0304(x,t)$ to the following one-dimensional diffusion equation:

*x*= 0:

*r*

_{1}> 0 and the patches on the particle are spherical caps with polar angle

*r*

_{2}/

*R*> 0. Assume that

*ɛ*≪ 1 and

*a*

_{1}> 0,

*a*

_{2}> 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

*a*

_{1}= 1 or

*a*

_{2}= 1, but we keep both parameters for pedagogical reasons). To derive the trapping rate

*κ*, we apply matched asymptotic analysis

^{19,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).

### A. Outer expansion

*t*> 0. As

*ɛ*→ 0, the surfaces of the particle and the left wall become perfectly reflecting, and thus,

*S*has a boundary layer in a neighborhood of each possible intersection of patches, so we introduce the outer expansion

*ɛ*

^{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

*f*

_{1}=

*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

*ɛ*].

*S*

_{1}must satisfy

### B. Inner expansion

*S*

_{1}as the

*n*th trap on the particle approaches the

*m*th trap on the left wall,

*m*∈ {1, 2, …} and

*n*∈ {1, …,

*N*

_{2}}:

*w*as a function of the stretched coordinates,

*w*satisfies

*β*is defined so that

*D*

_{rot}

*β*

^{−2}=

*R*

^{−2}

*D*

_{tr}. Therefore, plugging the inner expansion

*w*

_{0}is harmonic,

*w*

_{0}satisfies the following mixed boundary conditions at

*μ*= 0:

*b*

_{1}=

*a*

_{1}and

*b*

_{2}=

*a*

_{2}/

*β*.

### C. Matching

^{30}that

*w*

_{0}has the far-field behavior,

*α*is a constant determined by matching with the outer solution and

*a*

_{1}and

*a*

_{2}/

*β*. In particular,

*c*

_{0}(

*b*

_{1},

*b*

_{2}) is the electrostatic capacitance of a so-called duocylinder

^{27}embedded in $R5$,

*c*

_{0}(

*a*

_{1},

*a*

_{2}/

*β*) until Sec. IV and proceed in terms of

*c*=

*c*

_{0}(

*a*

_{1},

*a*

_{2}/

*β*).

*x*,

*y*,

*z*,

*θ*,

*φ*) → (0,

*y*

_{m},

*z*

_{m},

*θ*

_{n},

*φ*

_{n}) must agree with the far-field behavior of the inner expansion as

*ρ*→ ∞. That is,

*α*= 1 and that

*S*

_{1}has the following singular behavior as (

*x*,

*y*,

*z*,

*θ*,

*φ*) → (0,

*y*

_{m},

*z*

_{m},

*θ*

_{n},

*φ*

_{n}):

*n*∈ {1, …,

*N*

_{2}} and

*m*∈ {1, 2, …}, the boundary condition in (21) at

*x*= 0 becomes

*K*= 4

*π*

^{2}

*cRβ*

^{2}. The derivation of the distributional form in (29) of the singular behavior in (28) is given in Appendix B.

### D. Effective trapping rate *κ*

*N*

_{1}/(4

*πR*

^{2}), which means that

*N*

_{1}∈ (0, ∞). The integral in (30) is simply the number of patches in the square [−

*l*,

*l*]

^{2}, and thus,

*N*

_{1}is the average number of patches on the wall in an area equal to the surface area of the particle (note that

*N*

_{1}need not be an integer).

*ɛ*→ 0.

*ɛ*→ 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),

*c*=

*c*

_{0}(

*a*

_{1},

*a*

_{2}/

*β*), we may rewrite (A11) as

*κ*

_{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,

*χ*is the following function of

*r*

_{1}

*β*/

*r*

_{2}:

## IV. KINETIC MONTE CARLO CALCULATION OF CAPACITANCE *c*_{0}

The homogenized trapping rate *κ* in (37) involves the factor *χ* in (38), which depends on the electrostatic capacitance *c*_{0} 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}

### A. Summary of computational method

We now summarize this calculation of *c*_{0}. The details are given in Appendix C. The leading order inner solution *w*_{0} 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 *c*_{0} is determined by the far-field behavior of *w*_{0} [see (26)], we calculate *c*_{0} 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 *c*_{0} 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 *c*_{0} is very small.

*c*

_{0}. It is immediate that the capacitance is symmetric,

*b*∈ (0, 1]. The black solid curve in Fig. 2 is from the kinetic Monte Carlo algorithm described above with outer radius

*ρ*

_{∞}= 10

^{5}and

*M*= 10

^{7}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 $\rho \u221e\u22123$ as

*ρ*

_{∞}→ ∞. Hence, setting

*ρ*

_{∞}= 10

^{5}means that relative error in calculating the capacitance from the finite outer radius is on the order of 10

^{−15}.

## V. HEURISTICS

### A. Heuristic combined trapping rate

*κ*

_{1}in (3), then the survival probability

*S*

_{1}=

*S*

_{1}(

*x*,

*θ*,

*φ*,

*t*) satisfies

*S*

_{1}has the following probabilistic representation:

^{28}

*ℓ*

_{X,Θ,Φ}(

*t*) is the following local time:

_{A}denotes the indicator function on an event

*A*(i.e., 1

_{A}= 1 if

*A*occurs and 1

_{A}= 0 otherwise) and

*E*

_{1}is an independent exponential random variable with rate

*κ*

_{1}.

*κ*

_{2}in (6), then the survival probability

*S*

_{2}=

*S*

_{2}(

*x*,

*y*,

*z*,

*t*) satisfies

*S*

_{2}has the following probabilistic representation:

*ℓ*

_{X,Y,Z}(

*t*) is the following local time:

*E*

_{2}is an independent exponential random variable with rate

*κ*

_{2}.

*τ*

_{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

*ℓ*(

*t*) is the local time of

*X*at

*x*= 0,

*Y*(0),

*Z*(0)) yields

*τ*

_{0}in (9) except that we replace

*ℓ*

_{X,Θ,Φ}and

*ℓ*

_{X,Y,Z}by their averaged versions in (41) and (42),

*E*

_{1}/

*f*

_{2}is exponentially distributed with rate

*f*

_{2}

*κ*

_{1}(since

*E*

_{1}is exponentially distributed with rate

*κ*

_{1}). Similarly,

*E*

_{2}/

*f*

_{1}is exponentially distributed with rate

*f*

_{1}

*κ*

_{2}(since

*E*

_{2}is exponentially distributed with rate

*κ*

_{2}). Furthermore, the minimum of two exponential random variables with rates

*f*

_{1}

*κ*

_{2}and

*f*

_{2}

*κ*

_{1}is exponentially distributed with rate

*f*

_{1}

*κ*

_{2}+

*f*

_{2}

*κ*

_{1}. Therefore, $\tau 0\u0304$ is equal in distribution to

*E*

_{0}is exponentially distributed with rate

*κ*

_{0}in (43) at

*x*= 0,

### B. Comparing heuristic *κ*_{0} and asymptotic *κ*

*κ*

_{0}in (43) with the asymptotic trapping rate

*κ*in (37). Using (35) and (36), we can write

*κ*

_{0}in the following form:

*κ*

_{0}≈

*κ*if the following formula approximates the capacitance of the duocylinder:

*c*

_{0}(1,

*b*) computed from kinetic Monte Carlo simulations (see Sec. IV) and the red dashed curve is $c\u0304(1,b)=b(1+b)/\pi $. Clearly, $c0\u2260c\u0304$, but the approximation is reasonably accurate. Indeed, the relative error $|c0\u2212c\u0304|/c0$ is at most 16%, and thus, the relative error between

*κ*

_{0}and

*κ*is at most 16%,

*ζ*so that the resulting approximation agrees with

*c*

_{0}(1,

*b*) at

*b*= 1. Specifically, the black circles in Fig. 2 are

*b*∈ (0.1, 1].

## VI. KINETIC MONTE CARLO SIMULATIONS OF THE FULL STOCHASTIC SYSTEM

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*d*_{A}to the set of patches on the plane is calculated. The sphere is then propagated into the bulk to a hemisphere of radius*d*_{A}. 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*d*_{B}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*d*_{B}, 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.

### A. Numerically optimized the trapping rate

^{25}

*κ*

_{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

*t*

_{1}<

*t*

_{2}<⋯<

*t*

_{M}and define the empirical survival probability,

*t*

_{j}.

^{12}We then choose

*κ*

_{opt}so that the solution $S\u0304$ to the homogenized, one-dimensional problem in (16) and (17) approximates the empirical survival probability

*S*

_{emp}.

*κ*, the solution $S\u0304$ to (12), (14), (16), and (17) is

*u*) denotes the complementary error function and erfcx(

*u*) = exp(

*u*

^{2})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\u0304$ and

*S*

_{emp},

*κ*

_{opt}is then the value of the trapping rate, which minimizes $E$.

### B. Results for *ɛ* → 0

*M*= 10

^{6}stochastic realizations of the reaction times using the stochastic simulation algorithm described above for each of the following parameter combinations:

*r*

_{1}=

*r*

_{2}=

*ɛ*∈ {0.05, 0.10, …, 0.20},

*N*

_{2}∈ {11, 55, 99},

*N*

_{1}= 4

*π*, and

*D*

_{tr}=

*D*

_{rot}= 1. Recalling (35), the absorbing fractions thus range from

*f*

_{1}≈ 0.008 to

*f*

_{1}≈ 0.126 and

*f*

_{2}≈ 0.007 to

*f*

_{2}= 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

*N*

_{2}. The theoretical values of

*κ*are calculated from (35)–(37) using the numerically-calculated value of

*c*

_{0}from Sec. IV. In the bottom panel of Fig. 3, we plot the relative error,

*ɛ*for the three values of

*N*

_{2}. This plot shows excellent agreement between theoretical and numerically optimized trapping rate values.

In the top panel of Fig. 4, we plot the Kolmogorov–Smirnov distance $E(\kappa )$ in (47) (error in distribution) for these theoretical *κ* values as a function of *ɛ* for the three values of *N*_{2}. 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 *S*_{emp} and the theoretical survival probabilities *S*(*x*_{0}, *t*; *κ*) in (46) using the theoretical *κ* values for *ɛ* = 0.05 and the three values of *N*_{2}.

### C. Results for fixed *ɛ* and varying absorbing fractions *f*_{1}, *f*_{2}

*ɛ*→ 0, and thus, assumed that the fraction of the surface covered in patches,

*f*

_{1}, and the fraction of the particle covered in patches,

*f*

_{2}, are both small,

In Fig. 5, we fix *r*_{1} = *r*_{2} = *ɛ* = 0.1 and let *N*_{1}, *N*_{2} each take the values {11, 55, 99, 199, 299, 399} so that [recalling (35)] *f*_{1}, *f*_{2} each range from 2.75% up to 99.75% (we take *D*_{tr} = *D*_{rot} = 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(\kappa )$ in (47) (error in distribution) for the theoretical *κ* values.

*κ*and

*κ*

_{opt}) is fairly high even away from regime (48). The error in distribution is largest in the regimes

*κ*

_{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 *f*_{1} and *f*_{2} are both small (and simulations confirm its high accuracy in this regime), the theory is fairly accurate even in the intermediate regimes of (i) *f*_{1} ≪ 1 and *f*_{2} ≪̸ 1, *f*_{2} ≉ 1, (ii) *f*_{2} ≪ 1 and *f*_{1} ≪̸ 1, *f*_{1} ≉ 1, and (iii) *f*_{1} ≪̸ 1, *f*_{1} ≉ 1 and *f*_{2} ≪̸ 1, *f*_{2} ≉ 1. Furthermore, the regimes of (49) and (50) reduce to previously studied problems.

## VII. DISCUSSION

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 particles^{18–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, 36–42. This approach is related to the strong localized perturbation analysis initially developed in Refs. 43 and 44.

*R*

_{c}> 0, is much greater than the radius of the particle. That is, we require

*D*

_{rot}and (ii) the surface curvature becomes relevant on the timescale $Rc2/Dtr$, we thus require

*β*≫̸ 1 and

*β*≪̸ 1).

We also assumed that the absorbing fractions were both small (i.e., *f*_{1} ≪ 1, *f*_{2} ≪ 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 *f*_{1} ≪ 1, *f*_{2} ≪ 1 and the prior results for the regimes *f*_{1} ≈ 1, *f*_{2} ≪ 1 and *f*_{1} ≪ 1, *f*_{2} ≈ 1.

## ACKNOWLEDGMENTS

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

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

### APPENDIX A: DERIVATION OF *κ*_{1} IN (3) FOR A PATCHY SURFACE AND A PERFECTLY REACTIVE PARTICLE

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., *f*_{2} = 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*).

*t*> 0 and introduce an outer expansion, which is valid away from the patches,

*ɛ*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

*S*

_{1}satisfies

*S*

_{1}near the

*m*th trap on the left wall (0,

*y*

_{m},

*z*

_{m}), we introduce stretched coordinates for fixed

*m*∈ {1, 2, …},

*w*as

*w*satisfies

*w*=

*w*

_{0}+

*O*(

*ɛ*) into (A4) implies that

*w*

_{0}is harmonic,

*w*

_{0}satisfies the following boundary conditions:

^{45}For our purposes, we only need the far-field behavior of

*w*

_{0}, which is

*α*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,

*y*

_{m},

*z*

_{m}) must agree with the far-field behavior of the inner expansion as

*ρ*→ ∞. That is,

*α*= 1 and

*m*∈ {1, 2, …}, the boundary condition in (A2) becomes

*x*= 0. We now derive the effective trapping rate by determining the behavior of (A10) as

*ɛ*→ 0.

*ɛ*→ 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),

### APPENDIX B: DERIVATION OF DISTRIBUTIONAL FORM IN (29)

*f*satisfies

*x*= 0:

*f*as (

*x*,

*y*,

*z*,

*θ*,

*φ*) → (0,

*y*

_{m},

*z*

_{m},

*θ*

_{n},

*φ*

_{n}), define the inner solution analogous to (22),

*g*∼

*ɛ*

^{−3}

*g*

_{0}+⋯. By the same argument that yielded (24) and (25), we have that

*g*

_{0}is harmonic,

*g*

_{0}satisfies the following boundary condition at

*μ*= 0:

*δ*(

*αx*) =

*δ*(

*x*)/|

*α*|. Ignoring arbitrary additive constants, the solution to (B1) and (B2) is

^{26}

*g*with the near-field behavior of

*f*shows that

*f*indeed has the singular behavior in (28).

### APPENDIX C: SIMULATION ALGORITHM FOR THE ELECTROSTATIC CAPACITANCE *c*_{0}

*c*

_{0}. The method relies on a probabilistic representation of the solution

*w*

_{0}of (24) and (25). Let $Z(t)\u2208R5$ be a standard five-dimensional Brownian motion. Define the first time that this process reaches the duocylinder $D$ in (27),

*w*

_{0}satisfying (24) and (25) can be written as

*q*is the probability that

**Z**eventually reaches $D$, conditioned on the initial position of

**Z**,

*q*must be harmonic for

*μ*≠ 0 and satisfy the boundary conditions at

*μ*= 0,

*q*over the surface of the five-dimensional ball of radius $\rho =\mu 2+\nu 2+\eta 2+s2+p2>0$ centered at the origin. If

*μ*= 0 and

*ρ*> 0 is such that

*∂*

_{μ}

*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\u0304(\rho )$:

*w*

_{0}in (26) yield

Equation (C1) implies that we can calculate *c*_{0} by calculating the probability $q\u0304(\rho )$ 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\u0304(\rho )$ 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.

*μ*= 0, which is

^{12}

*U*is uniformly distributed on [0, 1]. Then, given

**Z**at position $(\mu ,\nu ,\eta ,s,p)\u2208R5$ at the beginning of stage I, the position at the end of stage I is

*ξ*

_{1},

*ξ*

_{2},

*ξ*

_{3},

*ξ*

_{4}are independent standard normal random variables.

*ν*,

*η*,

*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 $Z\u2209D$ and, therefore,

*b*

_{1}and

*b*

_{2}/

*β*are the two radii of the duocylinder in (27). If we define the distance

*d*≔ max{

*d*

_{1},

*d*

_{2}}, then the five-dimensional sphere of radius

*d*centered at $(0,\nu ,\eta ,s,p)\u2209D$ cannot intersect $D$. Therefore, stage II places the particle uniformly on the surface of this 5D sphere.

### APPENDIX D: SIMULATION ALGORITHM FOR THE FULL STOCHASTIC SYSTEM

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

*r*

_{1}> 0 with density

*N*

_{1}/(4

*π*) placed on a square lattice such that the centers are

*x*

_{0},

*y*

_{0},

*z*

_{0}) at

*t*= 0, where we fix starting distance from the left wall

*x*

_{0}and choose (

*y*

_{0},

*z*

_{0}) from a uniform distribution on $[\u22122\pi N1,2\pi N1]2$. The sphere has

*N*

_{2}(with

*N*

_{2}odd) reactive spherical caps of polar angle

*r*

_{2}> 0 with centers $(\theta \u0304k,\phi \u0304k)$ on a Fibonacci lattice,

^{20}

*θ*

_{0},

*ϕ*

_{0}) sampled from a uniform distribution on the unit sphere using

*ξ*

_{1},

*ξ*

_{2}) is uniformly distributed on [0, 1]

^{2}.

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

*et al.*

^{12}for the translational diffusion of the sphere plus rotational diffusion. Given a particle located at (

*x*

_{0},

*y*

_{0},

*z*

_{0}) 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

*t*

_{1}is (0,

*y*

_{1},

*z*

_{1}), where

*y*

_{1}and

*z*

_{1}are independently drawn from normal distributions with respective means

*y*

_{0}and

*z*

_{0}and variances 2

*D*

_{tr}

*t*

_{1}.

*θ*

_{1},

*φ*

_{1}), the orientation of the sphere following the propagation of the sphere to the wall, given time

*t*

_{1}. 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,

*f*(

*θ*,

*t*) from Ref. 46, which yields the following cumulative distribution function:

*θ*

_{d}in

*θ*, given

*t*

_{1}. 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}

*y*

_{1},

*z*

_{1}) 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

*d*

_{A}and

*d*

_{B}to the nearest patches on both the plane and the sphere, respectively, are calculated. The signed distance

*d*

_{A}to the nearest patch on the plane is given by

*d*

_{B}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 length

^{48}

*a*

_{i}to the center of the

*i*th patch and subtracting the patch radius

*r*

_{2}from the minimum arc length min

*a*

_{i},

*d*

_{A}< 0 and

*d*

_{B}< 0, then the trial is complete and the final elapsed time is recorded. If

*d*

_{A}> 0, the algorithm proceeds to case A of stage II. If

*d*

_{A}< 0 and

*d*

_{B}> 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, *y*_{1}, *z*_{1}) at some time *t* > 0 at the start of stage II, we propagate the particle to the hemisphere of radius *d*_{A} centered at (0, *y*_{1}, *z*_{1}), which ensures that the end condition could not have been met. We sample a time *t*_{2} from the well-known distribution for a three-dimensional diffusion process to exit a sphere of radius *d*_{A} (see, for example, Ref. 12). The time *t*_{2} is added to the cumulative elapsed time.

The particle is placed at (*x*_{2}, *y*_{2}, *z*_{2}), sampled uniformly on the hemisphere of radius *d*_{A} centered at (0, *y*_{1}, *z*_{1}) 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*(*θ*, *t*_{2}) in (D1), given *t*_{2} 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 (*x*_{2}, *y*_{2}, *z*_{2}) and (*θ*_{2}, *φ*_{2}) to (*x*_{0}, *y*_{0}, *z*_{0}) 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 *d*_{B} (the distance to the nearest patch on the particle).

*d*

_{B}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,

*d*

_{B}, 0 <

*θ*<

*d*

_{B}, and

*t*using the MATLAB pdepe function.

^{49}To sample an arrival time to the spherical cap of radius

*d*

_{B}, we sample a uniform random number

*ξ*∈ [0, 1] and numerically solve

*s*and then compute the exit time

*t*

_{2}=

*s*/

*D*

_{rot}. The time

*t*

_{2}is added to the elapsed time, and the angular position is updated from (

*θ*

_{1},

*φ*

_{1}) with rotational change (

*d*

_{B},

*φ*

_{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 (

*x*

_{2},

*y*

_{2},

*z*

_{2}) of the particle after time

*t*

_{2}is sampled by generating three standard normal random numbers (

*ξ*

_{1},

*ξ*

_{2},

*ξ*

_{3}) and computing

*x*

_{2},

*y*

_{2},

*z*

_{2}) and (

*θ*

_{2},

*φ*

_{2}) to (

*x*

_{0},

*y*

_{0},

*z*

_{0}) and (

*θ*

_{0},

*φ*

_{0}).

## REFERENCES

*The Visual Guide to Extra Dimensions: Visualizing the Fourth Dimension, Higher-Dimensional Polytopes, and Curved Hypersurfaces*

*Mixed Boundary Value Problems in Potential Theory*

*Quantum Computation and Quantum Information*