Denis S Grebenkov

Spectral theory of imperfect diffusion-controlled reactions on heterogeneous catalytic surfaces

Keywords: Propagator, Mixed boundary condition, Partial reactivity, Heterogeneity, Reaction rate, Dirichlet-to-Neumann operator

We propose a general theoretical description of chemical reactions occurring on a catalytic surface with heterogeneous reactivity. The propagator of a diffusion-reaction process with eventual absorption on the heterogeneous partially reactive surface is expressed in terms of a much simpler propagator toward a homogeneous perfectly reactive surface. In other words, the original problem with general Robin boundary condition that includes in particular mixed Robin-Neumann condition, is reduced to that with Dirichlet boundary condition. Chemical kinetics on the surface is incorporated as a matrix representation of the surface reactivity in the eigenbasis of the Dirichlet-to-Neumann operator. New spectral representations of important characteristics of diffusion-controlled reactions, such as the survival probability, the distribution of reaction times, and the reaction rate, are deduced. Theoretical and numerical advantages of this spectral approach are illustrated by solving interior and exterior problems for a spherical surface that may describe either an escape from a ball or hitting its surface from outside. The effect of continuously varying or piecewise constant surface reactivity (describing, e.g., many reactive patches) is analyzed.

I. INTRODUCTION

Marian von Smoluchowski first emphasized the importance of diffusive dynamics of reactant molecules and thus laid the foundations for the modern theory of diffusion-controlled reactions [START_REF] Smoluchowski | Versuch einer Mathematischen Theorie der Koagulations Kinetic Kolloider Lösungen[END_REF]. In the basic description, the concentration c(x, t) of molecules, diffusing toward a static catalytic surface with the diffusivity D, obeys the diffusion equation in a bulk domain Ω ∂ ∂t c(x, t) = D∆c(x, t)

(x ∈ Ω) (1)
(with ∆ being the Laplace operator), subject to the Dirichlet boundary condition on the surface ∂Ω: c(x, t) = 0 (x ∈ ∂Ω).

This condition describes a perfect sink, i.e., any molecule hitting the surface reacts with infinite reaction rate upon the first encounter. Since that seminal paper by Smoluchowski, diffusion-controlled reactions to perfect sinks and the related first-passage phenomena have been thoroughly investigated [START_REF] Rice | Diffusion-Limited Reactions[END_REF][START_REF] Redner | A Guide to First Passage Processes[END_REF][START_REF] Levitz | Brownian flights over a fractal nest and first passage statistics on irregular surfaces[END_REF][START_REF] Condamin | First-passage time in complex scaleinvariant media[END_REF][START_REF] Bénichou | Geometry-controlled kinetics[END_REF][START_REF] Bénichou | From first-passage times of random walks in confinement to geometry-controlled kinetics[END_REF][START_REF]First-Passage Phenomena and Their Applications[END_REF].

The assumption of infinite reaction rate is not realistic for most chemical reactions because a molecule that approached a catalytic surface, needs to overcome an activation energy barrier to react that results in a finite reaction rate [START_REF] Weiss | Overview of theoretical models for reaction rates[END_REF][START_REF] Hänggi | Reaction-rate theory: fifty years after Kramers[END_REF]. This effect was first incorporated by Collins and Kimball [START_REF] Collins | Diffusion-controlled reaction rates[END_REF] who replaced Dirichlet boundary condition [START_REF] Rice | Diffusion-Limited Reactions[END_REF] by Robin boundary condition (also known * Electronic address: denis.grebenkov@polytechnique.edu as Fourier, radiation, or third boundary condition):

-D ∂ ∂n c(x, t) = κ(x) c(x, t) (x ∈ ∂Ω), (3)
where ∂/∂n is the normal derivative oriented outward the bulk. This condition states that at each boundary point s ∈ ∂Ω, the diffusive flux density, j = (-D∇c) • n s , in the unit direction n s orthogonal to the surface, is proportional to the concentration at this point. The proportionality coefficient, κ(s), is called the reactivity (with units m/s) and can in general depend on the point s. In this formulation, Robin boundary condition is essentially a mass conservation law at each point of the boundary: the net influx of molecules diffusing toward the boundary is equal to the amount of reacted molecules. The limit κ(s) = 0 (for all s) describes an inert surface without any reaction (i.e., the net diffusive flux at the surface is zero), whereas the limit κ(s) = κ → ∞ reduces Eq. (3) to Eq.

(2) and describes an immediate reaction upon the first encounter. The reactivity is thus related to the probability of reaction event at the encounter [START_REF] Filoche | Can One Hear the Shape of an Electrode? II. Theoretical Study of the Laplacian Transfer[END_REF][START_REF] Grebenkov | Spectral Properties of the Brownian Self-Transport Operator[END_REF][START_REF] Grebenkov | Partially Reflected Brownian Motion: A Stochastic Approach to Transport Phenomena[END_REF]. Robin boundary condition with homogeneous (constant) reactivity κ was often employed to describe many chemical and biochemical reactions and permeation processes [START_REF] Lauffenburger | Receptors: Models for Binding, Trafficking, and Signaling[END_REF][START_REF] Sano | Partially diffusion-controlled recombination[END_REF][START_REF] Sano | Theory of diffusion-controlled reactions on spherical surfaces and its application to re-actions on micellar surfaces[END_REF][START_REF] Shoup | Role of diffusion in ligand binding to macromolecules and cell-bound receptors[END_REF][START_REF]General Formulation of Laplacian Transfer Across Irregular Surfaces[END_REF][START_REF] Sapoval | Smaller is better -but not too small: A physical scale for the design of the mammalian pulmonary acinus[END_REF][START_REF] Grebenkov | Diffusion-Reaction in Branched Structures: Theory and Application to the Lung Acinus[END_REF][START_REF] Qian | Time dependent diffusion in a disordered medium with partially absorbing walls: A perturbative approach[END_REF][START_REF] Traytak | Exact solution for anisotropic diffusion-controlled reactions with partially reflecting conditions[END_REF][START_REF] Bressloff | Diffusion of protein receptors on a cylindrical dendritic membrane with partially absorbing traps[END_REF][START_REF] Galanti | Theory of diffusion-influenced reactions in complex geometries[END_REF][START_REF] Grebenkov | Semi-analytical computation of Laplacian Green functions in threedimensional domains with disconnected spherical boundaries[END_REF], to model stochastic gating [START_REF] Bénichou | Kinetics of stochastically gated diffusion-limited reactions and geometry of random walk trajectories[END_REF][START_REF] Reingruber | Gated Narrow Escape Time for Molecular Signaling[END_REF][START_REF] Lawley | A New Derivation of Robin Boundary Conditions through Homogenization of a Stochastically Switching Boundary[END_REF][START_REF] Bressloff | Stochastic switching in biology: from genotype to phenotype[END_REF], or to approximate the effect of microscopic heterogeneities in a random distribution of reactive sites [START_REF] Berg | Physics of chemoreception[END_REF][START_REF] Shoup | Diffusioncontrolled bimolecular reaction rates. The effect of rotational diffusion and orientation constraints[END_REF][START_REF] Zwanzig | Time dependent rate of diffusion-influenced ligand binding to receptors on cell surfaces[END_REF] (see a recent overview in [START_REF] Grebenkov | Imperfect Diffusion-Controlled Reactions[END_REF]).

In spite of its practical importance, diffusion-controlled reactions with heterogeneous surface reactivity κ(s) remain much less studied. In fact, when κ(s) is not constant, the eigenfunctions of the Laplace operator with Robin boundary condition are not known explicitly even for simple domains (e.g. a ball) that prohibits using standard spectral decompositions [START_REF] Carslaw | Conduction of Heat in Solids[END_REF][START_REF] Crank | The Mathematics of Diffusion[END_REF], on which most of classical solutions are based. One needs therefore to resort to numerical tools such as finite element or finite difference methods for solving the diffusion equation, or to Monte Carlo simulations. A notable exception is the case of piecewise constant reactivity that describes a target Γ (or multiple targets Γ i ) with a constant reactivity κ on the otherwise inert surface. This situation corresponds to Robin-Neumann (for 0 < κ < ∞) or Dirichlet-Neumann (for κ = ∞) mixed boundary conditions [START_REF] Sneddon | Mixed Boundary Value Problems in Potential Theory[END_REF][START_REF] Duffy | Mixed boundary value problems[END_REF]. The Dirichlet-Neumann boundary value problem has been particularly well studied for the Poisson and Laplace equations determining the mean first-passage time to a small target and the reaction rate, respectively (see overviews in [START_REF] Holcman | The Narrow Escape Problem[END_REF][START_REF] Schuss | Brownian Dynamics at Boundaries and Interfaces in Physics, Chemistry and Biology[END_REF][START_REF] Holcman | Stochastic Narrow Escape in Molecular and Cellular Biology[END_REF] and references therein). On one hand, matched asymptotic analysis, dual series technique, and conformal mapping were applied to establish the behavior of the mean first-passage time in both twoand three-dimensional domains [START_REF] Singer | Narrow Escape, Part I[END_REF][START_REF] Singer | Narrow Escape, Part II The circular disk[END_REF][START_REF] Singer | Narrow Escape, Part III Riemann surfaces and non-smooth domains[END_REF][START_REF] Pillay | An Asymptotic Analysis of the Mean First Passage Time for Narrow Escape Problems: Part I: Two-Dimensional Domains[END_REF][START_REF] Cheviakov | An Asymptotic Analysis of the Mean First Passage Time for Narrow Escape Problems: Part II: The Sphere[END_REF][START_REF] Cheviakov | Mathematical modeling and numerical computation of narrow escape problems[END_REF][START_REF] Caginalp | Analytical and Numerical Results for an Escape Problem[END_REF][START_REF] Marshall | Analytical Solutions for an Escape Problem in a Disc with an Arbitrary Distribution of Exit Holes Along Its Boundary[END_REF][START_REF] Grebenkov | Universal formula for the mean first passage time in planar domains[END_REF]. On the other hand, homogenization techniques were used to substitute piecewise constant reactivity κ(s) by an effective homogeneous reactivity [START_REF] Berg | Physics of chemoreception[END_REF][START_REF] Shoup | Diffusioncontrolled bimolecular reaction rates. The effect of rotational diffusion and orientation constraints[END_REF][START_REF] Zwanzig | Time dependent rate of diffusion-influenced ligand binding to receptors on cell surfaces[END_REF][START_REF] Zwanzig | Diffusion-controlled ligand binding to spheres partially covered by receptors: an effective medium treatment[END_REF][START_REF] Berezhkovskii | Boundary homogenization for trapping by patchy surfaces[END_REF][START_REF] Berezhkovskii | Homogenization of boundary conditions for surfaces with regular arrays of traps[END_REF][START_REF] Muratov | Boundary homogenization for periodic arrays of absorbers[END_REF][START_REF] Bernoff | Boundary Homogenization and Capture Time Distributions of Semipermeable Membranes with Periodic Patterns of Reactive Sites[END_REF][START_REF] Dagdug | Boundary homogenization for a sphere with an absorbing cap of arbitrary size[END_REF][START_REF] Lindsay | First Passage Statistics for the Capture of a Brownian Particle by a Structured Spherical Target with Multiple Surface Traps[END_REF][START_REF] Bernoff | Numerical approximation of diffusive capture rates by planar and spherical surfaces with absorbing pores[END_REF]. More recent works investigated how the mean reaction time is affected by a finite lifetime of diffusing particles [START_REF] Yuste | Exploration and trapping of mortal random walkers[END_REF][START_REF] Meerson | Mortality, redundancy, and diversity in stochastic search[END_REF][START_REF] Grebenkov | The escape problem for mortal walkers[END_REF], by partial reactivity and interactions [START_REF] Grebenkov | Diffusive escape through a narrow opening: new insights into a classic problem[END_REF][START_REF] Agranov | Narrow Escape of Interacting Diffusing Particles[END_REF], by target aspect ratio [START_REF] Grebenkov | Effects of the target aspect ratio and intrinsic reactivity onto diffusive search in bounded domains[END_REF], by reversible target-binding kinetics [START_REF] Grebenkov | First passage times for multiple particles with reversible target-binding kinetics[END_REF][START_REF] Lawley | First passage time distribution of multiple impatient particles with reversible binding[END_REF] and surface-mediated diffusion [START_REF] Bénichou | Optimal Reaction Time for Surface-Mediated Diffusion[END_REF][START_REF] Bénichou | Mean First-Passage Time of Surface-Mediated Diffusion in Spherical Domains[END_REF][START_REF] Rupprecht | Kinetics of Active Surface-Mediated Diffusion in Spherically Symmetric Domains[END_REF][START_REF] Rupprecht | Exact mean exit time for surface-mediated diffusion[END_REF], by heterogeneous diffusivity [START_REF] Vaccario | First-Passage Times in d-Dimensional Heterogeneous Media[END_REF], and by rapid re-arrangments of the medium [START_REF] Jain | Diffusing Diffusivity: Survival in a Crowded Rearranging and Bounded Domain[END_REF][START_REF] Lanoiselée | Diffusion-limited reactions in dynamic heterogeneous media[END_REF][START_REF] Sposini | First passage statistics for diffusing diffusivity[END_REF]. Some of the related effects onto the whole distribution of reaction times were analyzed [START_REF] Godec | Universal Proximity Effect in Target Search Kinetics in the Few-Encounter Limit[END_REF][START_REF] Godec | First passage time distribution in heterogeneity controlled kinetics: going beyond the mean first passage time[END_REF][START_REF] Grebenkov | Towards a full quantitative description of single-molecule reaction kinetics in biological cells[END_REF][START_REF] Grebenkov | Strong defocusing of molecular reaction times results from an interplay of geometry and reaction control[END_REF][START_REF] Hartich | Duality between relaxation and first passage in reversible Markov dynamics: rugged energy landscapes disentangled[END_REF]. However, the current understanding of diffusion-controlled reactions on catalytic surfaces with continuously varying heterogeneous reactivity remains episodic.

In this paper, we propose a mathematical description of diffusion-controlled reactions on catalytic surfaces, in which chemical kinetics, characterized by heterogeneous surface reactivity κ(s), is disentangled from the firstpassage diffusive steps. In Sec. II, we express the propagator of the sophisticated diffusion-reaction process with multiple reflections on partially reactive surface in terms of a much simpler Dirichlet propagator toward a homogeneous perfectly reactive surface with Dirichlet boundary condition. Chemical kinetics is incorporated via a matrix representation of the heterogeneous surface reactivity in the eigenbasis of the Dirichlet-to-Neumann operator, which is also tightly related to the Dirichlet propagator. From the propagator, we deduce other important characteristics of diffusion-controlled reactions such as the survival probability, the distribution of reaction times and the reaction rate. This formalism provides a general description of such processes and brings conceptually new tools for its investigation. In Sec. III, this spectral approach is applied to an important example of a spherical surface for which the Dirichlet propagator and the Dirichlet-to-Neumann operator are known explicitly. We study both the interior and exterior problems that may describe either an escape from a ball or hitting its surface from outside. Semi-analytical solutions for the probability density of reaction times and for the reaction rate are derived. In Sec. IV, we discuss the advantages and limitations of the spectral approach, its possible extensions, and further applications, in particular, for analytical and numerical studies of mixed boundary value problems. Technical derivations are reported in Appendices.

II. GENERAL SPECTRAL DESCRIPTION

We consider a molecule diffusing with the diffusion coefficient D in an Euclidean domain Ω ⊂ R d toward a partially reactive catalytic boundary ∂Ω characterized by a prescribed heterogeneous (space-dependent) reactivity 0 ≤ κ(s) < ∞. Once the molecule hits the boundary at some point s, it may either react or be reflected back to resume its diffusion until the next encounter, and so on. The reaction probability at each encounter is characterized by the reactivity κ(s) at the encounter point. In this way, the molecule performs multiple diffusive excursions in the bulk until reaction occurs. The finite reactivity results therefore in a very sophisticated diffusive dynamics near the catalytic surface, which is much more intricate than just the first arrival to a homogeneous perfectly reactive surface. A probabilistic construction of this diffusive process (called partially reflected Brownian motion) was discussed in [START_REF] Grebenkov | Spectral Properties of the Brownian Self-Transport Operator[END_REF][START_REF] Grebenkov | Partially Reflected Brownian Motion: A Stochastic Approach to Transport Phenomena[END_REF][START_REF] Papanicolaou | The probabilistic solution of the third boundary value problem for second order elliptic equations[END_REF][START_REF] Bass | On the Robin problem in Fractal Domains[END_REF][START_REF] Grebenkov | Residence times and other functionals of reflected Brownian motion[END_REF][START_REF] Singer | Partially reflected diffusion[END_REF][START_REF] Grebenkov | Laplacian Eigenfunctions in NMR. II Theoretical Advances[END_REF] (see an overview in [START_REF] Grebenkov | Imperfect Diffusion-Controlled Reactions[END_REF]).

Without dwelling on the probabilistic aspects of the problem, we aim at characterizing such diffusion-reaction processes via the propagator G(x, t|x 0 ) (also known as heat kernel or Green's function). This is the probability density for a molecule that has not reacted until time t on the partially reactive boundary ∂Ω, to be in a vicinity of a point x at time t, given that it was started at a point x 0 at time 0. For any fixed starting point x 0 ∈ Ω = Ω ∪ ∂Ω, the propagator satisfies the following boundary value problem

∂G(x, t|x 0 ) ∂t -D∆G(x, t|x 0 ) = 0 (x ∈ Ω), (4a) G(x, t = 0|x 0 ) = δ(x -x 0 ), (4b)
D ∂ ∂n x + κ(x) G(x, t|x 0 ) = 0 (x ∈ ∂Ω), (4c)
where δ(x-x 0 ) is the Dirac distribution, and the Laplace operator ∆ acts on x. If the domain Ω is unbounded, these equations are completed by the regularity condition at infinity: G(x, t|x 0 ) → 0 as |x| → ∞. To avoid technicalities, we assume that the boundary ∂Ω is smooth. The following discussion extends our former results [START_REF] Grebenkov | Partially Reflected Brownian Motion: A Stochastic Approach to Transport Phenomena[END_REF][START_REF] Grebenkov | Imperfect Diffusion-Controlled Reactions[END_REF][START_REF] Grebenkov | Residence times and other functionals of reflected Brownian motion[END_REF][START_REF] Grebenkov | Mathematical Basis for a General Theory of Laplacian Transport towards Irregular Interfaces[END_REF] to heterogeneous reactivity and time-dependent diffusion equation. We consider the Laplace-transformed propagator,

G(x, p|x 0 ) = ∞ 0 dt e -pt G(x, t|x 0 ), (5)
which satisfies the modified Helmholtz equation for each fixed x 0 ∈ Ω:

(p -D∆) G(x, p|x 0 ) = δ(x -x 0 ), (6a)
D ∂ ∂n x + κ(x) G(x, p|x 0 ) = 0 (x ∈ ∂Ω) (6b)
(tilde will denote Laplace-transformed quantities).

Our goal is to express the propagator G(x, p|x 0 ) describing diffusion toward heterogeneous partially reactive surface ∂Ω in terms of the much simpler Dirichlet propagator G0 (x, p|x 0 ) that characterizes diffusion toward the homogeneous perfectly reactive surface and satisfies for each fixed x 0 ∈ Ω:

(p -D∆) G0 (x, p|x 0 ) = δ(x -x 0 ), (7a) G0 (x, p|x 0 ) = 0 (x ∈ ∂Ω). ( 7b
)
Due to the linearity of the problem (6), one can search its solution in the form G(x, p|x 0 ) = G0 (x, p|x 0 ) + g(x, p|x 0 ), [START_REF]First-Passage Phenomena and Their Applications[END_REF] where the unknown regular part g(x, p|x 0 ) satisfies

(p -D∆)g(x, p|x 0 ) = 0 (x ∈ Ω), (9a)
D ∂ ∂n x + κ(x) g(x, p|x 0 ) = j0 (x, p|x 0 ) (9b) for x ∈ ∂Ω, where j0 (s, p|x 0 ) = -D ∂ ∂n x G0 (x, p|x 0 ) x=s (s ∈ ∂Ω) (10
) is the Laplace transform of the diffusive flux density j 0 (s, t|x 0 ) at time t in a point s of the homogeneous perfectly reactive surface (i.e., the probability density of the first arrival in a vicinity of s at time t after starting from x 0 at time 0).

A. Dirichlet-to-Neumann operator

The solution of the boundary value problem ( 9) can be obtained with the help of the Dirichlet-to-Neumann operator M p (also known as Poincaré-Steklov operator) [START_REF] Yu | Pseudo-differential Operators, Singularities, Applications[END_REF][START_REF] Jacob | Pseudo-differential Operators and Markov Processes[END_REF][START_REF] Taylor | Pseudodifferential Operators[END_REF]. This is a pseudo-differential self-adjoint operator that associates to a function f on the boundary ∂Ω another function on that boundary:

[M p f ](s) = ∂ ũ(x, p) ∂n x=s (s ∈ ∂Ω), (11)
where ũ(x, p) is the solution of the Dirichlet boundary value problem:

(p -D∆)ũ(x, p) = 0 (x ∈ Ω), (12a) ũ(x, p) = f (x, p) (x ∈ ∂Ω) (12b)
(here we skip the usual regularity assumptions on ∂Ω, as well as the explicit description of the functional spaces involved in the rigorous definition of M p , see [START_REF] Yu | Pseudo-differential Operators, Singularities, Applications[END_REF][START_REF] Jacob | Pseudo-differential Operators and Markov Processes[END_REF][START_REF] Taylor | Pseudodifferential Operators[END_REF][START_REF] Marletta | Eigenvalue problems on exterior domains and Dirichlet to Neumann maps[END_REF][START_REF] Arendt | Spectral properties of the Dirichlet-to-Neumann operator on Lipschitz domains[END_REF][START_REF] Arendt | The Dirichlet-to-Neumann Operator on Exterior Domains[END_REF][START_REF] Hassell | Spectral asymptotics for the semiclassical Dirichlet to Neumann operator[END_REF] for details). For instance, if f is understood as a source of molecules on the boundary ∂Ω emitted into the reactive bulk, then the operator M p gives their flux density on that boundary. Note that there is a family of operators parameterized by p (or p/D).

As the solution of the Dirichlet boundary value problem [START_REF] Filoche | Can One Hear the Shape of an Electrode? II. Theoretical Study of the Laplacian Transfer[END_REF] can be expressed in terms of the Dirichlet propagator G0 (x, p|x 0 ) in a standard way,

ũ(x, p) = ∂Ω ds ′ j0 (s ′ , p|x) f (s ′ , p),
the Dirichlet-to-Neumann propagator acts formally as

[M p f ](s) =   ∂ ∂n ∂Ω ds ′ j0 (s ′ , p|x) f (s ′ , p)   x=s , (13)
and thus the Dirichlet propagator determines the Dirichlet-to-Neumann operator M p . In Appendix A, it is also shown how the Dirichlet propagator can be constructed from the operator M p . As a consequence, these two important objects are equivalent. As discussed in [START_REF] Grebenkov | Partially Reflected Brownian Motion: A Stochastic Approach to Transport Phenomena[END_REF] for the case p = 0, the Dirichlet-to-Neumann operator can also be interpreted as the continuous limit of the Brownian self-transport operator Q ij which was introduced in [START_REF] Filoche | Can One Hear the Shape of an Electrode? II. Theoretical Study of the Laplacian Transfer[END_REF][START_REF] Grebenkov | Spectral Properties of the Brownian Self-Transport Operator[END_REF] to describe the probability of the first arrival to a site j of a discretized boundary from another site i via bulk diffusion.

Let us now return to the boundary value problem [START_REF] Weiss | Overview of theoretical models for reaction rates[END_REF]. Suppose that we have solved this problem and found that the solution g(x, p|x 0 ) on the boundary ∂Ω is equal to some function f (x, p). Applying then the Dirichlet-to-Neumann operator to f (x, p), one can express the normal derivative of g(x, p|x 0 ), from which

g(s, p|x 0 ) = M p + K -1 j0 (s, p|x 0 ) D (s ∈ ∂Ω), ( 14
)
where K is the operator of multiplication by κ(s)/D. Knowing the restriction of g(x, p|x 0 ) on the boundary ∂Ω, one can reconstruct this function in the bulk Ω as the solution of the corresponding Dirichlet problem:

g(x, p|x 0 ) = ∂Ω ds j0 (s, p|x) g(s, p|x 0 ). ( 15
)
In this way, we obtain the desired representation of the propagator in the form of a scalar product between two functions on the boundary G(x, p|x 0 ) = G0 (x, p|x 0 ) ( 16)

+ 1 D j0 (•, p|x) • (M p + K) -1 j0 (•, p|x 0 ) L2(∂Ω)
, where (f • g) L2(∂Ω) denotes the standard scalar product between functions f and g on the boundary ∂Ω:

(f • g) L2(∂Ω) = ∂Ω ds f (s) g * (s),
and asterisk denotes the complex conjugate. Equation ( 16) is the first main result of the paper. Remarkably, all the "ingredients" of this formula correspond to the Dirichlet condition on a homogeneous perfectly reactive boundary, except for the operator K that keeps track of heterogeneous surface reactivity κ(s). We outline that Eq. ( 16) does not solve the original problem but reduces it to a much simpler and more thoroughly studied Dirichlet problem. When x and x 0 are boundary points, the identity j0 (s, p|s 0 ) = δ(ss 0 ) reduces Eq. ( 16) to

D G(s, p|s 0 ) = M p +K -1 δ(s-s 0 ) (s 0 , s ∈ ∂Ω), (17)
i.e., D G(s, p|s 0 ) is the kernel of the operator M p + K.

One can therefore rewrite Eq. ( 16) as

G(x, p|x 0 ) = G0 (x, p|x 0 ) (18)
G(s, t|x 0 ) = ∂Ω ds 1 t 0 dt 1 j 0 (s 1 , t 1 |x 0 ) G(s, t -t 1 |s 1 )
. [START_REF] Sapoval | Smaller is better -but not too small: A physical scale for the design of the mammalian pulmonary acinus[END_REF] This relation justifies the qualitative separation of the diffusion-reaction process into two steps: the first arrival step (described by j 0 (s 1 , t 1 |x 0 )) and the reaction step (described by G(s, t -t 1 |s 1 )). We stress, however, that the reaction step involves intricate diffusion process near the partially reactive catalytic surface. In addition to the new conceptual view onto partially reflected Brownian motion, the representations [START_REF]General Formulation of Laplacian Transfer Across Irregular Surfaces[END_REF][START_REF] Sapoval | Smaller is better -but not too small: A physical scale for the design of the mammalian pulmonary acinus[END_REF] can be helpful for a numerical computation of the propagator because only the boundary-to-boundary transport via G(s 2 , t|s 1 ) needs to be determined. This kernel significantly extends the Brownian self-transport operator introduced in [START_REF] Filoche | Can One Hear the Shape of an Electrode? II. Theoretical Study of the Laplacian Transfer[END_REF][START_REF] Grebenkov | Spectral Properties of the Brownian Self-Transport Operator[END_REF] (see below).

B. Other common diffusion characteristics

The propagator G(x, t|x 0 ) determines many quantities often considered in the context of diffusion-controlled reactions such as the survival probability up to time t, the reaction time distribution, the distribution of reaction points (at which reaction occurs), and the reaction rate. For instance, the diffusive flux density at a partially reactive point s ∈ ∂Ω is

j(s, t|x 0 ) = -D ∂G(x, t|x 0 ) ∂n x x=s = κ(s)G(s, t|x 0 ), (21)
where we used the Robin boundary condition (4c). This is the joint probability density for the reaction time and the reaction point on the catalytic surface. The integral over s yields the marginal probability density of reaction times, [START_REF] Qian | Time dependent diffusion in a disordered medium with partially absorbing walls: A perturbative approach[END_REF] whereas the integral over t gives the marginal probability density of reaction points:

H(t|x 0 ) = ∂Ω ds j(s, t|x 0 ) = ∂Ω ds κ(s) G(s, t|x 0 ),
ω(s|x 0 ) = ∞ 0 dt j(s, t|x 0 ) = j(s, 0|x 0 ) = κ(s) G(s, 0|x 0 ). (23)
The latter was called the spread harmonic measure density [START_REF] Grebenkov | Partially Reflected Brownian Motion: A Stochastic Approach to Transport Phenomena[END_REF][START_REF] Grebenkov | Mathematical Basis for a General Theory of Laplacian Transport towards Irregular Interfaces[END_REF][START_REF] Grebenkov | Scaling Properties of the Spread Harmonic Measures[END_REF][START_REF] Grebenkov | Analytical representations of the spread harmonic measure[END_REF]. This is a natural extension of the harmonic measure density j0 (s, 0|x 0 ) that characterizes the first arrival onto the perfectly reactive surface [START_REF] Garnett | Harmonic Measure[END_REF][START_REF] Grebenkov | Multifractal Properties of the Harmonic Measure on Koch Boundaries in Two and Three Dimensions[END_REF][START_REF] Grebenkov | What Makes a Boundary Less Accessible[END_REF]. As the probability density H(t|x 0 ) can be interpreted as the probability flux onto the surface for a molecule started from x 0 , its integral with the initial concentration of molecules, c 0 (x 0 ), yields the overall diffusive flux onto the surface, i.e., the reaction rate:

J(t) = Ω dx 0 c 0 (x 0 ) H(t|x 0 ). ( 24
)
In turn, the integral of H(t|x 0 ) from t to infinity gives the survival probability up to time t:

S(t|x 0 ) = 1 - t 0 dt ′ H(t ′ |x 0 ), (25)
while 1 -S(t|x 0 ) is the probability of reaction up to time t. All these quantities are expressed in terms of the propagator and thus determined from Eq. ( 16).

C. Spectral decompositions

When the boundary ∂Ω is bounded, the Dirichletto-Neumann operator M p has a discrete spectrum, with a set of nonnegative eigenvalues µ

M p v (p) n (s) = µ (p) n v (p) n (s) (n = 0, 1, . . .). (26)
We emphasize that both µ

(p) n and v (p)
n depend in general on p as a parameter. Expanding the scalar product in Eq. ( 16) over this basis, one gets

G(x, p|x 0 ) = G0 (x, p|x 0 ) (27) + 1 D ∞ n,n ′ =0 V (p) n (x 0 ) (M + K) -1 n,n ′ [V (p) n ′ (x)] * ,
where

V (p) n (x 0 ) = ∂Ω ds j0 (s, p|x 0 ) v (p) n (s) (28)
is the projection of the Laplace-transformed flux density j0 (s, p|x 0 ) onto the eigenfunction v

n (s), and

M n,n ′ = δ nn ′ µ (p) n , (29a)
K n,n ′ = ∂Ω ds [v (p) n (s)] * κ(s) D v (p) n ′ (s) (29b)
are infinite-dimensional matrices that represent the Dirichlet-to-Neumann operator M p and the reactivity multiplication operator K in the basis of eigenfunctions v

(p) n (s).
From the spectral representation [START_REF] Bénichou | Kinetics of stochastically gated diffusion-limited reactions and geometry of random walk trajectories[END_REF] and Eq. ( 21), we deduce

j(s, p|x 0 ) = ∞ n,n ′ =0 V (p) n (x 0 ) (M + K) -1 K n,n ′ [v (p) n ′ (s)] * , (30)
where we used the completeness of eigenfunctions v (p) n to represent κ(s)/D as multiplication by the matrix K. According to Eqs. [START_REF] Qian | Time dependent diffusion in a disordered medium with partially absorbing walls: A perturbative approach[END_REF][START_REF] Traytak | Exact solution for anisotropic diffusion-controlled reactions with partially reflecting conditions[END_REF], the spectral decomposition (30) yields immediately

ω(s|x 0 ) = ∞ n,n ′ =0 V (0) n (x 0 ) (M + K) -1 K (p=0) n,n ′ [v (0) n ′ (s)] * (31) and
H(p|x 0 ) = |∂Ω| 1/2 ∞ n=0 h (p) n V (p) n (x 0 ), (32)
where

h (p) n = |∂Ω| -1/2 ∞ n ′ =0 (M + K) -1 K n,n ′ ∂Ω ds [v (p) n ′ (s)] * (33)
are dimensionless coefficients. In particular, H(0|x 0 ) is the reaction probability (in Appendix B 1, we prove the expected identity H(0|x 0 ) = 1 for any bounded domain). According to Eq. ( 24), the Laplace-transformed reaction rate is then

J(p) = |∂Ω| 1/2 ∞ n,n ′ =0 h (p) n Ω dx 0 V (p) n (x 0 ) c 0 (x 0 ). ( 34
)
In Appendix B 2, we show how this expression can be further simplified in the case of the uniform initial concentration.

While we mainly focus on Laplace-transformed quantities, their representations in time domain can be obtained via Laplace transform inversion either analytically or numerically. For instance, the inversion in the case of bounded domains can be performed via the residue theorem by computing the poles {p n } ⊂ C of functions in Eqs. [START_REF] Bénichou | Kinetics of stochastically gated diffusion-limited reactions and geometry of random walk trajectories[END_REF][START_REF] Shoup | Diffusioncontrolled bimolecular reaction rates. The effect of rotational diffusion and orientation constraints[END_REF][START_REF] Grebenkov | Imperfect Diffusion-Controlled Reactions[END_REF], which are determined by the condition

det(M + K) = 0. (35)
In general, the spectral representation ( 27) is not simpler than Eq. ( 16) because all V (p) n , M and K depend on p as a parameter. However, in some domains, these "ingredients" can be evaluated explicitly, providing a semianalytical form of the Laplace-transformed propagator and related quantities. We will illustrate this point in Sec. III for a spherical boundary.

D. Homogeneous partial reactivity

In the particular case of homogeneous reactivity, κ(s) = κ, the operator K is proportional to the identity operator, and Eq. ( 17) implies that D G(s, p|s 0 ) is the resolvent of the Dirichlet-to-Neumann operator M p . Moreover, as

K n,n ′ = δ n,n ′ κ D , (36)
Eq. ( 27) is reduced to

Ghom (x, p|x 0 ) = G0 (x, p|x 0 ) + ∞ n=0 V (p) n (x 0 ) [V (p) n (x)] * Dµ (p) n + κ . (37)
In turn, the condition [START_REF] Carslaw | Conduction of Heat in Solids[END_REF] on the poles is reduced to a set of decoupled equations

µ (p) n + κ D = 0, (38)
showing how the eigenvalues µ (p) n of the Dirichlet-to-Neumann operator determine the eigenvalues of the associated Laplace operator with Robin boundary condition.

The other spectral decompositions are also simplified:

ω hom (s|x 0 ) = ∞ n=0 V (0) n (x 0 ) [v (0) n (s)] * D κ µ (0) n + 1 , (39)
Hhom (p|x 0 ) = ∞ n=0 V (p) n (x 0 ) ∂Ω ds [v (p) n (s)] * D κ µ (p) n + 1 . ( 40
)
and

Jhom (p) = c 0 D p ∞ n=0 µ (p) n D κ µ (p) n + 1 ∂Ω ds v (p) n (s) 2 , (41)
where we used Eq. (B8) for the uniform initial concentration c 0 . In the limit p → 0, one recovers the formula for the total steady-state flux derived in Ref. [START_REF] Grebenkov | Mathematical Basis for a General Theory of Laplacian Transport towards Irregular Interfaces[END_REF]; Eq. ( 41) is therefore its extension to time-dependent diffusion. To our knowledge, Eqs. [START_REF] Sneddon | Mixed Boundary Value Problems in Potential Theory[END_REF][START_REF] Holcman | The Narrow Escape Problem[END_REF][START_REF] Schuss | Brownian Dynamics at Boundaries and Interfaces in Physics, Chemistry and Biology[END_REF][START_REF] Holcman | Stochastic Narrow Escape in Molecular and Cellular Biology[END_REF] that are fully explicit in terms of the eigenvalues and eigenfunctions of the Dirichlet-to-Neumann operator, have not been earlier reported. While alternative spectral decompositions on the Laplace operator eigenfunctions are known for bounded domains, there is no such expansion for unbounded domains, for which the spectrum of the Laplace operator is continuous. The spectral formulation in terms of the eigenfunctions of the Dirichletto-Neumann operator opens therefore new perspectives for studying diffusion-reaction processes even for homogeneous reactivity. From the numerical point of view, the computation of the eigenfunctions of the Dirichletto-Neumann operator could in general be simpler due to the reduced dimensionality:

v (p)
n need to be found on the boundary ∂Ω, whereas the Laplace operator eigenfunctions have to be computed in the whole domain Ω.

III. SPHERICAL BOUNDARY

In this section, we apply our general spectral decompositions to the case of a spherical boundary for which the eigenbasis of the Dirichlet-to-Neumann operator is known explicitly. We first discuss in Sec. III A the interior problem that may describe, for instance, an escape from a ball, and then in Sec. III B we dwell on the exterior problem and related chemical kinetics. In both cases, we provide semi-analytical solutions for an arbitrary heterogeneous surface reactivity and then discuss some particular cases, e.g., a piecewise constant reactivity that describes single or multiple reactive targets on the otherwise inert boundary. Technical details of calculations are reported in Appendices C, D, and E.

A. Diffusion inside a ball

We consider a diffusion-reaction process inside a ball of radius R, Ω = {x ∈ R 3 : |x| < R}, with a prescribed heterogeneous surface reactivity κ(s). For this domain, the eigenvalues and eigenfunctions of the Dirichlet-to-Neumann operator are known explicitly (Appendix C),

µ (p) nm = p/D i ′ n (R p/D) i n (R p/D) , (42a)
v nm (θ, φ) = 1 R Y mn (θ, φ), (42b)
where i n (z) are the modified spherical Bessel functions of the first kind, Y mn (θ, φ) are the L 2 (∂Ω)-normalized spherical harmonics, prime denotes the derivative with respect to the argument, and we used spherical coordinates (r, θ, φ). Here we employ the double index nm to enumerate the eigenfunctions as well as the elements of the matrices M and K. Note that the eigenvalues do not depend on the index m and thus are of multiplicity 2n + 1, whereas the eigenfunctions v nm do not depend on the parameter p. The eigenvalues determine the matrix M via Eq. (29a), while Eq. (29b) for the matrix K reads

K nm,n ′ m ′ = π 0 dθ sin θ 2π 0 dφ κ(θ, φ) D Y * mn (θ, φ)Y m ′ n ′ (θ, φ).
(43) The calculation of the matrix K in Eq. ( 43) involves integrals with spherical harmonics that can often be evaluated explicitly. In Appendix D, we discuss several common situations such as a single target, multiple nonoverlapping targets of circular shape or multiple latitudinal stripes, axisymmetric reactivity κ(θ, φ) = κ(θ), and an expansion of κ(θ, φ) into a finite sum over spherical harmonics. Although cumbersome, resulting expressions for the matrix K are exact and do not involve numerical quadrature, providing a powerful computational tool. These cases can further be extended by adding another concentric surface with reflecting or absorbing boundary condition. This modification does not change the matrix K but affects the eigenvalues of the Dirichlet-to-Neumann operator and thus the matrix M.

As the Dirichlet propagator is also known, we deduce in Appendix C

V (p) nm (x 0 ) = R -1 i n (r 0 p/D) i n (R p/D) Y mn (θ 0 , φ 0 ) , (44)
so that the Laplace-transformed propagator G(x, p|x 0 ) is determined in the semi-analytical form [START_REF] Bénichou | Kinetics of stochastically gated diffusion-limited reactions and geometry of random walk trajectories[END_REF], in which the dependence on points x 0 and x is fully explicit, whereas the computation of the coefficients involves a numerical inversion of the matrix M + K. Similarly, one gets semianalytical expressions for the Laplace-transformed probability density of reaction times and the spread harmonic measure (see Appendix C), e.g.,

H(p|x 0 ) = √ 4π ∞ n=0 n m=-n h (p) nm i n (r 0 p/D) i n (R p/D) Y mn (θ 0 , φ 0 ), (45)
with

h (p) nm = (M + K) -1 K nm,00 . (46)
The Laplace-transformed survival probability is related to H(p|x 0 ) as

S(p|x 0 ) = 1 -H(p|x 0 ) p , (47)
whereas the mean reaction time is simply S(0|x 0 ). Figures 1(b,c) illustrate how the mean reaction time depends on the starting point x 0 for a particular choice of a continuously varying heterogeneous surface reactivity κ(θ, φ) shown in Fig. 1(a). When the mean reactivity is weak (κR/D = 1, Fig. 1(b)), S(0|x 0 ) is close to the mean reaction time Shom (0|x 0 ) = R/(3κ) corresponding to homogeneous reactivity κ. Here, multiple failed reaction attempts homogenize the mean reaction time, even though the starting point x 0 lies on the catalytic boundary. In turn, significant deviations from R/(3κ) are observed at a larger mean reactivity κR/D = 10. In this case, the mean reactivity is not representative and heterogeneities start to be more and more important.

B. Diffusion outside a ball

For diffusion in the unbounded domain Ω = {x ∈ R 3 : |x| > R} outside the spherical surface of radius R, the eigenfunctions of the Dirichlet-to-Neumann operator are still given by Eq. (42b), so that the matrix K remains unchanged. In turn, the matrix M is now determined by the eigenvalues

µ (p) nm = -p/D k ′ n (R p/D) k n (R p/D) , ( 48
)
where k n (z) are the modified spherical Bessel function of the second kind. From the known Dirichlet propagator, we compute in Appendix E

V (p) nm (x 0 ) = R -1 k n (r 0 p/D) k n (R p/D) Y mn (θ 0 , φ 0 ) . ( 49
)
As a consequence, our spectral decomposition [START_REF] Bénichou | Kinetics of stochastically gated diffusion-limited reactions and geometry of random walk trajectories[END_REF] fully determines the Laplace-transformed propagator G(x, p|x 0 ). The Laplace-transformed probability density of reaction times is again obtained from Eq. ( 32):

H(p|x 0 ) = √ 4π ∞ n=0 n m=-n h (p) nm k n (r 0 p/D) k n (R p/D) Y mn (θ 0 , φ 0 ), ( 50
) with h (p)
nm given by Eq. ( 46), in which the matrices K and M are determined by Eqs. [START_REF] Singer | Narrow Escape, Part II The circular disk[END_REF][START_REF] Caginalp | Analytical and Numerical Results for an Escape Problem[END_REF]. The Laplacetransformed survival probability is still given by Eq. ( 47), while the mean reaction time is infinite.

According to Eq. ( 34), the Laplace-transformed reaction rate can be obtained by integrating H(p|x 0 ) with the initial concentration of molecules, which for the uniform concentration, c 0 (x 0 ) = c 0 , yields

J(p) = 4πDRc 0 Rµ (p) 00 h (p) 00 p (51)
(note that the same formula with the appropriate matrix M holds for diffusion inside the ball). The prefactor 4πDRc 0 is the Smoluchowski rate to a homogeneous perfectly reactive ball of radius R, whereas the second factor describes the effect of heterogeneous surface reactivity. In the limit p → 0, this expression yields the steady-state reaction rate

J(∞) = 4πDRc 0 h (0) 00 . (52)
For the homogeneous reactivity, our formulas are reduced to that of Collins and Kimball [START_REF] Collins | Diffusion-controlled reaction rates[END_REF], see Appendix E, with

h (0) 00 = 1 1 + D/(κR) . ( 53
)
Figure 2 shows the reaction probability H(0|x 0 ) from Eq. ( 50) (see also Eq. (E8)) on the inert spherical surface covered by ten evenly distributed circular partially reactive targets of angular size ε = 0.2. Even though the starting point x 0 lies on the surface (r 0 = R), the reaction probability is not equal to 1 due to the partial reactivity and eventual failed attempts to react. When the reactivity is large (κR/D = 100, left panel), the reaction probability is close to 1 when the molecule starts at any target, and drops to 0.4 in between two targets. At intermediate reactivity (κR/D = 10, middle panel), the reaction probability is expectedly reduced, as more frequent failed reaction attempts give more chances for the molecule to escape to infinity. This effect is further enhanced at even smaller reactivity κR/D = 1 (right panel). In this regime, there is almost no distinction between weakly reactive targets and the remaining inert surface.

IV. DISCUSSION

We developed a general mathematical description of diffusion-controlled reactions on catalytic surfaces with heterogeneous reactivity κ(s). We showed how the propagator of the diffusion equation with Robin boundary condition can be expressed in terms of the Dirichlet propagator for homogeneous perfectly reactive surface. The latter involves much simpler and more studied Dirichlet boundary condition and thus describes exclusively the first-passage events to the boundary that are independent of the surface reactivity. In other words, the diffusive exploration of the bulk is disentangled from the chemical kinetics on the boundary. As a consequence, the Dirichlet propagator needs to be computed only once for a given geometric configuration, offering a powerful theoretical and numerical tool for investigating the effects of heterogeneous surface reactivity.

Numerical or eventually analytical inversion of the Laplace transform allows one to recover the propagator in time domain. Moreover, the Laplace-transformed propagator itself is important as it describes the steadystate diffusion of molecules which may spontaneously disappear in the bulk with the rate p [START_REF] Yuste | Exploration and trapping of mortal random walkers[END_REF][START_REF] Meerson | Mortality, redundancy, and diversity in stochastic search[END_REF][START_REF] Grebenkov | The escape problem for mortal walkers[END_REF]. Such "mortal walkers" may represent radioactive nuclei, photobleaching fluorophores, molecules in an excited state, metastable complexes, spermatozoa, and other particles subject to spontaneous decay, disintegration, ground state recovery, or death.

When the boundary of the domain is bounded, our general representation yields the spectral decompositions of the propagator and of other important quantities such as the survival probability, the probability density of reaction times, the spread harmonic measure, and the reaction rate. These decompositions involve the eigenvalues and eigenfunctions of the Dirichlet-to-Neumann operator, as well as the associated basis elements of the surface reactivity κ(s). This spectral description brings new insights onto imperfect diffusion-controlled reactions and creates a mathematical basis for formulating and solving optimization and inverse problems on surface reactivity κ(s) (see, e.g., [START_REF] Nguyen | A Spectral Approach to Survival Probability in Porous Media[END_REF][START_REF] Filoche | The toposcopy, a new tool to probe the geometry of an irregular interface by measuring its transfer impedance[END_REF]).

We highlight a similarity between the representation of surface reactivity in the eigenbasis of the Dirichlet-to-Neumann operator and the representation of the bulk reactivity in the Laplace operator eigenbasis studied in [START_REF] Nguyen | A Spectral Approach to Survival Probability in Porous Media[END_REF]. Such matrix representations have proved to be efficient for solving numerically the Bloch-Torrey equation that describes diffusion magnetic resonance imaging (see [START_REF] Grebenkov | Laplacian Eigenfunctions in NMR I. A Numerical Tool[END_REF][START_REF] Grebenkov | NMR Survey of Reflected Brownian Motion[END_REF][START_REF] Grebenkov | From the microstructure to diffusion NMR, and back[END_REF] and references therein). Note that the surface reactivity could also be incorporated via the Laplacian eigenbasis by introducing an infinitely thin reactive boundary layer as discussed in [START_REF] Grebenkov | Residence times and other functionals of reflected Brownian motion[END_REF][START_REF] Grebenkov | Laplacian Eigenfunctions in NMR. II Theoretical Advances[END_REF]. However, the eigenbasis of the Dirichlet-to-Neumann operator acting on the boundary seems to be more natural for dealing with surface reactivity. Most importantly, our spectral description is also valid for exterior problems, for which the spectrum of the Laplace operator is continuous and thus not suitable for such representations; in turn, the spectrum of the Dirichlet-to-Neumann operator on a bounded boundary remains discrete.

We applied the spectral approach to an important example of a spherical surface, for which both the Dirichlet propagator and the eigenbasis of the Dirichlet-to-Neumann operator are known explicitly. In this case, the Robin propagator and the related quantities (such as the probability density of reaction times) are obtained in a semi-analytical form, in which the dependence on the starting and arrival points is fully explicit, whereas the coefficients need to be computed by truncating and inverting an explicitly known matrix. However, the proposed approach is not limited to the spherical boundary. For instance, the case of a hyperplane was partly studied in [START_REF] Grebenkov | Imperfect Diffusion-Controlled Reactions[END_REF][START_REF] Sapoval | New Simple Properties of a Few Irregular Systems[END_REF]; apart from straightforward extensions to disks and cylinders, one can consider more complicated catalytic surfaces formed by multiple non-overlapping spheres, for which the Dirichlet propagator in the steadystate regime was recently investigated in [START_REF] Grebenkov | Semi-analytical computation of Laplacian Green functions in threedimensional domains with disconnected spherical boundaries[END_REF]. In general, the eigenbasis of the Dirichlet-to-Neumann operator M p can be constructed numerically; since M p is independent of the surface reactivity, this construction has to be performed only once for a given catalytic surface.

As mentioned earlier, most former studies focused on the mixed Dirichlet-Neumann boundary value problem describing perfectly reactive targets on an otherwise inert boundary. In spite of its oversimplified character from the chemical point of view, this problem may look simpler from the mathematical point of view. For instance, as both Dirichlet and Neumann boundary conditions are conformally invariant, conformal mapping results in a universal integral representation of the mean first-passage time for planar domains [START_REF] Grebenkov | Universal formula for the mean first passage time in planar domains[END_REF]. In addition, the technique of dual series is more developed for this case [START_REF] Sneddon | Mixed Boundary Value Problems in Potential Theory[END_REF][START_REF] Duffy | Mixed boundary value problems[END_REF]. At the same time, mixed Dirichlet-Neumann condition is the most problematic from the perspective of the present work. Even though the Dirichlet boundary condition can be formally implemented by setting κ(s) = κ on the target and then letting κ go to infinity, an infinitely large jump of reactivity at the border of the target requires elaborate asymptotic analysis. In fact, this limit is in general highly nontrivial because the unbounded Dirichlet-to-Neumann operator M p cannot be neglected as compared to the bounded operator K (representing the reactivity) even as κ → ∞. This situation resembles the asymptotic analysis of the Schrödinger operator -h 2 ∆+V in the semi-classical limit h → 0, where V is a bounded potential. While the application of asymptotic techniques from spectral theory and quantum mechanics to our setting presents an interesting mathematical perspective for future research, our spectral approach is not well suited for studying mixed Dirichlet-Neumann boundary value problems.

Similarly, in the narrow escape limit (when targets are very small), a large number of eigenfunctions of the Dirichlet-to-Neumann operator is needed to accurately represent the multiplication operator K by a truncated matrix K, making numerical computations timeconsuming. More generally, when the shape of the boundary is rather complex or not smooth enough (e.g., containing corners or cusps), the computation of the Dirichlet-to-Neumann eigenfunctions becomes difficult, whereas a large number of eigenfunctions may be needed to project even a smooth surface reactivity. In other words, when the surface reactivity has a substantial projection on a large number of eigenfunctions, the "effective" dimensionality of the matrix K can be large, making the proposed spectral approach less efficient from the numerical point of view. Nevertheless, the present approach can still be advantageous for exterior problems, which are particularly difficult to deal with by other numerical techniques. In this light, the present approach does not substitute conventional techniques but aims to complement them by addressing imperfect diffusioncontrolled reactions on catalytic surfaces with finite continuously varying heterogeneous reactivity.

Appendix A: Alternative representation based on the fundamental solution

In this Appendix, we describe an alternative scheme for representing the propagator in terms of the fundamental solution of the modified Helmholtz equation.

Dirichlet propagator and the Dirichlet-to-Neumann operator

The Laplace-transformed Dirichlet propagator G0 (x, p|x 0 ) and the Dirichlet-to-Neumann operator M p are closely related. On one hand, the action of M p onto a given function can be expressed via Eq. ( 13) in terms of the propagator G0 (x, p|x 0 ) by solving the corresponding Dirichlet boundary value problem. On the other hand, the Dirichlet propagator can be constructed explicitly from the Dirichlet-to-Neumann operator. For this purpose, one can first represent the propagator as

G0 (x, p|x 0 ) = Ĝf (x, p|x 0 ) + ĝ0 (x, p|x 0 ), (A1) where Ĝf (x, p|x 0 ) = exp(-|x -x 0 | p/D) 4πD|x -x 0 | (A2)
is the fundamental solution of the modified Helmholtz equation in three dimensions, (p -D∆) Ĝf (x, p|x 0 ) = δ(xx 0 ), (A3) whereas ĝ0 (x, p|x 0 ) is the regular part of the propagator satisfying, for any fixed

x 0 ∈ Ω, (p -D∆)ĝ 0 (x, p|x 0 ) = 0 (x ∈ Ω), (A4a) ĝ0 (x, p|x 0 ) + Ĝf (x, p|x 0 ) = 0 (x ∈ ∂Ω). (A4b)
Here we use hat instead of tilde in order to distinguish the involved quantities from those in Sec. II C.

The above problem can be solved in a standard way by using the Dirichlet propagator G0 (x, p|x 0 ):

ĝ0 (x, p|x 0 ) = ∂Ω ds ĝ0 (s, p|x 0 ) -D ∂ G0 (x ′ , p|x) ∂n x ′ x ′ =s .
Using the boundary condition (A4b) and substituting the representation (A1), one gets

ĝ0 (x, p|x 0 ) = ∂Ω ds -Ĝf (s, p|x 0 ) (A5) × ĵf (s, p|x) + [DM p Ĝf (•, p|x)](s) , where ĵf (s, p|x 0 ) = -D ∂ Ĝf (x, p|x 0 ) ∂n x x=s (s ∈ ∂Ω) (A6
) is also a fully explicit function, and we used the Dirichletto-Neumann operator M p , acting on Ĝf (s ′ , p|x) as a function of a boundary point s ′ , to represent the normal derivative of ĝ0 (x ′ , p|x). Combining Eqs. (A1, A5), we get the representation of the Dirichlet propagator in terms of the Dirichlet-to-Neumann operator M p and fully explicit functions Ĝf and ĵf .

General Robin boundary value problem

Similarly, for a given function f (s, p) on the boundary ∂Ω, the solution ũ(x, p) of a general Robin boundary value problem

(p -D∆)ũ = 0 (x ∈ Ω), (A7a)
D ∂ ∂n + κ(x) ũ = f (x ∈ ∂Ω), (A7b)
can be obtained by multiplying Eqs. (A3, A7a) by ũ(x, p) and Ĝf (x, p|x 0 ) respectively, subtracting them, integrating over x ∈ Ω, and applying the Green's formula:

ũ(x 0 , p) = ∂Ω ds D Ĝf (s, p|x 0 ) ∂ ũ(x, p) ∂n x x=s + ũ(s, p) ĵf (s, p|x 0 ) . (A8)
This is a standard representation of a solution of the modified Helmholtz equation in terms of the surface integral with the potential Ĝf (x, p|x 0 ) and its normal derivative ĵf (s, p|x 0 ). Here, ũ(x 0 , p) in a bulk point x 0 ∈ Ω is determined by its values and its normal derivative on the boundary. In turn, the Robin boundary condition (A7b) can be expressed in terms of the Dirichlet-to-Neumann operator M p and the operator K of multiplication by κ(x)/D as ũ(s, p)

= 1 D [(M p + K) -1 f ](s), from which Eq. (A8) yields ũ(x 0 , p) = ∂Ω ds Ĝf (s, p|x 0 ) [M p (M p + K) -1 f ](s) + 1 D [(M p + K) -1 f ](s) ĵf (s, p|x 0 ) . (A9)
Since the operator M p is self-adjoint, this solution can also be written as

ũ(x 0 , p) = 1 D ∂Ω ds ĵf (s, p|x 0 ) + [DM p Ĝf (•, p|x 0 )](s) × [(M p + K) -1 f ](s). (A10)
If the boundary ∂Ω is bounded, the spectrum of M p is discrete, and this solution can be written as a spectral decomposition:

ũ(x 0 , p) = 1 D ∞ n,n ′ =0 V (p) n (x 0 ) (M + K) -1 n,n ′ × ∂Ω ds [v (p) n ′ (s)] * f (s, p), (A11)
where

V (p) n (x 0 ) = ∂Ω ds v (p) n (s) ĵf (s, p|x 0 )+Dµ (p) n Ĝf (s, p|x 0 ) , (
A12) and the matrices M and K are defined in Eq. (29a). In particular, the Laplace-transformed propagator G(x, p|x 0 ) for the Robin boundary value problem [START_REF] Bénichou | Geometry-controlled kinetics[END_REF] can be written as

G(x, p|x 0 ) = Ĝf (x, p|x 0 ) + 1 D ∞ n,n ′ =0 V (p) n (x 0 ) × (M + K) -1 n,n ′ [ Û (p) n ′ (x)] * , ( A13
)
where

Û (p) n (x) = ∂Ω ds v (p) n (s) ĵf (s, p|x) -κ(s) Ĝf (s, p|x) .
(A14) In contrast to Eq. ( 27), this representation is based on the explicitly known fundamental solution Ĝf (x, p|x 0 ) and does not involve the Dirichlet propagator G0 (x, p|x 0 ). As a consequence, all the deduced spectral decompositions rely uniquely on the eigenbasis of the Dirichlet-to-Neumann operator. While the representations ( 27) and (A13) are equivalent and complementary to each other, we keep using the former one due to its simpler form and clearer probabilistic interpretation.

where the last integral reflects the normalization of the harmonic measure density j0 (s, 0|x 0 ) for a bounded domain.

For an unbounded domain, a nonzero constant cannot be a solution of the Laplace equation ∆u = 0 with u |∂Ω = 1 due to the regularity condition u(x) → 0 as |x| → ∞. The function M 0 1 is thus not zero, and the smallest eigenvalue µ (0) 0 is strictly positive. As a consequence, the second term in Eq. (B3) does not vanish, while the first term is not equal to 1. In other words, the reaction probability H(0|x 0 ) is in general less than 1 due to the possibility for a molecule to escape at infinity. In this case, H(t|x 0 ) should be renormalized by H(0|x 0 ) to get the conditional probability density of reaction times.

Laplace-transformed reaction rate

We briefly discuss how Eq. ( 34) for the Laplacetransformed reaction rate J(p) can be further simplified when the initial concentration is uniform: c 0 (x 0 ) = c 0 .

Integrating Eq. (6a) for the propagator G0 (x, p|x 0 ) over x ∈ Ω yields

p Ω dx 0 G0 (x, p|x 0 ) = 1 + D ∂Ω dx 0 ∂ G0 (x, p|x 0 ) ∂n x0 , (B4)
where we exchanged x and x 0 due to the symmetry of the propagator. Applying the normal derivative at a boundary point x = s ∈ ∂Ω and multiplying by -D, we get

Ω dx 0 j0 (s, p|x 0 ) = D p ∂Ω ds 0 ∂ j0 (s, p|x 0 ) ∂n x0 x0=s0 .
(B5) Multiplying this relation by a function f (s) and integrating over s ∈ ∂Ω, we have

Ω dx 0 ∂Ω ds j0 (s, p|x 0 ) f (s) = D p ∂Ω ds 0   ∂ ∂n x0 ∂Ω ds j0 (s, p|x 0 ) f (s)   x0=s0 = D p ∂Ω ds 0 [M p f ](s 0 ),
where the order of integrals was exchanged. As it is satisfied for any f (s), we conclude that

[M p 1](s) = p D Ω dx 0 j0 (s, p|x 0 ). (B6) Setting f (s) = v (p)
n (s), we also deduce

Ω dx 0 V (p) n (x 0 ) = D p µ (p) n ∂Ω ds v (p) n (s). ( B7
)
This expression allows us to compute the integral in Eq. ( 34), yielding

J(p) = c 0 D p ∞ n,n ′ =0   ∂Ω ds v (p) n (s)   (B8) × M(M + K) -1 K n,n ′   ∂Ω ds [v (p) n ′ (s)] *   .
Appendix C: Diffusion inside a ball Solutions of Dirichlet boundary value problems for the modified Helmholtz equation in a ball and the related operators are well known. For the sake of clarify and completeness, we summarize the main "ingredients" involved in our spectral decompositions.

To determine the eigenbasis of the Dirichlet-to-Neumann operator in a ball of radius R, Ω = {x ∈ R 3 : |x| < R}, one simply notes that a general solution of the modified Helmholtz equation (p -D∆)ũ = 0 can be written in spherical coordinates (r, θ, φ) as

ũ(x, p) = ∞ n=0 m m=-n a mn i n (r p/D) Y mn (θ, φ), (C1)
where a mn are unknown coefficients,

i n (z) = π/2 I n+1/2 (z) √ z (C2)
are the modified spherical Bessel functions of the first kind, and

Y mn (θ, φ) = c nm P m n (cos θ) e imφ (C3)
are the spherical harmonics, with P m n (x) being the associated Legendre functions and c nm the normalization coefficients:

c nm = 2n + 1 4π (n -m)! (n + m)! . ( C4
)
As the normal derivative of ũ on the boundary involves only the radial coordinate and does not affect Y mn (θ, φ), the eigenvalues and eigenfunctions of the Dirichlet-to-Neumann operator are

µ (p) nm = p/D i ′ n (R p/D) i n (R p/D) , (C5a)
v nm (θ, φ) = 1 R Y mn (θ, φ), (C5b)
where prime denotes the derivative with respect to the argument. As stated in the main text, the double index nm is employed to enumerate the eigenfunctions as well as the elements of the matrices M and K. Note that the eigenvalues that determine the matrix M in Eq. (29a), do not depend on the index m and thus are of multiplicity 2n + 1. In turn, the eigenfunctions v nm do not depend on the parameter p that will simplify further expressions. The explicit computation of the matrix K from Eq. ( 43) is discussed in Appendix D. For a ball, the Dirichlet propagator is known explicitly

G 0 (x, t|x 0 ) = 1 2πR 3 ∞ n=0 (2n + 1)P n (x • x 0 ) |x| |x 0 | (C6) × ∞ k=0 e -Dtα 2 nk /R 2 [j ′ n (α nk )] 2 j n (α nk r/R) j n (α nk r 0 /R),
where α nk are the positive zeros (enumerated by the index k = 0, 1, 2, . . .) of the spherical Bessel functions j n (z) of the first kind, and we used the addition theorem for spherical harmonics to evaluate the sum over the index m:

P n (x • x 0 ) |x| |x 0 | = 4π n m=-n Y mn (θ 0 , φ 0 )Y * mn (θ, φ) 2n + 1 (C7) = 4πR 2 2n + 1 n m=-n v nm (θ 0 , φ 0 ) v * nm (θ, φ),
where P n (z) are the Legendre polynomials. The Laplace transform of Eq. (C6) reads

G0 (x, p|x 0 ) = 1 2πR 3 ∞ n=0 (2n + 1)P n (x • x 0 ) |x| |x 0 | × ∞ k=0 j n (α nk r/R) j n (α nk r 0 /R) (Dα 2 nk /R 2 + p)[j ′ n (α nk )] 2 . (C8) As (x•x0) |x| |x0|
is the cosine of the angle between the vectors x and x 0 , it does not depend on the radial coordinates r and r 0 . We get thus j0 (s,

p|x 0 ) = - 1 2πR 2 ∞ n=0 (2n + 1)P n (s • x 0 ) |s| |x 0 | × ∞ k=0 α nk j n (α nk r 0 /R) (α 2 nk + pR 2 /D)j ′ n (α nk ) , (C9)
from which

V (p) nm (x 0 ) = -2v nm (θ 0 , φ 0 ) ∞ k=0 α nk j n (α nk r 0 /R) (α 2 nk + pR 2 /D)j ′ n (α nk ) = v nm (θ 0 , φ 0 ) i n (r 0 p/D) i n (R p/D) , (C10)
where we used the summation formula over zeros α nk (see Eq. (S9) from Table 3 of Ref. [START_REF] Grebenkov | A physicist's guide to explicit summation formulas involving zeros of Bessel functions and related spectral sums[END_REF]). Expressions (C8, C9, C10) determine the Laplace-transformed propagator G(x, p|x 0 ) of the Robin boundary value problem in the semi-analytical form [START_REF] Bénichou | Kinetics of stochastically gated diffusion-limited reactions and geometry of random walk trajectories[END_REF], in which the dependence on points x 0 and x is fully explicit, whereas the computation of the coefficients involves a numerical inversion of the matrix M + K. Similarly, we deduce semi-analytical expressions for the Laplace-transformed probability density of reaction times and the spread harmonic measure presented in the main text. For instance, as the eigenfunctions v nm are orthogonal to v 00 (θ, φ) = 1/( √ 4πR), the sum in Eq. ( 33) is reduced to a single term, yielding Eq. [START_REF] Cheviakov | An Asymptotic Analysis of the Mean First Passage Time for Narrow Escape Problems: Part II: The Sphere[END_REF].

We also compute the mean reaction time S(0|x 0 ) by evaluating the limit p → 0 of the Laplace-transformed survival probability from Eq. [START_REF] Cheviakov | Mathematical modeling and numerical computation of narrow escape problems[END_REF]. In this limit, one gets

µ (0) nm = n/R, V (0) nm (x 0 ) = v nm (θ 0 , φ 0 )(r 0 /R) n , (C11) so that S(0|x 0 ) = √ 4π ∞ n=0 n m=-n Y mn (θ 0 , φ 0 )(r 0 /R) n × R 2 -r 2 0 4D(n + 3/2) h (0) mn - h (p) mn dp p=0 . ( C12
)
The last term can be evaluated explicitly as

h (p) mn dp p=0 = -(M 0 + K) -1 M 1 (M 0 + K) -1 K nm,00 , (C13)
where M 0 and M 1 are diagonal matrices obtained by expanding the elements of M into powers p:

[M 0 ] nm,n ′ m ′ = δ n,n ′ δ m,m ′ n/R and [M 1 ] nm,n ′ m ′ = δ n,n ′ δ m,m ′ R/D
2n+3 . For a homogeneous reactivity, κ(θ, φ) = κ, Eqs. [START_REF] Crank | The Mathematics of Diffusion[END_REF][START_REF] Pillay | An Asymptotic Analysis of the Mean First Passage Time for Narrow Escape Problems: Part I: Two-Dimensional Domains[END_REF] imply

Hhom (p|x 0 ) = κ i 0 (r 0 p/D) √ pD i ′ 0 (R p/D) + κ i 0 (R p/D) (
C14) (since i 0 (z) = sinh(z)/z, one can further simplify this expression). In turn, Eq. [START_REF] Zwanzig | Diffusion-controlled ligand binding to spheres partially covered by receptors: an effective medium treatment[END_REF] gives

Jhom (p) = 4πDRc 0 i 0 (R p/D) R p/D i 1 (R p/D) + D κR -1
, (C15) and its inverse Laplace transform yields an infinite sum of exponentially decaying functions with the rates determined by the poles of this expression (see [START_REF] Redner | A Guide to First Passage Processes[END_REF][START_REF] Carslaw | Conduction of Heat in Solids[END_REF][START_REF] Crank | The Mathematics of Diffusion[END_REF]). Finally, Eq. (C12) yields the classical result

Shom (0|x 0 ) = R 2 -r 2 0 6D + R 3κ . ( C16
)

Numerical validation

To illustrate the quality of our semi-analytical solution, we look at the Laplace-transformed probability density H(p|x 0 ) that satisfies the boundary value problem

(p -D∆ x0 ) H(p|x 0 ) = 0 (x 0 ∈ Ω), D ∂ ∂n x0 + κ(x 0 ) H(p|x 0 ) = κ(x 0 ) (x 0 ∈ ∂Ω),
with ∆ x0 acting on x 0 . We set κ(θ, φ) = κ Θ(ε-θ) to describe a single partially reactive circular target of angular size ε and reactivity κ, located at the North pole (here Θ(z) is the Heaviside function). The axial symmetry of this geometric setting allows one to reduce the original three-dimensional problem to a two-dimensional one on the rectangle [0, R] × [0, π] in the coordinates (r, θ). We solve this problem by using a finite element method implemented in Matlab PDE toolbox. The computational domain was meshed with the constraint on the largest mesh size to be 0.01. For the sake of simplicity, we fix the starting point at the origin, in which case Eq. ( 45) is reduced to

H(p|0) = h (p) 00 R p/D sinh(R p/D) , (C17) where h (p)
00 is given by Eq. ( 46) and computed with the matrices K and M truncated to the size n max = 20 and constructed from Eqs. (D12, 29a). Figure 3 shows an excellent agreement between this semi-analytical form and the FEM solution for both a small target of angular size ε = 0.1 (with the surface fraction σ = (1 -cos ε)/2 ≈ 0.0025) and a large target of angular size ε = 1 (with σ ≈ 0.23), and different reactivities. The key element of the spectral approach is the possibility to disentangle first-passage diffusive steps from the heterogeneous reactivity which is incorporated via the matrix K. In this Appendix, we compute this matrix for several most common settings on the spherical boundary.

General setting

In general, the reactivity κ(θ, φ) can be expanded over the complete basis of spherical harmonics,

κ(θ, φ) = ∞ n=0 n m=-n κ nm Y mn (θ, φ), (D1)
with coefficients κ nm . When this expansion can be truncated at a low order n * , one can compute the elements of the matrix K explicitly, without numerical quadrature in Eq. [START_REF] Singer | Narrow Escape, Part II The circular disk[END_REF], by using the following identity

π 0 dθ sin θ 2π 0 dφ Y m1n1 (θ, φ) Y m2n2 (θ, φ) Y m3n3 (θ, φ) = (2n 1 + 1)(2n 2 + 1)(2n 3 + 1) 4π × n 1 n 2 n 3 m 1 m 2 m 3 n 1 n 2 n 3 0 0 0 , (D2)
where

n 1 n 2 n 3 m 1 m 2 m 3
is the Wigner 3j symbol. We note that the truncation order n max should significantly exceed n * to ensure accurate computations.

Axially symmetric problems

When the reactivity is axially symmetric, κ(θ, φ) = κ(θ), the integral over φ in Eq. ( 43) yields 2πδ m,m ′ , and the matrix K has a block structure. If in addition one is interested in axially symmetric quantities (e.g., H(p|x 0 ) which does not depend on φ 0 due to the axial symmetry), it is sufficient to construct a reduced version of the matrix K by eliminating repeated lines and rows and keeping only the elements with m = m ′ = 0:

K n0,n ′ 0 = (n + 1/2)(n ′ + 1/2) (D3) × π 0 dθ sin θ κ(θ) D P n (cos θ) P n ′ (cos θ).
From the numerical point of view, this drastically speeds up computations because the size of the matrix K, truncated to the order n max , becomes (n max + 1) × (n max + 1) instead of (n max + 1) 2 × (n max + 1) 2 in the general setting. Semi-analytical expressions also become simpler, e.g., Eqs. ( 45) and ( 50) read respectively

H(p|x 0 ) = ∞ n=0 √ 2n + 1 h (p) n0 i n (r 0 p/D) i n (R p/D) P n (cos θ 0 ) (D4) and H(p|x 0 ) = ∞ n=0 √ 2n + 1 h (p) n0 k n (r 0 p/D) k n (R p/D) P n (cos θ 0 ), (D5) with h (p
) n0 given by Eq. ( 46). If an expansion of the reactivity κ(θ) over the complete basis of Legendre polynomials is known,

κ(θ) = ∞ n=0 κ n P n (cos θ), (D6)
then the elements of the matrix K can be computed by using the identity

1 -1 dx P n1 (x) P n2 (x) P n (x) = 2 n 1 n 2 n 0 0 0 2 , (D7)
which follows from Eq. (D2). As the selection rule for Wigner 3j-symbols requires that |n 1 -n 2 | ≤ n ≤ n 1 + n 2 , the truncation of the expansion (D6) at the order n * implies that the matrix K has at most n * subdiagonals above and below the main diagonal that simplifies the construction of this matrix. One advantage of the representation (D6) is that the average reactivity is equal to κ 0 and is independent of κ n with n ≥ 1 due to the orthogonality of Legendre polynomials. We emphasize however that the above simplified construction is not sufficient for computing the Laplacetransformed propagator G(x, p|x 0 ) which is not axially symmetric. In fact, Eq. ( 27) involves the coefficients [(M + K) -1 ] nm,n ′ m ′ , whose computation requires all the elements K nm,n ′ m ′ even for axially symmetric reactivity, and it is not reducible to that with the elements K n0,n ′ 0 . In this case, the general scheme from Sec. D 1 should be used.

Single circular target

To model a single circular partially reactive target of angular size ε at the North pole (with the remaining inert boundary), one sets κ(θ, φ) = κ Θ(ε -θ), (D8) so that Eq. (D3) yields

K n0,n ′ 0 = κ D (n + 1/2)(n ′ + 1/2) 1 cos ε dx P n (x) P n ′ (x).
(D9) To compute explicitly the matrix K, one can use the Adams-Neumann's product formula (see [START_REF] Al-Salam | On the product of two Legendre polynomials[END_REF]):

P n (x)P n ′ (x) = min{n,n ′ } k=0 B k nn ′ P n+n ′ -2k (x), (D10)
where

B k nn ′ = A k A n-k A n ′ -k A n+n ′ -k 2n + 2n ′ -4k + 1 2n + 2n ′ -2k + 1 , (D11) with A k = Γ(k+1/2) √ πΓ(k+1) (with A 0 = 1). We get thus K n0,n ′ 0 = κ D (n + 1/2)(n ′ + 1/2) × (D12) min{n,n ′ } k=0 B k nn ′ P n+n ′ -2k-1 (cos ε) -P n+n ′ -2k+1 (cos ε) 2(n + n ′ -2k) + 1 ,
where we used the identity for n ≥ 0 An explicit formula for K is also easily deducible for multiple latitudinal stripes. The domain is still axially symmetric and one just needs to sum up contributions from each stripe, relying on the explicit integral of P n (x) in Eq. (D13).

Multiple targets of circular shape

The matrix K can also be computed explicitly for multiple partially reactive non-overlapping targets of circular shape. In fact, the additivity of the integral in Eq. [START_REF] Singer | Narrow Escape, Part II The circular disk[END_REF] implies that contributions for all targets are just summed up. We consider thus the contribution of the i-th target Γ i of angle ε i , reactivity κ i , and the angular coordinates (θ i , φ i ) for its center. It is convenient to apply the rotational addition theorem for spherical harmonics to rotate the coordinate system [START_REF] Steinborn | Rotation and Translation of Regular and Irregular Solid Spherical Harmonics[END_REF]:

Y mn (θ ′ , φ ′ ) = n m ′ =-n [D n mm ′ (φ i , θ i , φ i )] * Y m ′ n (θ, φ),
(D14) where D n mm ′ (α, β, γ) is the Wigner D-matrix describing the rotation by Euler angles (α, β, γ). As a consequence, the i-th contribution to the matrix K reads

K (i) n1m1,n2m2 = κ i DR 2 Γi ds ′ Y * m1n1 (θ ′ , φ ′ ) Y m2n2 (θ ′ , φ ′ ) = κ i DR 2 n1 m ′ 1 =-n1 D n1 m1m ′ 1 (φ i , θ i , φ i ) n2 m ′ 2 =-n2 ×[D n2 m2m ′ 2 (φ i , θ i , φ i )] * Γ0 ds Y * m ′ 1 n1 (θ, φ)Y m ′ 2 n2 (θ, φ),
where Γ 0 is the i-th target rotated to be centered around the North pole. To proceed, one can express the product of two spherical harmonics as

Y m ′ 1 n1 (θ, φ) Y m ′ 2 n2 (θ, φ) (D15) = n1+n2 n=|n1-n2| B n m ′ 1 n1m ′ 2 n2 Y (m ′ 1 +m ′ 2 )n (θ, φ)
(which follows from Eq. (D2)), where

B n m ′ 1 n1m ′ 2 n2 = (2n + 1)(2n 1 + 1)(2n 2 + 1) 4π (D16) ×(-1) m ′ 1 +m ′ 2 n 1 n 2 n m ′ 1 m ′ 2 -m ′ 1 -m ′ 2 n 1 n 2 n 0 0 0 , with n 1 n 2 n m 1 m 2 m
being again the Wigner 3-j symbols [START_REF] Brink | Angular Momentum[END_REF], and we employ the convention that Y mn (θ, φ) ≡ 0 if |m| > n. Using the identity

Y * mn (θ, φ) = (-1) m Y (-m)n (θ, φ),
the above formula yields

Y * m ′ 1 n1 (θ, φ) Y m ′ 2 n2 (θ, φ) (D17) = (-1) m ′ 2 n1+n2 n=|n1-n2| B n (-m ′ 1 )n1m ′ 2 n2 Y (m ′ 2 -m ′ 1 )n (θ, φ), from which
K (i) n1m1,n2m2 = √ π κ i D n1+n2 n=|n1-n2| P n-1 (cos ε i ) -P n+1 (cos ε i ) √ 2n + 1 × min{n1,n2} m=-min{n1,n2} (-1) m D n1 m1m (φ i , θ i , φ i )[D n2 m2m (φ i , θ i , φ i )] * B n (-m)n1mn2 , (D18)
where the integral over φ yielded 2πδ m ′ 1 ,m ′ 2 that removed one sum, while the integral of P n (x) was evaluated from Eq. (D13). We get therefore a fully explicit expression for the contribution of the i-th target to the matrix K. One can thus compute the Laplace-transformed propagator in the semi-analytical form, as for a single target. where k n (z) are the modified spherical Bessel functions of the second kind:

k n (z) = 2/π K n+1/2 (z) √ z . ( E2
)
As for the interior problem, the eigenvalue µ (E5) where we used the Wronskian i ′ n (z)k n (z) -k ′ n (z)i n (z) = 1/z 2 . We compute then

V (p) nm = v nm (θ 0 , φ 0 ) k n (r 0 p/D) k n (R p/D) . (E6)
According to our spectral decomposition [START_REF] Bénichou | Kinetics of stochastically gated diffusion-limited reactions and geometry of random walk trajectories[END_REF], Eqs. (E4, E6), together with the matrices M and K, fully determine the Laplace-transformed propagator G(x, p|x 0 ). Similarly, we deduce the Laplace-transformed probability density of reaction times and the reaction time presented in the main text.

In the limit p → 0, one has µ (0) nm = (n + 1)/R, V (0) nm (x 0 ) = v nm (θ 0 , φ 0 )(R/r 0 ) n+1 , (E7) from which

H(0|x 0 ) = √ 4π ∞ n=0 n m=-n h (0)
nm (R/r 0 ) n+1 Y mn (θ 0 , φ 0 ) (E8) is the probability of reaction on the ball, and h (p) nm are defined by Eq. [START_REF] Cheviakov | An Asymptotic Analysis of the Mean First Passage Time for Narrow Escape Problems: Part II: The Sphere[END_REF]. In contrast to bounded domains, for which this probability was equal to 1 (see Eq. (B3)), the transient character of Brownian motion in three dimensions makes this probability less than 1. Rewriting Eq. [START_REF] Cheviakov | An Asymptotic Analysis of the Mean First Passage Time for Narrow Escape Problems: Part II: The Sphere[END_REF] as

h (0) nm = δ n,0 δ m,0 - 1 R (M p=0 + K) -1 nm,00 , (E9)
one can split H(0|x 0 ) into two terms, in which the first term R/r 0 is the hitting probability to a perfectly reactive ball, while the second term accounts for partial heterogeneous reactivity. The limit p → 0 also determines the long-time behavior of Eq. ( 51) that yields Eq. ( 52) for the steady-state reaction rate J(∞).

For homogeneous reactivity, κ(s) = κ, M + K is a diagonal matrix and thus only the term with m = n = 0 survives in Eq. (E8), yielding the classical result for the hitting probability of a homogeneous partially reactive ball: Hhom (0|x 0 ) = R r 0 where we used k 0 (z) = e -z /z. The Laplace inversion recovers the result by Collins and Kimball [START_REF] Collins | Diffusion-controlled reaction rates[END_REF] H hom (t|x

0 ) = κ r 0 exp - (r 0 -R) 2 4Dt R √ πDt (E12) -1 + κR D erfcx r 0 -R √ 4Dt + 1 + κR D √ Dt R ,
where erfcx(x) = e x 2 erfc(x) is the scaled complementary error function (see also discussion in [START_REF] Grebenkov | Strong defocusing of molecular reaction times results from an interplay of geometry and reaction control[END_REF]). In the limit κ → ∞, this expression reduces to

H hom (t|x 0 ) = R r 0 r 0 -R √ 4πDt 3 exp - (r 0 -R) 2 4Dt . (
In the short-time limit t → 0, the reaction rate approaches a constant, J hom (0) = 4πκR 2 c 0 , which corresponds to reaction-limited kinetics (note that J hom (0) > J hom (∞)). In the limit κ → ∞, one retrieves the Smoluchowski result:

J hom (t) = 4πDRc 0 1 + √ R √ πDt . ( E16
)
We note that the above analysis can be easily extended to the case when the spherical target of radius R is surrounded by an outer reflecting concentric sphere of radius R o . The eigenfunctions of the Dirichlet-to-Neumann operator remain unchanged, whereas the eigenvalues become . (E17)

In the limit R o → ∞, one retrieves Eq. (E1) for the exterior of a ball. As p → 0, one also gets

µ (0) n = n + 1 R 1 -(R/R o ) 2n+1 1 + (1 + 1/n)(R/R o ) 2n+1 . ( E18
)
As a consequence, one can easily extend the former results to this setting.

  orthogonal basis in L 2 (∂Ω):

FIG. 1 :

 1 FIG. 1: Diffusion inside a ball of radius R. (a) Heterogeneous surface reactivity κ(θ, φ) = κ(1 + cY2,3(θ, φ) + cY-2,3(θ, φ)), with κR/D = 1 and c = 1.2728 (this value ensures the positivity of κ(θ, φ)). (b,c) Mean reaction time S(0|x0), rescaled by Shom (0|x0) = R/(3κ), as a function of the starting point x0 = (r0, θ0, φ0) with r0 = R, for such a reactivity κ(θ, φ), with κR/D = 1 (b) and κR/D = 10 (c). The matrix K was computed with the truncation order nmax = 20 as described in Appendix D 1.

FIG. 2 :

 2 FIG. 2: Diffusion outside a ball of radius R. The reaction probability H(0|x0) on the spherical surface of radius R as a function of the starting point x0 = (r0, θ0, φ0) with r0 = R. Ten circular targets of angular size ε = 0.2 (shown by thin black circles) are evenly distributed on the surface, with three values of reactivity: κR/D = 100 (left), κR/D = 10 (middle), and κR/D = 1 (right). The matrix K was computed with the truncation order nmax = 20 as described in Appendix D 4.

  Appendix D: Computation of the matrix K

FIG. 3 :

 3 FIG.3: Laplace-transformed probability density H(p|0) of reaction times on a partially reactive circular target of reactivity κ and angular size ε, located on the inert spherical surface of radius R, for a molecule started from the origin, with ε = 0.1 (a) and ε = 1 (b). Lines show the semi-analytical solution (C17), in which h (p) 00 was found from Eq. (46) with the matrices M and K truncated at nmax = 20. Symbols present a FEM numerical solution with the maximal mesh size of 0.01.

  n (x) = P n+1 (b) -P n-1 (b) -P n+1 (a) + P n-1 (a)2n + 1 (D13) (with the convention P -1 (x) = 1).

  Appendix E: Diffusion outside a ballFor diffusion in the unbounded domain Ω = {x ∈ R 3 : |x| > R} outside a ball of radius R, the eigenfunctions of the Dirichlet-to-Neumann operator are still given by Eq. (42b), whereas the eigenvalues areµ (p) nm = -p/D k ′ n (R p/D) k n (R p/D) ,(E1)

  depend on m and has thus the multiplicity 2n + 1. Note that µ (p) nm are just polynomials of R p/D, e.g., µ×P n (x • x 0 ) |x| |x 0 | i n (R p/D) k n (R p/D) k n (r p/D)k n (r 0 p/D).Note that the fundamental solution (the first term) also admits the decomposition: n(x • x 0 ) |x| |x 0 | × k n (r 0 p/D) i n (r p/D) (E3)for r < r 0 (and r 0 is exchanged with r for r > r 0 ), where we applied the addition theorem (C7) for spherical harmonics. One gets thenG0 (x, p|x 0 ) = p/D 4πD ∞ n=0 (2n + 1)P n (x • x 0 ) |x| |x 0 | × (E4) k n (r 0 p/D) i n (r p/D) -k n (r p/D) i n (R p/D) k n (R p/D)for r < r 0 . In particular, one deducesj0 (s, p|x 0 ) = ∞ n=0 2n + 1 4πR 2 P n (s • x 0 ) |s| |x 0 | k n (r 0 p/D) k n (R p/D) ,

  retrieves in time domain the reaction rate derived by Collins and Kimball[START_REF] Collins | Diffusion-controlled reaction rates[END_REF] J hom (t) J hom (∞)

µ

  (p) n =p/D k ′ n (R o p/D) i ′ n (R p/D) -i ′ n (R o p/D) k ′ n (R p/D) k ′ n (R o p/D) i n (R p/D) -i ′ n (R o p/D) k n (R p/D)

Appendix B: Technical derivations 1. Reaction probability

The reaction probability can be obtained by integrating the probability density H(t|x 0 ) of reaction times over t from 0 to infinity, giving H(0|x 0 ). For any bounded domain, a diffusing molecule cannot avoid the reaction event so that H(0|x 0 ) = 1, ensuring the correct normalization of the probability density H(t|x 0 ). This property can be checked directly from our spectral representation [START_REF] Shoup | Diffusioncontrolled bimolecular reaction rates. The effect of rotational diffusion and orientation constraints[END_REF]. Setting p = 0 yields the reaction probability

For any bounded domain, the Laplace equation ∆u = 0 with u |∂Ω = 1 on the boundary has the constant solution, u ≡ 1, so that a constant function 1 on the boundary is an eigenfunction of the Dirichlet-to-Neumann operator, v (0) 0 (s) = |∂Ω| -1/2 , corresponding to µ (0) 0 = 0. As a consequence, the second sum over n ′ in Eq. (B1) vanishes due to the orthogonality of eigenfunctions, yielding