title: Surface Hopping Propagator: An Alternative Approach to Diffusion-Influenced Reactions
authors: Denis S. Grebenkov
text: 
  
  
  
  
   denis.grebenkov@polytechnique.edu
  
  Laboratoire de Physique de la Matière Condensée (UMR 7643), 
   
  CNRS – Ecole Polytechnique, IP Paris, 91128 Palaiseau, France
  
  
  
  
  
  
  
  
  Dynamics of a particle diffusing in a confinement can be seen a
  sequence of bulk-diffusion-mediated hops on the confinement surface.
  Here, we investigate the surface hopping propagator that describes the
  position of the diffusing particle after a prescribed number of
  encounters with that surface.  This quantity plays the central role in
  diffusion-influenced reactions and determines their most common
  characteristics such as the propagator, the first-passage time
  distribution, and the reaction rate.  We derive explicit formulas for
  the surface hopping propagator and related quantities for several
  Euclidean domains: half-space, circular annuli, circular cylinders,
  and spherical shells.  These results provide the theoretical ground
  for studying diffusion-mediated surface phenomena.  The behavior of
  the surface hopping propagator is investigated for both “immortal”
  and “mortal” particles.
  
  
  
  02.50.-r, 05.40.-a, 02.70.Rr, 05.10.Gg
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  Surface Hopping Propagator: 
   An Alternative Approach to Diffusion-Influenced Reactions
      Denis S. Grebenkov
      December 30, 2023
  =======================================================================================
  
  
  
  
  
  
  
  
  
  § INTRODUCTION
  
  
  In many natural phenomena, particles diffuse in a confinement towards
  its surface where they can react, permeate, relax their activity or be
  killed.  Examples include heterogeneous catalysis, permeation across
  cell membranes, filtering in porous media, surface relaxation in
  nuclear magnetic resonance, and animal foraging
  <cit.>.
  These phenomena are conventionally described by diffusion equation (or
  more general Fokker-Planck equation) with appropriate boundary
  conditions <cit.>.  In particular, most common
  properties of diffusion-influenced reactions are derived from the
  propagator G_q(,t|_0), i.e., the probability density of the
  event that a particle, started from a bulk point _0 at time 0,
  has not reacted on the surface and located at a bulk point  at
  time t.  For normal diffusion, the propagator satisfies the
  diffusion equation inside a confining domain Ω
  
  
      ∂_t G_q(,t|_0) = D  Δ_ G_q(,t|_0)    (∈Ω),
  
  subject to the initial condition G_q(,t=0|_0) = δ(-_0)
  and the Robin boundary condition on the boundary :
  
  
      - ∂_ G_q(,t|_0) = q   G_q(,t|_0)    (∈),
  
  where Δ_ is the Laplace operator acting on ,
  δ() is the Dirac distribution, and ∂_ is the
  normal derivative on the boundary  oriented outwards the domain
  Ω.
  The parameter q = κ/D is the ratio between the surface
  reactivity (or permeability, or relaxivity, etc.) κ and bulk
  diffusivity D.  In chemical physics, the Robin boundary condition
  was put forward by Collins and Kimball <cit.> and later
  explored by many researchers
  <cit.>
  (see an overview in <cit.>).  The major disadvantage of
  the conventional description is that the surface reactivity κ
  (or q) enters implicitly as a parameter of the Robin boundary
  condition (<ref>).
  
  In a recent work <cit.>, we proposed an alternative
  description based on the concept of boundary local time.  The boundary
  local time ℓ_t characterizes the fraction of time that a
  diffusing particle spends in a close vicinity of the reflecting
  boundary, as well as the number of encounters with that boundary
  <cit.>, see Eqs. (<ref>, <ref>) below.  This is
  a fundamental concept in the theory of stochastic processes
  <cit.>, which remains largely unknown and almost
  unemployed in physics, chemistry and biology.  To incorporate
  ℓ_t, we introduced the full propagator P(,ℓ,t|_0), i.e.,
  the joint probability density of finding a particle at point  at
  time t with its boundary local time ℓ, given that it started
  from _0 at time 0.  The crucial advantage of this alternative
  description is that P(,ℓ,t|_0) characterizes diffusion in
  confinement with reflecting (inert) boundary.  In turn, the
  surface reactivity is introduced via a stopping condition on the
  boundary local time.  In particular, we derived
  
  
      G_q(,t|_0) = ∫_0^∞ dℓ  e^-qℓ  P(,ℓ,t|_0) ,
  
  where the surface reactivity q appears explicitly as a
  parameter of the Laplace transform with respect to the boundary local
  time ℓ.  In this way, the single full propagator
  P(,ℓ,t|_0) describes the whole family of partially reactive
  surfaces (characterized by q).  Moreover, one can replace the
  exponential factor e^-qℓ by a more general function to
  implement other surface reaction mechanisms far beyond the
  conventional partial reactivity described by the Robin boundary
  condition (<ref>), see
  <cit.> for details.  In this light, the full propagator
  P(,ℓ,t|_0) turns out to be the intrinsic key quantity that
  describes all sorts of diffusion-mediated surface phenomena in a given
  confinement.
  
  A successful implementation of this new paradigm requires efficient
  methods for accessing the full propagator.  In <cit.>, the
  Laplace transform of the full propagator was expressed in terms of the
  so-called “surface hopping propagator” Σ_p(,ℓ|_0),
  i.e., the probability density of the event that a particle, started
  from a boundary point _0, has survived against a “bulk killing”
  with the rate p and located at a boundary point  at the boundary
  local time ℓ.  The rate p ≥ 0 accounts for eventual
  disappearance of the particles during its diffusion in the domain
  Ω due to a bulk reaction or spontaneous disintegration,
  relaxation, photobleaching or death.  In this scheme, one can consider
  both “mortal” (p >0) and “immortal” (p = 0) particles
  <cit.>.  In other words, the surface
  hopping propagator describes bulk-diffusion-mediated displacements
  between two encounters with the boundary, separated by the boundary
  local time ℓ.  The concept of such a surface exploration by
  successive hops through the bulk was formulated by Bychuk and
  O'Shaugnessy <cit.> and later confirmed by
  single-particle tracking experiments <cit.>.
  Former theoretical descriptions of surface hopping diffusion in terms
  of effective surface propagators were based on coupled bulk-surface
  diffusion equations with adsorption/desorption kinetics
  <cit.>.
  In turn, the surface hopping propagator Σ_p(,ℓ|_0) is a
  conceptually different quantity, which characterizes surface
  displacements not in terms of physical time t (as earlier) but
  in terms of the boundary local time ℓ (the number of encounters).
  To our knowledge, the surface hopping propagator, introduced in
  <cit.> as an efficient way to access the full propagator,
  is a new object, and the present paper aims at uncovering its
  properties.
  
  The paper is organized as follows.  In Sec. <ref>, we
  formulate the theoretical framework for diffusion-mediated surface
  phenomena, build an intuitive ground for the surface hopping
  propagator, and recall some general relations from <cit.>.
  Main results are reported in Sec. <ref>, in which the
  surface hopping propagator is computed and investigated for several
  domains.  Section <ref> summarizes and concludes the
  paper.
  
  
  
  
  
  
  § SURFACE HOPPING PROPAGATOR
  
  
  
  How many reflections does a particle undertake up to a given time t
  or during its lifetime?  Where is the particle after n reflections?
  For the common continuous-time Brownian motion, these natural
  questions have old but disappointing (trivial) answers.  In fact,
  Brownian motion crossing a smooth surface is known to return
  infinitely many times to that surface within an infinitely short time
  period <cit.>.  To get more satisfactory answers, one needs to
  reformulate these questions in a regularized way.  For instance, one
  can substitute Brownian motion by a sequence of independent jumps
  (e.g., a random walk on a lattice).  However, it is more convenient to
  keep considering continuous stochastic process _t but to introduce
  a thin surface layer of width a, _a = {∈Ω :  | -
  | < a}, and to count the number _t^a of crossings of this
  layer by reflected Brownian motion up to time t.  As a → 0, the
  number of crossings diverges but a _t^a converges to the random
  process, introduced by Lévy and called the boundary local time
  <cit.>:
  
  
      ℓ_t = lim_a→ 0 a _t^a .
  
  While the continuous time t represents the number of jumps of
  duration δ that the particle undertakes in the bulk, the
  boundary local time ℓ is the proxy of the number of encounters
  with the boundary (reflections of amplitude a).  Equivalently,
  ℓ_t is proportional to the fraction of time that a particle spent
  in the surface layer of width a up to time t:
  
  
      ℓ_t = lim_a→ 0D/a∫_0^t dt'  __a(_t') ,
  
  where the integral is the residence time of reflected Brownian motion
  _t in _a, and __a() is the indicator function of
  that layer: __a() = 1 for ∈_a, and 0
  otherwise.  While ℓ_t is historically called “local time”, it
  has units of length (see also <cit.>).
  In some definitions, the diffusion coefficient D is removed
  from Eq. (<ref>), yielding the boundary local time in units
  of time per length, i.e., the time spent in the surface layer rescaled
  by its width.  We also stress that the boundary local time should not
  be confused with a closely related notion of the point local time,
  i.e., a fraction of time spent in an infinitesimal vicinity of a fixed
  bulk point.  The latter was thoroughly investigated, in particular,
  for Brownian motion and Bessel processes (see
  <cit.> and references therein).
  
  
  The former two questions should thus be reformulated in terms of the
  boundary local time: How large is the boundary local time ℓ_t up
  to a given time t or during the lifetime of a particle?  Where is
  the particle after a boundary local time ℓ?  Answers to both
  these questions are given by the surface hopping propagator
  Σ_p(,ℓ|_0), as discussed below.
  
  
  
  
  
   §.§ Intuitive picture
  
  
  Before presenting main results for general domains, it is instructive
  to provide the motivation and intuition for the surface hopping
  propagator.  Let us consider a particle diffusing in the upper
  half-plane, Ω = ×_+, with reflecting horizontal axis
  = {(x,0)  :  x∈}.  If the particle started from a
  boundary point _0 = (x_0,0), its “next” encounter with the
  boundary would occur exactly at _0, as discussed above.  To
  overcome this problem, we introduce a thin surface layer of width a
  (Fig. <ref>).  Now, one can ask what is the position of the
  next encounter with the boundary after crossing the horizontal line y
  = a (i.e., after exiting from the surface layer).  This
  regularization eliminates too short Brownian trajectories that remain
  within the layer.  As the first crossing of the line y = a typically
  occurs near the starting point _0, one can move the starting point
  from (x_0,0) to (x_0,a) and then search for the probability
  density of the first arrival onto the horizontal axis.  This is the
  harmonic measure density, which for the upper half-plane takes the
  form of the Cauchy density,
  
  
      p_1(x|(x_0,a)) = a/π [(x-x_0)^2 + a^2] ,
  
  and can thus describe the first encounter position x after leaving
  the boundary from x_0 and crossing the surface layer of width a.
  After this encounter, the particle continues diffusion, independently
  of its past, so that the second encounter position is determined by
  the convolution of two Cauchy densities:
  
  
      p_2(x|(x_0,a))     = ∫_ dx_1  p_1(x|(x_1,a))   p_1(x_1|(x_0,a)) 
          = 2a/π [(x-x_0)^2 + (2a)^2] .
  
  Similarly, the position of the n-th encounter is determined by
  
  
      p_n(x|(x_0,a)) = na/π [(x-x_0)^2 + (na)^2] .
  
  In the limit a → 0 with any fixed n, this density converges to
  the Dirac distribution, p_n(x|(x_0,a)) →δ(x-x_0),
  illustrating the above statement that (reflected) Brownian motion
  returns infinitely many times to the first hitting point within an
  infinitely short period.  As the right-hand side of Eq. (<ref>)
  depends on a via the product na, a nontrivial result can only be
  obtained in the limit a→ 0 when na is fixed.  Setting ℓ =
  na, one obtains the surface hopping propagator for the upper
  half-plane:
  
  
      Σ_0(,ℓ|_0) = ℓ/π [(x-x_0)^2 + ℓ^2] ,
  
  with boundary points = (x,0) and _0 = (x_0,0).  As eventual
  death of the particle during its bulk diffusion was ignored, we set p
  = 0 in the subscript.
  
  
  
  
  
  While the above construction can be performed in any confining domain,
  its practical realization involves numerous convolutions of the
  harmonic measure density which in general are difficult to compute
  (the above explicit computation was possible due to the explicit form
  of the Cauchy distribution and its specific “infinite divisibility”
  property, i.e., the invariance of its form upon convolutions).  In the
  next subsection, we present a general approach to access the surface
  hopping propagator.
  
  
  
  
  
  
  
  
   §.§ General approach
  
  
  We consider a particle diffusing in an Euclidean domain Ω⊂^d with a smooth boundary .  In <cit.>, the
  surface hopping propagator Σ_p(,ℓ|_0) (with p ≥ 0)
  was shown to be the kernel of the semi-group exp(-_p ℓ)
  generated by the Dirichlet-to-Neumann operator _p.  This is a
  pseudo-differential self-adjoint operator acting on functions on the
  boundary  (see rigorous definitions and mathematical details in
  <cit.>).  For a given
  function f on , this operator associates another function g
  on , _p  :  f→ g = (∂_ w)_|, where w
  satisfies the Dirichlet boundary value problem:
  
  
      (p - D Δ) w = 0   in Ω,      w_| = f
  
  (if Ω is unbounded, one also needs to impose the regularity
  condition: w→ 0 as ||→∞, see below).  For instance, if
  f describes a concentration of particles maintained on the boundary
  , then _p f = (∂_ w)_| is proportional to the
  diffusive flux of these particles into the bulk.  Note that there is a
  family of Dirichlet-to-Neumann operators parameterized by p (or
  p/D).  As the kernel of the semi-group exp(-_pℓ), the
  surface hopping propagator satisfies
  
  
      ∂_ℓΣ_p(,ℓ|_0) = - _p Σ_p(,ℓ|_0),
  
  subject to the initial condition Σ_p(,ℓ=0|_0) = δ(
  - _0).  This equation resembles the diffusion equation
  (<ref>), in which the physical time t is replaced by the
  boundary local time ℓ, and the Laplace operator Δ is
  replaced by -_p.
  
  When the boundary  of the domain is bounded, the
  Dirichlet-to-Neumann operator has a discrete spectrum, i.e., a
  countable set of positive eigenvalues μ_n^(p) and eigenfunctions
  v_n^(p)() forming a complete orthonormal basis in the space
  L_2():
  
  
      _p   v_n^(p)() =  μ_n^(p)  v_n^(p)() .
  
  The surface hopping propagator admits thus the spectral expansion:
  
  
      Σ_p(,ℓ|_0) = ∑_n [v_n^(p)(_0)]^*   v_n^(p)()   e^-μ_n^(p)ℓ ,
  
  where asterisk denotes complex conjugate.  In other words, finding
  Σ_p(,ℓ|_0) is equivalent to studying the spectral
  properties of the Dirichlet-to-Neumann operator _p.
  
  When p > 0, all eigenvalues μ_n^(p) are strictly positive, and
  the surface hopping propagator vanishes as ℓ increases.  This is
  a direct consequence of bulk reaction that may lead to eventual death
  or disappearance of the diffusing particle during its motion.  In
  particular, 
  
  
      ∫_ d Σ_p(,ℓ|_0) < 1       (ℓ > 0),
  
  i.e., this density is not normalized to 1, in the same way as the
  conventional propagator G_q(,t|_0) is not normalized to 1 in
  the presence of reactive boundary (q > 0).  
  
  In turn, the case p = 0 is more subtle.  For restricted diffusion in
  a bounded domain, the ground eigenfunction is constant, v_0^(0)()
  = ||^-1/2, whereas the associated eigenvalue is zero:
  μ_0^(0) = 0.  Due to the orthogonality of other eigenfunctions
  to v_0^(0), the surface hopping propagator is normalized to 1:
  
  
      ∫_ d Σ_0(,ℓ|_0) = 1       (ℓ > 0).
  
  In contrast, if diffusion is transient, all eigenvalues μ_n^(0)
  are strictly positive, and the normalization is lost again, here, due
  to the possibility of escaping at infinity.  This is the case of
  diffusion in the exterior of a bounded domain in ^d with d ≥
  3 (for d = 2, see a short discussion in Sec. <ref>).
  
  
  
  
  
  
  
  
   §.§ Relation to other quantities
  
  
  
  As shown in <cit.>, the surface hopping propagator opens
  the door to access most common diffusion-reaction characteristics such
  as the full propagator, the conventional propagator, the first-passage
  time distribution, and the reaction rate.  In particular, the Laplace
  transform of the full propagator reads
  
  
      P̃(,ℓ,p|_0) = G̃_∞(,p|_0)  δ(ℓ) 
          + ∫_ d_0
      ∫_ d j̃_∞(,p|)  Σ_p(,ℓ|_0)/D j̃_∞(_0,p|_0) ,
  
  where G_∞(,t|_0) is the propagator for perfectly reactive
  boundary (with Dirichlet boundary condition
  G_∞(,t|_0)_| = 0),
  
  
      j_∞(,t|_0) = - D(∂_ G_∞(,t|_0))_=
  
  is the probability flux density on that boundary, and tilde denotes
  the Laplace transform with respect to time t, e.g.,
  
  
      P̃(,ℓ,p|_0) = ∫_0^∞ dt   e^-pt  P(,ℓ,t|_0).
  
  Substituting the spectral expansion (<ref>) into
  Eq. (<ref>), one also gets
  
  
      P̃(,ℓ,p|_0)     = G̃_∞(,p|_0)  δ(ℓ) 
          + 1/D∑_n [V_n^(p)(_0)]^*   V_n^(p)()   e^-μ_n^(p)ℓ ,
  
  with
  
  
      V_n^(p)() = ∫_ d j̃_∞(,p|)   v_n^(p)().
  
  
  In turn, the full propagator determines most common quantities of
  diffusion-influenced reactions, in particular, the conventional
  propagator via Eq. (<ref>).  Moreover, one gets the marginal
  probability density of the boundary local time ℓ_t (see also
  <cit.>):
  
  
      ρ(ℓ,t|_0) = ∫_Ω d  P(,ℓ,t|_0),
  
  and the probability density of the first-crossing time T_ℓ =
  inf{ t>0  : ℓ_t > ℓ} of a level ℓ by ℓ_t:
  
  
      U(ℓ,t|_0) = D ∫_ d  P(,ℓ,t|_0).
  
  The latter determines the probability density of the conventional
  first-passage time to a partial reactive boundary as:
  
  
      H_q(t|_0) = ∫_0^∞ dℓ  q  e^-qℓ  U(ℓ,t|_0).
  
  In the Laplace domain, one can use the spectral expansion
  (<ref>) to write
  
  
      Ũ(ℓ,p|_0) = ∑_n [V_n^(p)(_0)]^*   e^-μ_n^(p)ℓ∫_ d  v_n^(p)().
  
  Inverting the Laplace transform, one gets U(ℓ,t|_0) and thus
  gains access via Eq. (<ref>) to the whole family of
  probability densities H_q(t|_0).
  
  The surface hopping propagator also determines the spread harmonic
  measure density, ω_q(|_0), which characterizes the boundary
  point on a partially reactive boundary, at which reaction occurs
  <cit.>.  For a particle started
  from _0, one has
  
  
      ω_q(|_0) = ∫_ d_0  ω_q(|_0)  j̃_∞(_0,0|_0),
  
  where
  
  
      ω_q(|_0) = ∫_0^∞ dℓ  q   e^-qℓ Σ_0(,ℓ|_0) .
  
  More generally, the Laplace transform of q Σ_p(,ℓ|_0)
  with respect to the boundary local time ℓ yields the probability
  density of the reaction point  on the boundary in the presence of
  bulk reactions with the rate p.
  
  
  
  
  
  
  
  
   §.§ Extension
  
  
  
  In the above construction of the surface hopping propagator
  Σ_p(,ℓ|_0) and the full propagator P(,ℓ,t|_0),
  the boundary local time ℓ is counted on the whole boundary 
  of the confining domain Ω.  This boundary local time is then
  used to incorporate surface reactions, like in Eq. (<ref>).
  In certain applications, however, only a subset of the boundary,
  Γ⊂, is reactive, whereas the remaining part
  \Γ is inert and just passively confines the
  diffusing particle inside the domain.  This is a typical case of an
  escape through a hole Γ, or of a target Γ surrounded by
  a reflecting surface.  As encounters of the particle with the passive
  part of the boundary do not matter, one needs to count the boundary
  local time only on the reactive part Γ.
  
  An extension to this setting is straightforward.  In fact, one can
  re-define ℓ_t through the residence time in a close vicinity of
  the reactive part: Γ_a = {∈Ω :  | - Γ| <
  a}:
  
  
      ℓ_t = lim_a→ 0D/a∫_0^t dt'  _Γ_a(_t') .
  
  The associated surface hopping propagator can be constructed as
  earlier by modifying the definition of the Dirichlet-to-Neumann
  operator.  In fact, as the passive part of the boundary,
  \Γ, is irrelevant, one can define the
  Dirichlet-to-Neumann operator _p^Γ, acting on functions on
  Γ as _p^Γ :  f → g = (∂_ w)_|Γ,
  where w is the solution of the mixed Dirichlet-Neumann boundary
  value problem:
  
  
      (p - DΔ) w = 0   in Ω,     {[       w_|Γ = f,; (∂_ w)_|\Γ = 0.;                 ].
  
  Here, the Neumann boundary condition on \Γ
  implements explicitly the reflecting character of the passive part of
  the boundary.
  
  In the following, we will present several examples of the surface
  hopping propagator for a subset of the boundary.
  
  
  
  
  
  
  
  
  
  
  
  
  § EXPLICIT RESULTS FOR SEVERAL CONFINING DOMAINS
  
  
  
  In this section, we illustrate the properties of the surface hopping
  propagator and related quantities for several confining domains, for
  which the eigenbasis of the Dirichlet-to-Neumann operator can be
  explicitly derived.  Even though these spectral properties are known
  to experts, we will provide some clarifications to guide readers.  We
  start with the half-space, for which all quantities, including the
  full propagator, will be derived in closed explicit forms.  Then we
  consider two-dimensional circular annuli between two concentric
  circles that include as limiting cases the interior and the exterior
  of a disk.  An extension of these results to three-dimensional
  circular cylinders is briefly presented.  Similarly, we discuss
  spherical shells between two concentric spheres that include the
  interior and the exterior of a ball.  A numerical computation of
  the eigenmodes of the Dirichlet-to-Neumann operator in non-concentric
  perforated spherical domains is discussed in
  <cit.>.  Technical derivations are moved to
  Appendices.
  
  
  
  
  
   §.§ Half-space
  
  
  
  It is instructive to start with the case of the half-space Ω =
  { = (x_1,…,x_d)∈^d :  x_d > 0}, for which all
  quantities of interest can be obtained in a closed analytic form.
  Even though the boundary  of the half-space is not bounded, the
  derived formulas can be adapted.  In this case, the spectrum of the
  Dirichlet-to-Neumann operator is continuous, and sums over eigenmodes
  should be replaced by integrals.  Moreover, the solutions of the
  eigenvalue equation (<ref>), v_n^(p), are not
  L_2()-normalized and thus cannot be called “eigenfunctions”.
  Nevertheless, we will keep this term for v_n^(p), bearing in mind
  its limitations.
  
  To clarify the ideas, we start again with the planar case (d = 2).
  The translational symmetry of the boundary implies that v_n(s) =
  e^ins/√(2π) are the “eigenfunctions” of the
  Dirichlet-to-Neumann operator _p.  In fact, since w(x,y) = e^inx
  - y √(n^2 + p/D) satisfies the modified Helmholtz equation
  (<ref>) in the upper half-plane, e^ins is an
  “eigenfunction”, associated to the “eigenvalue” μ_n^(p) =
  √(n^2 + p/D).  Note that this “eigenfunction” does not depend
  on p due to the above symmetry.  The prefactor 1/√(2π) comes
  from the orthogonality of “eigenfunctions”:
  
  
      ∫_-∞^∞ ds   v_n(s)   [v_n'(s)]^* = δ(n-n').
  
  Skipping technical details, we formally rewrite the spectral expansion
  (<ref>) of the surface hopping propagator as
  
  
      Σ_p(s,ℓ|s_0) = ∫_-∞^∞ dn   [v_n(s_0)]^*   v_n(s)  exp(-μ_n^(p)ℓ) ,
  
  where the former summation index n now takes real values in .
  As a consequence, we get
  
  
      Σ_p(s,ℓ|s_0)     = ∫_-∞^∞dn/2π  e^in(s_0 - s) -  ℓ√(n^2 + p/D)
          = ℓ/π (ℓ^2 + (s-s_0)^2) ζ K_1(ζ) ,
  
  with 
  
  
      ζ = √(p/D) √(ℓ^2 + (s-s_0)^2) ,
  
  and K_ν(z) is the modified Bessel function of the second kind.  At
  p = 0, we retrieve the Cauchy density (<ref>).  In this
  case, the surface exploration up to the boundary local time ℓ is
  equivalent to the first arrival onto that surface of Brownian motion
  started from (s_0,ℓ), i.e. the distance ℓ above the surface.
  Figure <ref> illustrates the behavior of the surface
  hopping propagator.  Changing progressively the boundary local time
  ℓ, one observes the spreading of the surface hopping propagator.
  
  For “immortal” particles (p = 0), the surface hopping propagator
  exhibits heavy tails, Σ_0(s,ℓ|s_0) ∝ |s-s_0|^-2, in
  particular, the variance of the arrival point s is infinite (the
  mean is s_0 due to the symmetric form of this propagator).  Such
  displacements with infinite variance resemble Lévy flights
  <cit.>.  This is the consequence of unbounded exploration
  region that allows for very long and far-reaching trajectories.  The
  situation is drastically different for “mortal” particles (p > 0),
  for which long trajectories are penalized by tiny chances of survival.
  In fact, the central part of this distribution (when ℓ^2 +
  |s-s_0|^2 ≪ D/p) resembles again the Cauchy density (Fig.
  <ref>(b)), which is however truncated by exponential
  tails at large s (when ζ≫ 1): Σ_p(s,ℓ|s_0) ∝
  e^-ζ∼ e^-|s-s_0|√(p/D).  The bulk rate p (or, more
  precisely, p/D) controls this truncation.  As a consequence, all the
  positive moments of the arrival boundary point are finite.
  Expectedly, the surface hopping propagator is not normalized to 1
  for p > 0:
  
  
      ∫_ ds  Σ_p(s,ℓ|s_0) = e^-ℓ√(p/D) .
  
  Interestingly, even the conditional surface hopping propagator that
  accounts only for the survived particles after renormalization by
  e^-ℓ√(p/D), shows an exponential decay with s.
  
  
  
  
  
  
  
  In higher dimensions, the “eigenfunctions” and “eigenvalues” of
  the Dirichlet-to-Neumann operator have similar form:
  
  
      v_() = e^i(·)/(2π)^(d-1)/2 ,     μ_^(p) = √(||^2 + p/D),
  
  with the “multi-index” = (n_1,…,n_d-1) ∈^d-1.  The surface propagator can thus be written as
  
  
      Σ_p(,ℓ|_0) = ∫_^d-1d/(2π)^d-1  e^i· (_0 -) - ℓ√(||^2 + p/D) .
  
  In spherical coordinates, the integral over all orientations gives
  
  
      Σ_p(,ℓ|_0)     = |-_0|^3-d/2/(2π)^d-1/2∫_0^∞ dn   n^d-1/2 J_d-3/2(n|-_0|) 
         × e^-ℓ√(n^2 + p/D)
          = Σ_0(,ℓ|_0)  ζ^d/2  K_d/2(ζ)/Γ(d/2) 2^d/2-1 ,
    
  where ζ is given by Eq. (<ref>), and
  
  
      Σ_0(,ℓ|_0) = Γ(d/2)/π^d/2 ℓ/(ℓ^2 + |-_0|^2)^d/2
  
  is again the harmonic measure density on the hyperplane ^d-1
  (see also <cit.>), and J_ν(z) is the Bessel function
  of the first kind.  In the right-hand side of
  Eq. (<ref>), one can recognize the Laplace-transformed
  probability flux density j̃_∞(,p|(_0,ℓ)) onto a
  perfectly absorbing hyperplane from the bulk point _0 = (_0,ℓ).
  
  The inverse Laplace transform with respect to p yields then
  
  
      Ł^-1_t{Σ_p(,ℓ|_0)}    = j_∞(,t|(_0,ℓ)) 
          = exp(-| - _0|^2/4Dt)/(4π D t)^(d-1)/2 ℓ exp(-ℓ^2/4Dt)/√(4π Dt^3) .
  
  We emphasize that this relation is specific to the case of the
  half-space.
  
  Using this relation and the representation (<ref>), one
  can easily compute the double integral over _1 and _2 in
  Eq. (<ref>) to get
  
  
      P̃(,ℓ,p|_0) = G̃_∞(,p|_0) δ(ℓ) + Σ_p(,z+z_0+ℓ|_0)/D ,
  
  where = (,z) and _0 = (_0,z_0).  The inverse Laplace
  transform with respect to p yields then
  
  
      P(,ℓ,t|_0) = G_∞(,t|_0) δ(ℓ) + j_∞(,t|(_0,z+z_0+ℓ))/D
          = exp(-|-_0|^2/4Dt)/(4π Dt)^d/2{( e^-(z-z_0)^2/(4Dt) - e^-(z+z_0)^2/(4Dt)) δ(ℓ) 
          + ℓ + z + z_0/Dt e^-(z+z_0+ℓ)^2/(4Dt)}.
  
  This is the explicit exact form of the full propagator for the
  half-space.
  
  One can easily check that the integral over ∈^d_+ gives the
  marginal probability density of the boundary local time ℓ_t:
  
  
      ρ(ℓ,t|_0) = (z_0/√(4Dt)) δ(ℓ) + exp(- (z_0+ℓ)^2/4Dt)/√(π D t) .
  
  This expression does not depend on the dimension d and the lateral
  coordinate _0 of the starting point _0, given that the
  boundary local time is independent of lateral displacements and
  determined by the transverse motion (on the half-line).  The
  distribution of the boundary local time ℓ_t was studied in
  <cit.>.  
  
  In turn, the integral of Eq. (<ref>) over ℓ∈_+
  yields the marginal probability density of the position, i.e.,
  conventional propagator G_0(,t|_0) in the half-space with
  reflecting boundary:
  
  
      G_0(,t|_0)     = exp(-|-_0|^2/4Dt)/(4π Dt)^d/2
         ×( e^-(z-z_0)^2/(4Dt) + e^-(z+z_0)^2/(4Dt)).
  
  Moreover, with the general expression (<ref>), one retrieves
  the propagator G_q(,t|_0) for reactive boundary:
  
  
      G_q(,t|_0) = exp(-|-_0|^2/4Dt)/(4π Dt)^d/2{ e^-(z-z_0)^2/(4Dt) + 
          e^-(z+z_0)^2/(4Dt)(1 - 2q√(π Dt) (z+z_0/√(4Dt) + q√(Dt)))},
  
  where (z) = e^z^2(z) is the scaled complementary error
  function.  Expectedly, the propagators in Eqs. (<ref>,
  <ref>) exhibit translational invariance along 
  coordinate and are factored into the lateral Gaussian (free)
  propagator and the transverse propagator on the half-line
  (0,∞).  Note that its integral over the arrival point 
  yields the survival probability
  
  
      S_q(t|_0) = (z_0/√(4Dt)) + 
      e^-z_0^2/4Dt(z_0/√(4Dt) + q√(Dt)),
  
  which does not depend on _0 and coincides with the survival 
  probability for the semi-axis with the partially reactive endpoint.
  The classical expression for the associated probability density of 
  the reaction time is then retrieved:
  
  
      H_q(t|_0)     = qD e^-z_0^2/(4Dt){1/√(π Dt)
          - q (z_0/√(4Dt) + q√(Dt)) } .
  
  Finally, Eq. (<ref>) yields
  
  
      U(ℓ,t|_0) = (ℓ+z_0) e^-(ℓ+z_0)^2/(4Dt)/√(4π D t^3) ,
  
  i.e., we retrieved the classical formula for the probability density
  of the first crossing time of a level ℓ by the boundary local
  time of reflected Brownian motion on the half-line (see, e.g.,
  <cit.>).
  
  
  To complete this section, we briefly mention that the above
  computations can be easily extended to a slab domain between parallel
  hyperplanes, one of which is reflecting (see <cit.> for
  more details).  In other words, one can consider Ω =
  {∈^d  :  0 < x_d < L } and study the Dirichlet-to-Neumann
  operator on the hyperplane x_d = 0 in the presence of the reflecting
  hyperplane x_d = L.  The “eigenfunctions” of such
  Dirichlet-to-Neumann operator remain unchanged due to the
  translational symmetry, whereas the “eigenvalues” are
  
  
      μ_^(p) = √(||^2 + p/D) tanh(L √(||^2 + p/D)).
  
  As L→∞, one retrieves the former case of the half-space.  The
  former integral representations for the surface hopping propagator and
  related quantities remain valid if Eq. (<ref>) is used for
  μ_^(p).  In contrast, the presence of tanh(z) prevents
  from getting simple closed formulas for these quantities in the case
  of a slab.
  
  
  
  
  
  
  
   §.§ Circular annuli and spherical shells
  
  
  In Appendices <ref> and <ref>, we provide
  explicit formulas for the eigenfunctions and eigenvalues of the
  Dirichlet-to-Neumann operator in several rotationally invariant
  domains: circular annuli, the interior and the exterior of a disk,
  circular cylinders, spherical shells, the interior and the exterior of
  a ball.  These formulas allow one to get the surface hopping
  propagator Σ_p(,ℓ|_0) via the spectral expansion
  (<ref>), as well as the full propagator P(,ℓ,t|_0)
  and all the related quantities, as discussed in Sec. <ref>.
  From these basic results, one can thoroughly investigate various
  diffusion-mediated surface phenomena in the above domains.  In this
  paper, we keep our focus on the surface hopping propagator and
  illustrate its properties for these domains.
  
  
  
  Figure <ref>(a,b) shows the surface hopping propagator
  for the interior of a disk of radius R.  At p = 0 (panel (a)), the
  surface hopping propagator Σ_0(s,ℓ|s_0) coincides with the
  harmonic measure density on the circle, see Eq. (<ref>).
  Expectedly, it evolves from the Dirac distribution at ℓ = 0 to
  the uniform distribution 1/(2π R) as ℓ→∞.  When p >
  0 (panel (b)), eventual death of the particle during its bulk
  explorations affects this propagator.  At small ℓ, the particle
  spends short time in the bulk so that there is almost no effect of the
  bulk rate p (at moderate p): blue curves on (a) and (b) panels are
  almost identical.  As ℓ increases, the effect of p becomes more
  prominent.  At large ℓ, the particle has enough time to explore
  the interior of the disk, leading again to the uniform distribution of
  the arrival point.  This can also be seen from the spectral expansion
  (<ref>), in which the contribution from higher-order
  eigenmodes of the Dirichlet-to-Neumann operator vanished.  The surface
  hopping propagator becomes almost flat again, Σ_p(s,ℓ|s_0)
  ≃ e^-μ_0^(p)ℓ/(2π R), but its level is now attenuated
  by bulk reactions.  As μ_0^(p=0.1)≈ 0.0494, the
  attenuation factor e^-μ_0^(p)ℓ is not strong even at ℓ
  = 10 in the considered example.  However, one can still see a
  qualitative difference between p = 0 and p = 0.1 cases: in the
  former case, the curves approach the limit 1/(2π R), whereas in
  the latter case, they are progressively shifted downward.
  
  Figure <ref>(c,d) shows the surface hopping propagator
  for the exterior of a disk of radius R.  We note that
  Eq. (<ref>) for the surface hopping propagator at p =
  0 remains the same for diffusion outside the disk.  This reflects the
  conformal invariance of the harmonic measure with respect to an
  inversion mapping of the interior of the disk to its exterior.  In
  contrast, for p > 0, the behavior of the surface hopping propagator
  is different for diffusion inside and outside the disk, especially for
  large ℓ.  In fact, the particle diffusing outside the disk
  undertakes much longer bulk excursions between successive encounters
  with the boundary and thus has higher chances to be killed by a bulk
  reaction.  This leads to much smaller values of the surface hopping
  propagator (we recall that the propagator is not normalized to 1 for
  p > 0).
  
  
  
  
  
  Figure <ref> presents the surface hopping propagator
  Σ_0(s,ℓ|s_0) for a circular annulus with both reactive
  circles of radii R and L.  As the boundary is composed of two
  circles, there are two distinct choices of the starting point s_0
  (shown by red dot): either on the inner circle of radius R, or on
  the outer circle of radius L (and one can set s_0 = 0 in both
  cases due to rotational invariance).  For convenience of presentation,
  the curvilinear coordinate s runs here from -2π R to 0 for the
  inner circle, and from 0 to 2π L for the outer circle.  On the
  horizontal axis, s is rescaled by R for negative values and by L
  for positive values so that the horizontal axis varies from -2π
  and 2π.  Let us first consider the particle started on the inner
  circle (Fig. <ref>(a)).  At small ℓ (blue
  curves), this particle does not move far away from the starting point
  s_0, so that Σ_0(s,ℓ|s_0) rapidly decays when s/R varies
  from 0 to -π (its later increase for s/R ranging from -π
  to -2π is due to the symmetry).  Similarly, bulk excursions of the
  particle rarely terminate at the outer circle so that the surface
  hopping propagator remains small on that boundary (the range of
  positive s).  Clearly, the minimum corresponds to the boundary point
  s = π L on the outer boundary which is located behind the starting
  point s_0.  As the boundary local time increases, the particle
  explores further boundary regions, both on the inner and outer
  circles.  In the limit ℓ→∞, the surface hopping
  propagator approaches the uniform density, 1/(2π(R+L)), on both
  circles, as expected.
  
  When the starting point is on the outer circle
  (Fig. <ref>(b)), the picture is very similar, i.e.,
  the particle remains on the outer circle (close to the starting point)
  at small ℓ but then spreads away.  Note that the stronger decay
  of the surface hopping propagator is caused by the fact that the outer
  circle is much longer than the inner one.  For this reason, a large
  truncation order was needed to accurately compute
  Σ_0(s,ℓ|s_0) in this case.
  
  For “mortal” particles (not shown), eventual death due to the bulk
  rate p > 0 penalizes long trajectories, as in the case of diffusion
  inside a disk.  Moreover, in highly reactive media, the particle has
  tiny chances to move from one circle to the other, and these parts of
  the boundary become decoupled.  In other words, as there is almost no
  survived particles that crossed the annulus, the properties of the
  boundary far away from the starting point do not matter.
  
  
  
  
  
  
  
  
  
  
  Figure <ref> shows similar results for the interior and
  exterior of a ball.  As previously, Σ_p(,ℓ|_0) behaves
  differently for diffusion inside and outside the ball.  However, this
  difference is considerably enhanced in three dimensions due to the
  recurrent versus transient character of diffusion.  In fact, the
  particle diffusing outside a ball can escape to infinity with a finite
  probability.  As a consequence, the surface hopping propagator is not
  normalized to 1 here even for p = 0.  Actually, the normalization
  constant is e^-ℓ/R, i.e., the escape probability is 1 -
  e^-ℓ/R.  On panel (c), one can see that the surface hopping
  propagator becomes again uniform (as the contribution of higher-order
  eigenfunctions vanishes) but attenuated by the factor e^-ℓ/R.
  
  
  
  
  
  Finally, Fig. <ref> illustrates the surface hopping
  propagator on a spherical target of radius R, surrounded by an outer
  reflecting sphere of radius L.  This setting is qualitatively in
  between the interior and the exterior of a ball.  On one hand, as this
  domain is bounded and diffusion is recurrent, the surface hopping
  propagator evolves towards the uniform density 1/(4π R^2), as for
  the interior case.  On the other hand, the outer reflecting sphere is
  located relatively far from the target and thus allows for long
  trajectories, as for the exterior case.  For p = 0, the panel (a) of
  Fig. <ref> resembles the panel (a) of
  Fig. <ref>, even so diffusion occurs in different
  regions in these two settings.  In contrast, when p = 0.1, the
  panel (b) of Fig. <ref> is much closer to the panel
  (d) of Fig. <ref>.  In fact, in both cases, the
  particle is allowed to undertake long trajectories between successive
  encounters with the target.
  
  
  
  
  
  
  
  
  
  
  
  § DISCUSSION AND CONCLUSION
  
  
  
  In this paper, we investigated the properties of the surface hopping
  propagator Σ_p(,ℓ|_0) recently introduced in
  <cit.>.  This is a conceptually new quantity that
  describes bulk-diffusion-mediated exploration of a surface.  In
  contrast to former works
  <cit.>,
  which relied on coupled bulk-surface diffusion equations and aimed to
  characterize the position of the particle on a surface after some
  physical time t, here we operate with the boundary local time
  ℓ, which is a proxy of the number of encounters with that
  surface.  The surface hopping propagator turns out to be dual to
  the conventional propagator G_q(,t|_0).  In fact, as
  G_q(,t|_0) characterizes displacements between bulk points
  _0 and  in physical time t (i.e., after a number of bulk
  jumps), Σ_p(,ℓ|_0) characterizes effective displacements
  between boundary points _0 and  in boundary local time ℓ
  (i.e., after a number of reflections on the boundary).  While
  G_q(,t|_0) is the semi-group of the Laplace operator -Δ
  (acting in the bulk), Σ_p(,ℓ|_0) is the semi-group of
  the Dirichlet-to-Neumann operator _p (acting on the boundary).  In
  this light, the spectral expansion (<ref>) is dual to the
  spectral expansion for the conventional propagator:
  
  
      G_q(,t|_0) = ∑_n [u_n^(q)(_0)]^*   u_n^(q)()   e^-λ_n^(q) t ,
  
  where λ_n^(q) and u_n^(q)() are the eigenvalues and
  L_2(Ω)-normalized eigenfunctions of the diffusion operator
  -DΔ:
  
  
  
      -D Δ u_n^(q)    = λ_n^(q)  u_n^(q)     (∈Ω), 
      ∂_ u_n^(q) + q   u_n^(q)    = 0      (∈) ,
  
  
  where we highlighted the dependence on the reactivity parameter q
  through the Robin boundary condition (<ref>).  Similarity
  and duality of Eqs. (<ref>, <ref>) are
  remarkable.  We recall that the spectral expansion (<ref>)
  is valid for a bounded domain Ω, whereas the spectral expansion
  (<ref>) is valid for a bounded boundary .  As a
  consequence, Eq. (<ref>) seems to be more general, as it is
  also valid when Ω is the exterior of a bounded domain, for
  which Eq. (<ref>) is not applicable anymore.  In spite of
  this considerable advantage, the Dirichlet-to-Neumann operator and its
  eigenbasis were not earlier employed to describe diffusion-influenced
  reactions and other diffusion-mediated surface phenomena.  The present
  paper, along with Refs. <cit.>, aim to shift
  the theoretical description of these phenomena towards a new
  fundamental ground.  Using the formulas derived in this paper, one can
  access directly not only the surface hopping propagator, but all the
  related quantities, including the full propagator, first-passage time
  distribution, the survival probability and the reaction rate.
  
  
  
  
  
  
  
  
  The author acknowledges a partial financial support from the Alexander von
  Humboldt Foundation through a Bessel Research Award. 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § CIRCULAR ANNULUS
  
  
  
  Here, we consider a circular annulus between two concentric circles of
  radii R < L: Ω = {∈^2  :  R < || < L}.  There are
  four possible combinations of surface reactivity: (i) both circles are
  reactive, (ii) the inner circle is reflecting while the outer circle
  is reactive, (iii) the inner circle is reactive while the outer circle
  is reflecting, and (iv) both circles are reflecting.  As surface
  reaction is not possible in the last case, it is excluded.  The first
  case corresponds directly to our general setting when the whole
  boundary is reactive.  However, as this case involves two disjoint
  parts of the boundary (the inner and the outer circles), its analysis
  is the most complicated.  For this reason, we start with the case
  (ii), then briefly discuss the case (iii), and finally give the
  solution for the case (i).
  
  
  
  
  
  
   §.§ Reactive outer circle
  
  
  In order to determine the spectrum of the Dirichlet-to-Neumann
  operator associated with the outer circle Γ = {∈^2  : 
  || = L}, one needs to solve the mixed boundary value problem
  (<ref>).  Its general solution can be searched in polar
  coordinates (r,ϕ) as:
  
  
      w(r,ϕ) = ∑_n=-∞^∞ c_n   e^inϕ  g_n(r),
  
  where c_n are unknown coefficients (to be fixed by the boundary
  condition), 
  
  
      g_n(r) = K'_n(α R) I_n(α r) - I'_n(α R) K_n(α r)/K'_n(α R) I_n(α L) - I'_n(α R) K_n(α L) .
  
  are radial functions with α = √(p/D), prime denotes the
  derivative with respect to the argument, and I_n(z) is the modified
  Bessel functions of the first kind.  One can easily check that
  g'_n(R) = 0 by construction.  For convenience, we have chosen a
  particular normalization g_n(L) = 1.  As the normal derivative is
  equal to the radial derivative, the action of the Dirichlet-to-Neumann
  operator _p does not affect the angular part.  In other words, the
  rotational symmetry of this domain implies that Fourier harmonics are
  the eigenfunctions of _p, defined on the outer circle Γ:
  
  
      v_n(s) = e^in s/L/√(2π L)     (n ∈),
  
  where the curvilinear coordinate s is related to the polar angle
  ϕ as s = ϕ L.  The associated eigenvalues are
  
  
      μ_n^(p) = g'_n(L)       (n ∈).
  
  The eigenfunctions do not depend on p, whereas the eigenvalues are
  twice degenerate, except for n = 0.  Here, the index n runs over
  all integer numbers.  In the limit p→ 0, one gets
  
  
      μ_n^(0) = |n|/L 1 - (R/L)^2|n|/1 + (R/L)^2|n| .
  
  
  The spectral decomposition (<ref>) of the surface hopping
  propagator reads
  
  
      Σ_p(s,ℓ|s_0) = ∑_n=-∞^∞e^in(s_0-s)/L/2π Lexp( - μ_n^(p)ℓ).
  
  
  
  To access the full propagator, one also needs to compute
  V_n^(p)(_0) from Eq. (<ref>).  Using the summation
  formulas from <cit.>, the Laplace-transformed propagator
  G̃_∞(,p|_0) and thus j̃_∞(,p|_0)
  for a circular annulus with Dirichlet boundary condition on the outer
  circle and Neumann boundary condition on the inner circle read
  
  
      G̃_∞(,p|_0)    = 1/2π D∑_n=-∞^∞ e^in(ϕ-ϕ_0)  g_n(r_0) 
         ×[K_n(α L) I_n(α r) - I_n(α L) K_n(α r)]  , 
      j̃_∞(s,p|_0)    = 1/2π L∑_n=-∞^∞ e^in(ϕ-ϕ_0)  g_n(r_0) ,
  
  where = (r,ϕ) and _0 = (r_0,ϕ_0) in polar coordinates,
  R ≤ r_0 ≤ r ≤ L, s = ϕ L, and we used the Wronskian
  
  
      I'_n(z) K_n(z) - K'_n(z) I_n(z) = 1/z .
  
  The projection of j̃_∞(,p|_0) onto an eigenfunction
  v_n(s) from Eq. (<ref>) reads then
  
  
      V_n^(p)(_0) = v_n(ϕ_0)   g_n(r_0)  .
  
  The orthogonality of Fourier harmonics to a constant function reduces
  Eq. (<ref>) to
  
  
      Ũ(ℓ,p|_0) =  g_0(r_0) exp(-μ_0^(p)ℓ)  ,
  
  from which Eq. (<ref>) gives the Laplace-transformed
  probability density of the reaction time as
  
  
      H̃_q(p|_0) = g_0(r_0)  1/1 + μ_0^(p)/q .
  
  
  
  
  
  
  
  
  
  
   §.§ Interior of a disk
  
  
  In the limit R→ 0, the inner boundary shrinks to a point, and one
  gets the solution for the interior of a disk of radius L: Ω =
  {∈^2 :  || < L} (see also <cit.>).  The radial
  functions become
  
  
      g_n(r) = I_n(r√(p/D))/I_n(L√(p/D)) ,
  
  while the eigenvalues and eigenfunctions of the Dirichlet-to-Neumann
  operator are still given by Eqs. (<ref>,
  <ref>).  Other expressions are also valid; in particular, 
  Eqs. (<ref>, <ref>) are applicable.
    
  At p = 0, Eq. (<ref>) yields μ_n^(0) = |n|/L
  and thus
  
  
      Σ_0(s,ℓ|s_0) = 1 - e^-2ℓ/L/2π L(1 - 2cos(s-s_0/L) e^-ℓ/L + e^-2ℓ/L) .
  
  As expected, this propagator evolves from the Dirac distribution
  δ(s-s_0) at ℓ = 0 to the uniform distribution 1/(2π L)
  as ℓ→∞.  Setting ρ = L e^-ℓ/L, one can recognize
  in this form the Poisson kernel in the disk of radius L.  The
  Poisson kernel describes the harmonic measure density on the disk,
  i.e., the probability density of the first arrival onto the circle of
  radius R at point (L,s/L) for Brownian motion started from a point
  (ρ,s_0/L) (written in polar coordinates).  As in the planar case
  discussed in Sec. <ref>, the distribution of the position of
  the diffusing particle at the boundary local time ℓ (i.e., after
  a prescribed number of encounters with the reflecting circle) is
  identical to the distribution of the first arrival point on the fully
  absorbing circle, where the boundary local time ℓ determines the
  starting point in the latter setting.  Figure
  <ref>(a,b) illustrates the behavior of the surface
  hopping propagator.
  
  
  
  
  
  
  
  
  
  
   §.§ Reactive inner circle
  
  
  We briefly discuss the case (iii) when the inner circle is reactive
  and surrounded by a reflecting outer circle.  This is a typical
  setting of a small reactive target confined in a domain surrounded by
  an outer reflecting boundary <cit.>.  Here, we search for
  the spectrum of the Dirichlet-to-Neumann operator associated with the
  inner circle: Γ = {∈^2  :  || = R}.  Repeating the
  above construction step by step, one realizes that the eigenfunctions
  are the Fourier harmonics on the inner circle
  
  
      v_n(s) = e^in s/R/√(2π R)     (n ∈),
  
  where the curvilinear coordinate s is related to the polar angle
  ϕ as s = ϕ R.  The associated eigenvalues are
  
  
      μ_n^(p) = - g'_n(R)       (n ∈),
  
  where sign minus appears due to the direction of the normal
  derivative, ∂_ = -∂_r, and
  
  
      g_n(r) = K'_n(α L) I_n(α r) - I'_n(α L) K_n(α r)/K'_n(α L) I_n(α R) - I'_n(α L) K_n(α R)
  
  are the radial functions satisfying g_n(R) = 1 and g'_n(L) = 0.
  In the limit p→ 0, one gets
  
  
      μ_n^(0) = |n|/R 1 - (R/L)^2|n|/1 + (R/L)^2|n| .
  
  
  
  The spectral decomposition (<ref>) of the surface hopping
  propagator reads
  
  
      Σ_p(s,ℓ|s_0) = ∑_n=-∞^∞e^in(s_0-s)/R/2π Rexp( - μ_n^(p)ℓ).
  
  
  The Laplace-transformed propagator G̃_∞(,p|_0) and
  thus j̃_∞(,p|_0) for a circular annulus with
  Dirichlet boundary condition on the inner circle and Neumann boundary
  condition on the outer circle read
  
  
      G̃_∞(,p|_0)    = 1/2π D∑_n=-∞^∞ e^in(ϕ-ϕ_0)  g_n(r_0) 
         ×[K_n(α R) I_n(α r) - I_n(α R) K_n(α r)]  , 
      j̃_∞(s,p|_0)    = 1/2π R∑_n=-∞^∞ e^in(ϕ-ϕ_0)  g_n(r_0) ,
  
  where = (r,ϕ), _0 = (r_0,ϕ_0), R ≤ r ≤ r_0 ≤
  L, s = ϕ R.  In turn, Eqs. (<ref>,
  <ref>, <ref>) remain unchanged.  Some
  first-passage properties in this setting were studied in
  <cit.>.
  
  
  
  
  
  
  
   §.§ Exterior of a disk
  
  
  
  In the limit L→∞, the outer boundary is pushed away to
  infinity, and one deals with diffusion in the exterior of a disk of
  radius R: Ω = {∈^2  :  || > R}.  In this limit, the
  radial functions from Eq. (<ref>) are reduced to
  
  
      g_n(r) = K_n(r√(p/D))/K_n(R√(p/D)) .
  
  The eigenvalues and eigenfunctions of the Dirichlet-to-Neumann
  operator are still given by Eqs. (<ref>,
  <ref>), with g_n(r) from Eq. (<ref>).
  Other earlier expressions are as well applicable; in particular,
  Eqs. (<ref>, <ref>) are valid.
  Figure <ref>(c,d) illustrates the behavior of the
  surface hopping propagator.
  
  The exterior of a disk presents a convenient example to illustrate
  subtle points of recurrent diffusion outside a planar bounded domain.
  As the radial functions g_n(r) in Eq. (<ref>) for p >
  0 vanish exponentially fast in the limit r→∞, a general
  solution w(r,ϕ) of the modified Helmholtz equation
  (<ref>) also vanishes, in agreement with the regularity
  condition.  In the case p = 0, the radial functions for n  0
  become g_n(r) = (R/r)^|n| and vanish again.  However, the limit of
  g_0(r) as p→ 0 is equal to 1 that does not vanish at infinity,
  thus violating the regularity condition.  This is a consequence of the
  simple fact that the rotationally invariant Laplace equation in the
  plane, w” + 1/r w' = 0, has a general solution c_1 +
  c_2ln r that does not vanish as r→∞, except for the trivial
  choice c_1 = c_2 = 0.  This is a well-known problem for such planar
  domains, for which, in particular, there is no steady-state reaction
  rate <cit.>.  In the remaining part of the
  paper, we do not discuss this subtle case.
  
  
  
  
  
  
  
  
  
  
   §.§ Both reactive circles
  
  
  When both circles are reactive, one needs to consider the
  Dirichlet-to-Neumann operator on the whole boundary composed of two
  disjoint circles: Γ_1 = {∈^2  :  || = R} and
  Γ_2 = {∈^2  :  || = L}.  A general solution of
  Eq. (<ref>) can be searched in the form
  
  
      w(r,ϕ) = ∑_n=-∞^∞(c_n,1  g_n,1(r) + c_n,2  g_n,2(r)) e^inϕ,
  
  where the unknown coefficients c_n,1 and c_n,2 are set by
  boundary conditions, and the radial functions
  
  
      g_n,1(r)     = K_n(α L) I_n(α r) - I_n(α L) K_n(α r)/K_n(α L) I_n(α R) - I_n(α L) K_n(α R) , 
      
      g_n,2(r)     = K_n(α R) I_n(α r) - I_n(α R) K_n(α r)/K_n(α R) I_n(α L)  - I_n(α R) K_n(α L)
  
  satisfy g_n,1(L) = 0, g_n,1(R) = 1, and g_n,2(R) = 0,
  g_n,2(L) = 1 for convenience (other linear combinations of
  I_n(α r) and K_n(α r) could also be used).  Note that
  g_n,1(r) monotonously decreases, whereas g_n,2(r)
  monotonously increases on the interval (R,L).
  
  As the inner and outer circles are concentric, one can expect that an
  eigenfunction of the operator _p can be written as 
  
  
      v_n^(p)(s) = {[      a_n^(p)  e^in ϕ s ∈Γ_1    (ϕ = s/R),;      b_n^(p)  e^in ϕ s ∈Γ_2    (ϕ = s/L),;                      ].
  
  where a_n^(p) and b_n^(p) are some coefficients.  In fact, as
  the space of functions on the whole boundary = Γ_1 ∪Γ_2 is the direct product of spaces of functions on the inner
  (Γ_1) and outer (Γ_2) circles, the eigenfunction
  v_n^(p)(s) can be thought of being composed of two components.
  Here, these two components are proportional to e^inϕ due to the
  rotational symmetry (note that this claim can be shown rigorously by
  representing each component as a Fourier series and then using the
  orthogonality of Fourier harmonics).  Substituting such v_n^(p)(s)
  as the Dirichlet condition to Eq. (<ref>), one gets its
  solution as
  
  
      w(r,ϕ) = (a_n^(p)   g_n,1(r) + b_n^(p)  g_n,2(r)) e^inϕ,
  
  while its normal derivative on  reads
  
  
      _p v_n     = (∂_n w)_|
          = {[ - (a_n^(p)  g'_n,1(R) + b_n^(p)  g'_n,2(R)) e^in ϕ                                            s ∈Γ_1,;   (a_n^(p)  g'_n,1(L) + b_n^(p)  g'_n,2(L)) e^in ϕ                                            s ∈Γ_2.;                                                    ].
  
  If the right-hand side is proportional to v_n^(p), then
  v_n^(p) is indeed an eigenfunction of _p.  In other words, we
  get two equations
  
  
      {[ - (a_n^(p)  g'_n,1(R) + b_n^(p)  g'_n,2(R)) = μ_n^(p)  a_n^(p)  ,;    (a_n^(p)  g'_n,1(L) + b_n^(p)  g'_n,2(L)) = μ_n^(p)  b_n^(p) ,;                                                                   ].
  
  where the proportionality coefficient μ_n^(p) is the associated
  eigenvalue, and we used the particular form of functions g_n,1(r)
  and g_n,2(r).  These equations can be written in a matrix form as
  
  
      ( [ - g'_n,1(R) - g'_n,2(R);   g'_n,1(L)   g'_n,2(L);             ])
      ( [ a_n^(p); b_n^(p);         ]) = μ_n^(p)( [ a_n^(p); b_n^(p);         ]) .
  
  Solving this eigenvalue problem for the 2× 2 matrix on the
  left-hand side, one determines the eigenvalue μ_n^(p), as well
  as one of the coefficients (e.g., a_n^(p)).  The other coefficient
  (e.g., b_n^(p)) is fixed by imposing the L_2()-normalization
  of the eigenfunction v_n^(p).
  
  The Wronskian (<ref>) yields
  
  
      g'_n,1(L) = A/L ,      g'_n,2(R) = - A/R  ,
  
  where
  
  
      A = 1/K_n(α L) I_n(α R) - I_n(α L) K_n(α R) .
  
  The eigenvalue μ_n^(p) is then obtained as a solution of the
  quadratic equation 
  
  
      μ^2 + μ( g'_n,1(R) - g'_n,2(L) ) + B = 0 ,
  
  where 
  
  
      B     = - g'_n,1(R) g'_n,2(L) + g'_n,1(L) g'_n,2(R) 
          = - α^2 K'_n(α L) I'_n(α R) - I'_n(α L) K'_n(α R)/K_n(α L) I_n(α R) - I_n(α L) K_n(α R) .
  
  With the help of Eq. (<ref>), it is easy to check that the
  determinant of this equation is positive, so that there are two real
  roots.  Moreover, as B > 0, both roots are positive:
  
  
      μ_n,±^(p)    = 1/2( g'_n,2(L) - g'_n,1(R)  
         ±√((g'_n,1(R) + g'_n,2(L))^2 - 4 g'_n,1(L) g'_n,2(R))) .
  
  As a consequence, for each index n, there are two distinct
  eigenmodes.  For each of them, the coefficients are determined by the
  corresponding eigenvalue.  To avoid round-off errors in practical
  implementation, it is convenient to use slightly different (but
  formally equivalent) representation for + and - modes.  In fact,
  we set
  
  
      a_n,+^(p) = - c_n,+^(p) (g'_n,2(L) - μ_n,+^(p)),      b_n,+^(p) = c_n,+^(p)  g'_n,1(L)  ,
  
  where c_n,+^(p) is fixed by the normalization of the
  eigenfunction:
  
  
      1     = ∫_ ds   |v_n,+^(p)(s)|^2 
          = [c_n,+^(p)]^2 (2π R [g'_n,2(L) - μ_n,+^(p)]^2 + 2π L [g'_n,1(L)]^2 ).
  
  In turn, for the mode with μ_n,-^(p), one can use
  
  
      a_n,-^(p) = c_n,-^(p)  g'_n,2(R) ,      b_n,-^(p) = - c_n,-^(p) (g'_n,1(R) + μ_n,-^(p))  ,
  
  with
  
  
      1 = [c_n,-^(p)]^2 ( 2π R   [g'_n,2(R)]^2 + 2π L [g'_n,1(R) + μ_n,-^(p)]^2).
  
  
  
  In the limit p→ 0, one has
  
  
      g'_n,1(R)    → - |n|/R 1+(R/L)^2|n|/1-(R/L)^2|n| , 
      
      g'_n,1(L)    → - |n|/L 2(R/L)^|n|/1-(R/L)^2|n| , 
      
      g'_n,2(R)    →|n|/R 2(R/L)^|n|/1-(R/L)^2|n| , 
      
      g'_n,2(L)    →|n|/L 1+(R/L)^2|n|/1-(R/L)^2|n| ,
  
  so that
  
  
      μ_n,±^(0)    = |n|/2 1+γ_n/1-γ_n(1/L + 1/R. 
         . ±√((1/L + 1/R)^2 - 4(1-γ_n)^2/LR(1+γ_n)^2))  ,
  
  where γ_n = (R/L)^2|n|.  In the case n = 0, one can take
  the limit n→ 0 to get
  
  
      μ_0,-^(0) = 0,      μ_0,+^(0) = 1/L + 1/R/ln(L/R) .
  
  We also get
  
  
      a_0,-^(0)    = b_0,-^(0) = 1/√(2π(R+L)) , 
      
      a_0,+^(0)    = 1/√(2π R(1+R/L)) ,      b_0,+^(0) = -1/√(2π L(1 + L/R)) .
  
  
  Using the asymptotic behavior of the modified Bessel functions, one
  can check that 
  
  
      lim_R→ 0μ_n,-^(p)    = √(p/D) I'_n(L√(p/D))/I_n(L√(p/D)) , 
      lim_R→ 0μ_n,+^(p)    = +∞ , 
      lim_L→∞μ_n,-^(p)    = - √(p/D) K'_n(R√(p/D))/K_n(R√(p/D)) ,  
      lim_L→∞μ_n,+^(p)    = √(p/D) ,
  
  where we used that K'_n(z)/K_n(z) ≤ - 1.  In the limit R→ 0,
  μ_n,-^(p) approach the eigenvalues of _p for the interior
  of a disk of radius L, whereas μ_n,+^(p) diverge and thus do
  not contribute.  The opposite limit L→∞ is more subtle:
  μ_n,+^(p) approach the eigenvalues of _p for the exterior
  of a disk of radius R; however, μ_n,-^(p) accumulate near
  √(p/D).
  
  Finally, the Dirichlet propagator in the Laplace domain is
  
  
  
      G̃_∞(,p|_0) = 1/2π D∑_n=-∞^∞ e^in(ϕ-ϕ_0)  g_n,2(r_0) 
         ×[K_n(α L) I_n(α r) - I_n(α L) K_n(α r)]     (r_0 ≤ r),  
          19mm = 1/2π D∑_n=-∞^∞ e^in(ϕ-ϕ_0)  g_n,1(r_0) 
         ×[K_n(α R) I_n(α r) - I_n(α R) K_n(α r)]     (r ≤ r_0),
  
  
  where = (r,ϕ) and _0 = (r_0,ϕ_0).  The probability flux
  density reads then
  
  
  
      j̃_∞(s,p|_0)    = 1/2π R∑_n=-∞^∞ e^in(ϕ-ϕ_0)  g_n,1(r_0)     (s∈Γ_1) , 
         = 1/2π L∑_n=-∞^∞ e^in(ϕ-ϕ_0)  g_n,2(r_0)     (s∈Γ_2) .
  
  
  As a consequence, one gets
  
  
      V_n^(p)(_0) = (a_n^(p)  g_n,1(r_0) + b_n^(p)  g_n,2(r_0) ) e^inϕ_0
  
  that gives access to the full propagator P(,ℓ,t|_0).
  
  
  
  
  
  
   §.§ Cylindrical domains
  
  
  The above analysis can also be extended to cylindrical domains.  Let
  us first consider an infinite cylinder, Ω = Ω_0 ×,
  where Ω_0 is a disk of radius L.  As the boundary  is
  unbounded, the spectrum of the Dirichlet-to-Neumann operator _p is
  not discrete anymore.  Nevertheless, the symmetries of this domain
  admit the separation of variables and allow for getting
  “eigenfunctions” and “eigenvalues” in cylindrical coordinates as
  
  
  
      v_nk(ϕ,z)     = e^inϕ + ikz/2π √(L)     (n∈,  k∈) , 
      μ_nk^(p)    = α I'_n(α L)/I_n(α L) ,   α = √(p/D + k^2) .
  
  
  The surface hopping propagator reads then
  
  
      Σ_p(,ℓ|_0) = ∑_n=-∞^∞e^in(ϕ_0-ϕ)/2π L∫_dk/2π e^ik(z_0-z)  e^-μ_nk^(p)ℓ ,
  
  where = (L,ϕ,z) and _0 = (L,ϕ_0,z_0) in cylindrical
  coordinates.  If Ω_0 is a circular annulus, one has to use the
  appropriate radial function g_n(r), with α = √(p/D +
  k^2).  Other related quantities can also be obtained.
  
  When the cylinder is finite, Ω = Ω_0 × (0,b), the
  spectrum of _p is discrete again, but the analysis is more subtle.
  In fact, as in the case of a circular annulus, different combinations
  of reactivity patterns are possible: all the boundary is reactive;
  only the lateral boundary is reactive but the top and bottom disks are
  reflecting; only the top disk is reactive but the remaining boundary
  is reflecting; etc.  When only one part of the boundary is reactive,
  the analysis is rather simple.  For instance, if only the lateral
  boundary is reactive, the eigenfunctions and eigenvalues are 
  
  
  
      v_nk(ϕ,z)     = e^inϕ/√(2π L) √(2-δ_n,0)/√(b)cos(π k z/b) , 
      μ_nk^(p)    = α I'_n(α L)/I_n(α L) ,   α = √(p/D + (π k/b)^2) ,
  
  
  with n∈ and k=0,1,2,….  However, the analysis is more
  involved when the whole boundary is reactive.  
  
  
  
  
  
  
  
  
  § SPHERICAL SHELL
  
  
  
  In three dimensions, one can consider a spherical shell between two
  concentric spheres of radii R < L: Ω = {∈^3  :  R <
  || < L}, with three combinations of boundary conditions.  As the
  analysis is rather similar to the two-dimensional setting, the results
  are presented in a concise form.
  
  
  
  
  
  
   §.§ Reactive outer sphere
  
  
  As previously, we start with the case of the Dirichlet-to-Neumann
  operator _p associated with the reactive outer sphere Γ = {∈^3  :  || = L}.  The rotational invariance implies that the
  eigenfunctions of _p are the (normalized) spherical harmonics,
  
  
      v_nm() = 1/L  Y_mn(θ,ϕ)     (n=0,1,2,…,   |m| ≤ n).
  
  The eigenvalues are obtained by solving the mixed boundary value
  problem (<ref>):
  
  
      μ_n^(p) = g'_n(R)       (n=0,1,2,…),
  
  where
  
  
      g_n(r) = k'_n(α R) i_n(α r) - i'_n(α R) k_n(α r)/k'_n(α R) i_n(α L) - i'_n(α R) k_n(α L) ,
  
  α = √(p/D), and 
  
  
      i_n(z)     = √(π/(2z))  I_n+1/2(z) , 
      
      k_n(z)     = √(2/(π z))  K_n+1/2(z)
  
  are the modified spherical Bessel functions of the first and second
  kind, respectively.  The eigenfunctions do not depend on p, whereas
  the eigenvalues μ_n^(p) do not depend on the second index m
  and are thus (2n+1) times degenerate.  In the limit p→ 0, one
  gets
  
  
      μ_n^(0) = n(n+1)/L 1 - (R/L)^2n+1/n + 1 + n (R/L)^2n+1 .
  
  
  The surface hopping propagator from Eq. (<ref>) reads
  
  
      Σ_p(,ℓ|_0) = 1/L^2∑_n=0^∞∑_m=-n^n Y_mn^*(θ_0,ϕ_0)   
      Y_mn(θ,ϕ)   e^-μ_n^(p)ℓ .
  
  Since the eigenvalues do not depend on the index m, one can apply
  the addition theorem for spherical harmonics to evaluate the sum over
  m:
  
  
      Σ_p(,ℓ|_0) = 1/4π L^2∑_n=0^∞ (2n+1) P_n((·_0)/||  |_0|)  e^-μ_n^(p)ℓ ,
  
  where P_n(z) are Legendre polynomials.  
  
  
  
  
  One also needs to compute V_n^(p)(_0) from Eq. (<ref>).
  Using the summation formulas from <cit.>, the
  Laplace-transformed quantities G̃_∞(,p|_0) and thus
  j̃_∞(,p|_0) for a spherical shell with Dirichlet
  boundary condition on the outer sphere and Neumann boundary condition
  on the inner sphere read
  
  
      G̃_∞(,p|_0)    = ∑_n=0^∞α (2n+1)/4π D  P_n((·_0)/||   |_0|)
      g_n(r_0)  
         ×[k_n(α L) i_n(α r) - i_n(α L) k_n(α r)] ,
      j̃_∞(,p|_0)    = ∑_n=0^∞2n+1/4π L^2  P_n((·_0)/||   |_0|)
      g_n(r_0) ,
  
  where r = ||, r_0 = |_0|, R ≤ r_0 ≤ r ≤ L, and we
  used the Wronskian
  
  
      i'_n(z) k_n(z) - k'_n(z) i_n(z) = 1/z^2  .
  
  The projection of j̃_∞(,p|_0) onto an eigenfunction
  v_nm() from Eq. (<ref>) reads then
  
  
      V_nm^(p)(_0) = v_mn(θ_0,ϕ_0)   g_n(r_0) ,
  
  where _0 = (r_0,θ_0,ϕ_0) in spherical coordinates.  The
  orthogonality of spherical harmonics reduces Eq. (<ref>) to
  
  
      Ũ(ℓ,p|_0) =  g_0(r_0)  exp(-μ_0^(p)ℓ)  ,
  
  while the probability density of the reaction time reads
  
  
      H̃_q(p|_0) = g_0(r_0)  1/1 + μ_0^(p)/q .
  
  
  
  
  
  
  
  
  
  
  
   §.§ Interior of a ball
  
  
  In the limit R→ 0, the inner boundary shrinks to a point, and one
  gets the solution for the interior of a ball of radius L: Ω =
  {∈^3  :  || < L}, with radial functions
  
  
      g_n(r) = i_n(r√(p/D))/i_n(L√(p/D)) .
  
  The eigenvalues and eigenfunctions of the Dirichlet-to-Neumann
  operator are still given by Eqs. (<ref>,
  <ref>), and other earlier expressions remain valid;
  in particular, Eqs. (<ref>, <ref>) are
  applicable.  At p = 0, the eigenvalues are simply μ_n^(0) =
  n/L, for which Eq. (<ref>) can be evaluated explicitly
  using the generating function of Legendre polynomials:
  
  
      Σ_0(,ℓ|_0) = L/4π 1 - e^-2ℓ/L/|e^-ℓ/L_0 - |^3 .
   
  This expression coincides with the harmonic measure density on the
  sphere when the starting point is _0 e^-ℓ/L.  It can also be
  written in terms of the angle θ between vectors _0 and
  :
  
  
      Σ_0(,ℓ|_0) = 1/4π L^2 1 - e^-2ℓ/L/[1 - 2e^-ℓ/Lcosθ + e^-2ℓ/L]^3/2 .
  
  Figure <ref>(a,b) illustrates the behavior of the
  surface hopping propagator. 
  
  The orthogonality of spherical harmonics reduces Eq. (<ref>) to
  
  
      Ũ(ℓ,p|_0) = R/r_0 sinh(r_0√(p/D))/sinh(R√(p/D)) exp(-μ_0^(p)ℓ) ,
  
  with μ_0^(p) = √(p/D) (R√(p/D)) - 1/R.
  Similarly, one has
  
  
      H̃_q(p|_0) = R/r_0sinh(r_0 √(p/D))/sinh(R √(p/D)) 1/1 + μ_0^(p)/q ,
  
  from which the inverse Laplace transform yields the standard spectral
  expansion for H_q(t|_0).
  
  
  
  
  
  
  
  
  
  
   §.§ Reactive inner sphere
  
  
  The analysis for the reactive inner sphere is very similar.  The
  eigenfunctions of the Dirichlet-to-Neumann operator associated with
  the inner sphere Γ = {∈^3  :  || = R} are again the
  spherical harmonics but with the prefactor 1/R for a proper
  normalization:
  
  
      v_nm() = 1/R  Y_mn(θ,ϕ)     (n=0,1,2,…,   |m| ≤ n).
  
  The eigenvalues are 
  
  
      μ_n^(p) = -g'_n(R)       (n=0,1,2,…),
  
  where
  
  
      g_n(r) = k'_n(α L) i_n(α r) - i'_n(α L) k_n(α r)/k'_n(α L) i_n(α R) - i'_n(α L) k_n(α R) .
  
  In the limit p→ 0, one gets
  
  
      μ_n^(0) = n(n+1)/R 1 - (R/L)^2n+1/n + (n+1) (R/L)^2n+1 .
  
  The expression for the surface hopping propagator is almost identical
  to Eq. (<ref>):
  
  
      Σ_p(,ℓ|_0) = 1/4π R^2∑_n=0^∞ (2n+1) P_n((·_0)/||  |_0|)  e^-μ_n^(p)ℓ .
  
  
  
  The Laplace-transformed propagator G̃_∞(,p|_0) and
  thus j̃_∞(,p|_0) for a spherical shell with
  Dirichlet boundary condition on the inner sphere and Neumann boundary
  condition on the outer sphere read
  
  
      G̃_∞(,p|_0)    = ∑_n=0^∞α (2n+1)/4π D  P_n((·_0)/||   |_0|)
      g_n(r_0)  
         ×[k_n(α R) i_n(α r) - i_n(α R) k_n(α r)] ,
      j̃_∞(,p|_0)    = ∑_n=0^∞2n+1/4π R^2  P_n((·_0)/||   |_0|)
      g_n(r_0) ,
  
  where r = ||, r_0 = |_0|, R ≤ r ≤ r_0 ≤ L, from
  which
  
  
      V_nm^(p)(_0) = v_mn(θ_0,ϕ_0)   g_n(r_0) .
  
  These quantities determine the full propagator P(,ℓ,t|_0).
  
  
  
  
  
  
  
   §.§ Exterior of a ball
  
  
  In the limit L→∞, the outer boundary is pushed away to
  infinity, and one deals with diffusion in the exterior of a ball of
  radius R: Ω = {∈^3  :  ||> R}.  The radial
  functions are reduced to
  
  
      g_n(r) = k_n(r√(p/D))/k_n(R√(p/D)) ,
  
  while the eigenvalues and eigenfunctions of the Dirichlet-to-Neumann
  operator are still given by Eqs. (<ref>,
  <ref>).  Interestingly, the eigenvalues are just
  polynomials of √(p/D), e.g., μ_0^(p) = (1 +
  R√(p/D))/R.  The above expressions are as well applicable; in
  particular, Eqs. (<ref>, <ref>) are valid.
  
  At p = 0, the eigenvalues in Eq. (<ref>) are
  simplified as μ_n^(0) = (n+1)/R, and the surface hopping
  propagator can be computed explicitly as
  
  
      Σ_0(,ℓ|_0) = R   e^-ℓ/R/4π 1 - e^-2ℓ/R/|e^-ℓ/R_0 - |^3 .
   
  If L is replaced by R, this expression coincides with
  Eq. (<ref>), except for an extra factor e^-ℓ/R
  that accounts for the possibility of escaping to infinity.  Figure
  <ref>(c,d) illustrates the behavior of the surface
  hopping propagator.  
  
  The orthogonality of spherical harmonics reduces Eq. (<ref>) to
  
  
      Ũ(ℓ,p|_0) = R/r_0 exp(-(r_0-R+ℓ)√(p/D) - ℓ/R),
  
  from which the inverse Laplace transform yields
  
  
      U(ℓ,t|_0) = R e^-ℓ/R/r_0 r_0-R+ℓ/√(4π Dt^3) e^-(r_0 - R+ℓ)^2/(4Dt) .
  
  This is a rare example when the probability density U(ℓ,t|_0)
  is found in a simple closed form.  Setting ℓ = 0, one retrieves
  the probability density of the first-passage time for a perfectly
  absorbing sphere <cit.>.  In turn, the integral
  (<ref>) yields the probability density of the first-passage
  time to a partially reactive sphere <cit.>
  
  
      H_q(t|_0)     = qD/r_0 e^-(r_0-R)^2/(4Dt){R/√(π Dt)
          - (1 + qR) (r_0-R/√(4Dt) + (1+qR) √(Dt)/R) }.
  
  
  
  
  
  
  
  
  
  
  
   §.§ Both reactive spheres
  
  
  Finally, the analysis for both reactive spheres is the most involved
  but very similar to the planar case.  For this reason, we just
  reproduce the main formulas adapted to the three-dimensional case.
  Here, one employs two families of radial functions,
  
  
      g_n,1(r)     = k_n(α L) i_n(α r) - i_n(α L) k_n(α r)/k_n(α L) i_n(α R) - i_n(α L) k_n(α R) , 
      
      g_n,2(r)     = k_n(α R) i_n(α r) - i_n(α R) k_n(α r)/k_n(α R) i_n(α L)  - i_n(α R) k_n(α L) ,
  
  which satisfy g_n,1(R) = 1, g_n,1(L) = 0 and g_n,2(L) = 1,
  g_n,2(R) = 0.  The eigenfunctions are searched in the form
  
  
      v_nm^(p)() = {[ a_n^(p)  Y_mn(θ,ϕ)             ∈Γ_1 ,; b_n^(p)  Y_mn(θ,ϕ)             ∈Γ_2 ,;                    ].
  
  where Γ_1 and Γ_2 are the inner and the outer spheres
  forming the boundary .  The Wronskian (<ref>)
  implies
  
  
      g'_n,1(L) = A/α L^2 ,      g'_n,2(R) = - A/α R^2  ,
  
  where
  
  
      A = 1/k_n(α L) i_n(α R) - i_n(α L) k_n(α R) .
  
  In this case, one also has
  
  
      B     = - g'_n,1(R) g'_n,2(L) + g'_n,1(L) g'_n,2(R) 
          = - α^2 k'_n(α L) i'_n(α R) - i'_n(α L) k'_n(α R)/k_n(α L) i_n(α R) - i_n(α L) k_n(α R) .
  
  Using these expressions, one deduces again Eqs. (<ref>,
  <ref>) for the eigenvalue μ_n,±^(p) and the
  coefficients a_n,±^(p) and b_n,±^(p).  The
  normalization coefficient c_n,±^(p) is fixed by normalization:
  
  
      1     = [c_n,+^(p)]^2 (4π R^2 [g'_n,2(L) - μ_n,+^(p)]^2 + 4π L^2 [g'_n,1(L)]^2 ) , 
      
      1     = [c_n,-^(p)]^2 (4π R^2  [g'_n,2(R)]^2 + 4π L^2 [g'_n,1(R) + μ_n,-^(p)]^2) .
  
  
  
  In the limit p→ 0, one gets
  
  
      g'_n,1(R)    → - n+1 + n(R/L)^2n+1/R(1 - (R/L)^2n+1) ,
      
      g'_n,1(L)    → -  (2n+1) (R/L)^n+1/L(1 - (R/L)^2n+1) ,
      
      g'_n,2(R)    →(2n+1) (R/L)^n/R(1 - (R/L)^2n+1) ,
      
      g'_n,2(L)    →n + (n+1)(R/L)^2n+1/L(1 - (R/L)^2n+1) ,
  
  (and B → n(n+1)/(LR)), from which
  
  
      μ_n,±^(0) = 1/2R(1-γ_n)(n(β+1)+1 + γ_n(n+(n+1)β) 
         ±√((n(β+1)+1 + γ_n(n+(n+1)β))^2 - 4β n(n+1))) ,
  
  where β = R/L and γ_n = (R/L)^2n+1.  We also get
  
  
      μ_0,-^(0) = 0,      μ_0,+^(0) = 1 + β^2/R(1-β)
  
  and
  
  
      a_0,-^(0)    = b_0,-^(0) = 1/√(4π(R^2+L^2)) , 
      
      a_0,+^(0)    = 1/R√(4π (1+β^2)) ,      b_0,+^(0) = -β^2/R √(4π (1 + β^2)) .
  
  In the limit L→∞, one retrieves μ_n,+^(0)→
  (n+1)/R and μ_n,-^(0)→ 0.  In turn, as R→ 0, one has
  μ_n,-^(0)→ n/L, whereas μ_n,+^(0)→∞.
  
  Using the asymptotic behavior of the modified spherical Bessel
  functions, one can also check that
  
  
      lim_R→ 0μ_n,-^(p)    = √(p/D) i'_n(L√(p/D))/i_n(L√(p/D)) , 
      lim_R→ 0μ_n,+^(p)    = ∞ , 
      lim_L→∞μ_n,+^(p)    = - √(p/D) k'_n(R√(p/D))/k_n(R√(p/D)) , 
      lim_L→∞μ_n,-^(p)    = √(p/D) ,
  
  where we used that k'_n(z)/k_n(z) ≤ - 1.  As a consequence, in
  the limit R→ 0, one retrieves the eigenvalues of _p for the
  interior of a ball of radius L.  In turn, in the limit L→∞, μ_n,+^(p) approach the eigenvalues for the exterior
  of a ball of radius R, while μ_n,-^(p) accumulate near
  √(p/D).
  
  Finally, the Dirichlet propagator in the Laplace domain is
  
  
  
      G̃_∞(,p|_0) = ∑_n=0^∞α (2n+1)/4π D  P_n((·_0)/||   |_0|)
         g_n,2(r_0) 
         ×[k_n(α L) i_n(α r) - i_n(α L) k_n(α r)]     (r_0 ≤ r),  
          19mm = ∑_n=0^∞α (2n+1)/4π D  P_n((·_0)/||   |_0|)
         g_n,1(r_0) 
         ×[k_n(α R) i_n(α r) - i_n(α R) k_n(α r)]     (r ≤ r_0),
  
  
  where = (r,θ,ϕ) and _0 = (r_0,θ_0,ϕ_0).  The
  probability flux density reads then
  
  
  
      j̃_∞(,p|_0) 
         = ∑_n=0^∞2n+1/4π R^2  
      P_n((·_0)/||   |_0|)   g_n,1(r_0)     (∈Γ_1) , 
         = ∑_n=0^∞2n+1/4π R^2  
      P_n((·_0)/||   |_0|)   g_n,2(r_0)     (∈Γ_2) .
  
  
  As a consequence, one gets
  
  
      V_n^(p)(_0) = (a_n^(p)  g_n,1(r_0) + b_n^(p)  g_n,2(r_0) ) Y_mn(θ_0,ϕ_0)
  
  that gives access to the full propagator P(,ℓ,t|_0).
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  40
  
  Rice				S. Rice,
  					Diffusion-Limited Reactions
  					(Elsevier, Amsterdam, 1985).
  
  Redner			S. Redner, 
  					A Guide to First Passage Processes 
  					(Cambridge: Cambridge University press, 2001).
  
  Schuss			Z. Schuss,
  					Brownian Dynamics at Boundaries and Interfaces in Physics, Chemistry and Biology
  					(Springer, New York, 2013).
  
  Metzler			R. Metzler, G. Oshanin, and S. Redner (Eds.) 
  					First-Passage Phenomena and Their Applications
  					(Singapore: World Scientific, 2014).
  
  Oshanin			G. Oshanin, R. Metzler, K. Lindenberg (Eds.)
  					Chemical Kinetics: Beyond the Textbook
  					(New Jersey: World Scientific, 2019). 
  
  
  Grebenkov07			D. S. Grebenkov, 
  					“NMR Survey of Reflected Brownian Motion”, 
  					Rev. Mod. Phys. 79, 1077-1137 (2007).
  
  Benichou11			O. Bénichou, C. Loverdo, M. Moreau, and R. Voituriez,
  					“Intermittent search strategies”,
  					Rev. Mod. Phys. 83, 81-130 (2011).
  
  Bressloff13			P. C. Bressloff and J. M. Newby,
  					“Stochastic models of intracellular transport”,
  					Rev. Mod. Phys. 85, 135-196 (2013).
  
  Benichou14			O. Bénichou and R. Voituriez,
  					“From first-passage times of random walks in confinement to geometry-controlled kinetics”,
  					Phys. Rep. 539, 225-284 (2014).
  
  
  
  
  
  Gardiner			C. W. Gardiner,
  					Handbook of stochastic methods for physics, chemistry and the natural sciences
  					(Springer: Berlin, 1985).
  
  Risken			H. Risken,
  					The Fokker-Planck equation: methods of solution and applications,
  					3rd Ed. (Berlin: Springer, 1996).
  
  
  
  
  
  
  Collins49			F. C. Collins and G. E. Kimball,
  					“Diffusion-controlled reaction rates”,
  					J. Coll. Sci. 4, 425 (1949).
  
  Lauffenburger			D. A. Lauffenburger and J. Linderman, 
  					Receptors: Models for Binding, Trafficking, and Signaling
  					(Oxford University Press, 1993).
  
  Sano79			H. Sano and M. Tachiya,
  					“Partially diffusion-controlled recombination”,
  					J. Chem. Phys. 71, 1276-1282 (1979).
  
  Sano81			H. Sano and  M. Tachiya,
  					“Theory of diffusion-controlled reactions on spherical surfaces and 
  					its application to reactions on micellar surfaces”,
  					J. Chem. Phys. 75, 2870-2878 (1981).
  
  Shoup82			D. Shoup and A. Szabo, 
  					“Role of diffusion in ligand binding to macromolecules and cell-bound receptors”, 
  					Biophys. J. 40, 33-39 (1982).
  
  
  
  Zwanzig90			R. Zwanzig, 
  					“Diffusion-controlled ligand binding to spheres partially covered by receptors: an effective medium treatment”,
  					Proc. Natl. Acad. Sci. USA 87, 5856 (1990).
  
  
  
  Sapoval94			B. Sapoval,
  					“General Formulation of Laplacian Transfer Across Irregular Surfaces”,
  					Phys. Rev. Lett. 73, 3314-3317 (1994).
  
  Filoche99			M. Filoche and B. Sapoval, 
  					“Can One Hear the Shape of an Electrode? II. Theoretical Study of the Laplacian Transfer”,
  					Eur. Phys. J. B 9, 755-763 (1999).
  
  Sapoval02			B. Sapoval, M. Filoche, and E. Weibel,
  					“Smaller is better – but not too small: A physical scale for the design of the mammalian pulmonary acinus”,
  					Proc. Nat. Ac. Sci. USA 99, 10411-10416 (2002).
  
  Grebenkov03			D. S. Grebenkov, M. Filoche, and B. Sapoval, 
  					“Spectral Properties of the Brownian Self-Transport Operator”, 
  					Eur. Phys. J. B 36, 221-231 (2003).
  
  Berezhkovskii04		A. Berezhkovskii, Y. Makhnovskii, M. Monine, V. Zitserman, and S. Shvartsman, 
  					“Boundary homogenization for trapping by patchy surfaces”,
  					J. Chem. Phys. 121, 11390 (2004).
  
  Grebenkov05			D. S. Grebenkov, M. Filoche, B. Sapoval, and M. Felici,
  					“Diffusion-Reaction in Branched Structures: Theory and Application to the Lung Acinus”,
  					Phys. Rev. Lett. 94, 050602 (2005).
  
  Grebenkov06a			D. S. Grebenkov, M. Filoche, and B. Sapoval, 
  					“Mathematical Basis for a General Theory of Laplacian Transport towards Irregular Interfaces”, 
  					Phys. Rev. E 73, 021103 (2006).
  
  Traytak07			S. D. Traytak and W. Price,
  					“Exact solution for anisotropic diffusion-controlled reactions with partially reflecting conditions”,
  					J. Chem. Phys. 127, 184508 (2007).
  
  Bressloff08			P. C. Bressloff, B. A. Earnshaw, and M. J. Ward,
  					“Diffusion of protein receptors on a cylindrical dendritic membrane with partially absorbing traps”,
  					SIAM J. Appl. Math. 68, 1223-1246 (2008).
  
  Lawley15			S. D. Lawley and J. P. Keener,
  					“A New Derivation of Robin Boundary Conditions through Homogenization of a Stochastically Switching Boundary”,
  					SIAM J. Appl. Dyn. Sys. 14, 1845-1867 (2015).
  
  Galanti16			M. Galanti, D. Fanelli, S. D. Traytak, and F. Piazza, 
  					“Theory of diffusion-influenced reactions in complex geometries”,
  					Phys. Chem. Chem. Phys. 18, 15950-15954 (2016).
  
  Lindsay17			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, 74-109 (2017).
  
  Grebenkov17			D. S. Grebenkov and G. Oshanin, 
  					“Diffusive escape through a narrow opening: new insights into a classic problem,”
  					Phys. Chem. Chem. Phys. 19, 2723-2739 (2017).
  
  Bernoff18b			A. Bernoff, A. Lindsay, and D. Schmidt, 
  					“Boundary Homogenization and Capture Time Distributions of Semipermeable Membranes 
  					with Periodic Patterns of Reactive Sites”,
  					Multiscale Model. Simul. 16, 1411-1447 (2018).
  
  
  Grebenkov19d			D. S. Grebenkov and S. Traytak, 
  					“Semi-analytical computation of Laplacian Green functions in three-dimensional
  					domains with disconnected spherical boundaries,” 
  					J. Comput. Phys. 379, 91-117 (2019).
  
  
  
  
  
  
  
  
  
  
  
  Grebenkov19b			D. S. Grebenkov,
  					“Imperfect Diffusion-Controlled Reactions”, 
  					in Chemical Kinetics: Beyond the Textbook, Eds. K. Lindenberg, R. Metzler, and G. Oshanin 
  					(World Scientific, New Jersey, 2019).
  
  
  
  
  
  
  Grebenkov20			D. S. Grebenkov,
  					“Paradigm shift in diffusion-mediated surface phenomena”,
  					Phys. Rev. Lett. 125, 078102 (2020).
  
  
  
  Levy				P. Lévy,
  					Processus Stochastiques et Mouvement Brownien
  					(Paris, Gauthier-Villard, 1965).
  
  Ito				K. Ito and H. P. McKean,
  					Diffusion Processes and Their Sample Paths
  					(Springer-Verlag, Berlin, 1965).
  
  Freidlin			M. Freidlin,
  					Functional Integration and Partial Differential Equations
  					(Annals of Mathematics Studies, Princeton University Press, Princeton, New Jersey, 1985).
  
  
  
  Yuste13			S. B. Yuste, E. Abad, and K. Lindenberg, 
  					“Exploration and trapping of mortal random walkers”, 
  					Phys. Rev. Lett. 110, 220603 (2013).
  
  Meerson15			B. Meerson and S. Redner, 
  					“Mortality, redundancy, and diversity in stochastic search”, 
  					Phys. Rev. Lett. 114, 198101 (2015).
  
  Grebenkov17d			D. S. Grebenkov and J.-F. Rupprecht, 
  					“The escape problem for mortal walkers”, 
  					J. Chem. Phys. 146, 084106 (2017).
  
  
  
  
  Bychuk94			O. V. Bychuk and B. O'Shaugnessy, 
  					“Anomalous surface diffusion: A numerical study”,
  					J. Chem. Phys. 101, 772 (1994).
  
  Bychuk95			O. V. Bychuk and B. O'Shaughnessy,
  					“Anomalous Diffusion at Liquid Surfaces”,
  					Phys. Rev. Lett. 74, 1795 (1995).
  
  Walder11			R. Walder, N. Nelson, and D. K. Schwartz,
  					“Single Molecule Observations of Desorption-Mediated Diffusion at the Solid-Liquid Interface”,
  					Phys. Rev. Lett. 107, 156102 (2011).
  
  Skaug13			M. J. Skaug, J. Mabry, and D. K. Schwartz,
  					“Intermittent Molecular Hopping at the Solid-Liquid Interface”,
  					Phys. Rev. Lett. 110, 256101 (2013).
  
  Wang17			D. Wang, H. Wu, and D. K. Schwartz,
  					“Three-Dimensional Tracking of Interfacial Hopping Diffusion”,
  					Phys. Rev. Lett. 119, 268001 (2017).
  
  
  Chechkin09			A. V. Chechkin, I. M. Zaid, M. A. Lomholt, I. M. Sokolov, and R. Metzler,
  					“Bulk-mediated surface diffusion along a cylinder: Propagators and crossovers”,
  					Phys. Rev. E 79, 040105(R) (2009).
  
  Chechkin11			A. V. Chechkin, I. M. Zaid, M. A. Lomholt, I. M. Sokolov, and R. Metzler,
  					“Effective surface motion on a reactive cylinder of particles that perform intermittent bulk diffusion”,
  					J. Chem. Phys. 134, 204116 (2011).
  
  Chechkin12			A. V. Chechkin, I. M. Zaid, M. A. Lomholt, I. M. Sokolov, and R. Metzler,
  					“Bulk-mediated diffusion on a planar surface: Full solution”,
  					Phys. Rev. E 86, 041101 (2012).
  
  Berezhkovskii15		A. M. Berezhkovskii, L. Dagdug, and S. M. Bezrukov,
  					“A new approach to the problem of bulk-mediated surface diffusion”,
  					J. Chem. Phys. 143, 084103 (2015).
  
  Berezhkovskii17		A. M. Berezhkovskii, L. Dagdug, and S. M. Bezrukov,
  					“Bulk-mediated surface transport in the presence of bias”,
  					J. Chem. Phys. 147, 014103 (2017).
  
  
  
  Morters			P. Mörters and Y. Peres,
  					Brownian Motion
  					(Cambridge Series in Statistical and Probabilistic Mathematics,
  					Cambridge University Press, Cambridge, 2010).
  
  
  Grebenkov07a			D. S. Grebenkov,
  					“Residence times and other functionals of reflected Brownian motion”,
  					Phys. Rev. E 76, 041139 (2007). 
  
  Grebenkov19c			D. S. Grebenkov,
  					“Probability distribution of the boundary local time of reflected Brownian motion in Euclidean domains”,
  					Phys. Rev. E 100, 062110 (2019).
  
  
  
  Borodin			A. N. Borodin and P. Salminen,
  					Handbook of Brownian Motion: Facts and Formulae
  					(Birkhauser Verlag, Basel-Boston-Berlin, 1996).
  
  
  Takacs95			L. Takacs, 
  					“On the local time of the Brownian motion”, 
  					Ann. Appl. Probab. 5, 741 (1995).
  
  Randon18			J. Randon-Furling and S. Redner, 
  					“Residence time near an absorbing set”, 
  					J. Stat. Mech. 103205 (2018).
  
  
  
  
  
  Arendt14			W. Arendt, A. F. M. ter Elst, J. B. Kennedy, and M. Sauter,
  					“The Dirichlet-to-Neumann operator via hidden compactness”,
  					J. Funct. Anal. 266, 1757-1786 (2014).
  
  Daners14			D. Daners,
  					“Non-positivity of the semigroup generated by the Dirichlet-to-Neumann operator”,
  					Positivity 18, 235-256 (2014).
  
  Arendt15			W. Arendt and A. F. M. ter Elst,
  					“The Dirichlet-to-Neumann Operator on Exterior Domains”,
  					Potential Anal. 43, 313-340 (2015).
  
  Hassell17			A. Hassell and V. Ivrii,
  					“Spectral asymptotics for the semiclassical Dirichlet to Neumann operator”,
  					J. Spectr. Theory 7, 881-905 (2017).
  
  Girouard17			A. Girouard and I. Polterovich,
  					“Spectral geometry of the Steklov problem”,
  					J. Spectr. Theory 7, 321-359 (2017).
  
  
  
  
  
  
  Grebenkov06			D. S. Grebenkov, 
  					Partially Reflected Brownian Motion: A Stochastic Approach to Transport Phenomena, 
  					in “Focus on Probability Theory”, Ed. L. R. Velle, pp. 135-169 (Nova Science Publishers, 2006). 
  
  Grebenkov06b			D. S. Grebenkov, 
  					“Scaling Properties of the Spread Harmonic Measures”, 
  					Fractals 14, 231-243 (2006).
  
  Grebenkov15			D. S. Grebenkov,
  					“Analytical representations of the spread harmonic measure density”,
  					Phys. Rev. E 91, 052108 (2015).
  
  
  
  Grebenkov20c			D. S. Grebenkov, 
  					“Diffusion toward non-overlapping partially reactive spherical traps: 
  					fresh insights onto classic problems,” 
  					J. Chem. Phys. 152, 244108 (2020).
  
  
  
  
  
  Zaburdaev15			V. Zaburdaev, S. Denisov, and J. Klafter,
  					“Levy walks”,
  					Rev. Mod. Phys. 87, 483 (2015).
  
  
  
  
  Grebenkov20b			D. S. Grebenkov, 
  					“Joint distribution of multiple boundary local times and related first-passage 
  					time problems with multiple targets”
  					(accepted to J. Stat. Mech.)
  
  
  Grebenkov19			D. S. Grebenkov,
  					“Spectral theory of imperfect diffusion-controlled reactions on heterogeneous catalytic surfaces”,
  					J. Chem. Phys. 151, 104108 (2019).
  
  Grebenkov19g			D. S. Grebenkov,
  					“A physicist's guide to explicit summation formulas involving zeros of Bessel functions 
  					and related spectral sums”
  					(submitted; available online arXiv:1904.11190v2)
  
  Grebenkov18			D. S. Grebenkov, R. Metzler, and G. Oshanin, 
  					“Strong defocusing of molecular reaction times results from an interplay 
  					of geometry and reaction control”, 
  					Commun. Chem. 1, 96 (2018).
  
  
  
  
  
  
  
  Torney83			D. C. Torney and H. M. McConnel, 
  					“Diffusion-limited reaction rate theory for two-dimensional systems”, 
  					Proc. R. Soc. London A 387, 147-170 (1983).
  
  Grebenkov19j			D. S. Grebenkov, 
  					“Reversible reactions controlled by surface diffusion on a sphere”, 
  					J. Chem. Phys. 151, 154103 (2019). 
  
  Smoluchowski17		M. Smoluchowski, 
  					“Versuch einer Mathematischen Theorie der Koagulations Kinetic Kolloider Lösungen”,
  					Z. Phys. Chem. 92U, 129-168 (1917).
  
  
  
  
  
  
  
  
  
  