Denis S Grebenkov

Sergey D Traytak

Semi-analytical computation of Laplacian Green functions in three-dimensional domains with disconnected spherical boundaries

Keywords: Green function, Laplace operator, boundary value problem, diffusion-reaction, semi-analytical solution

Semi-analytical computation of Laplacian Green

functions in three-dimensional domains with disconnected spherical boundaries

Introduction

Diffusion-reaction processes in porous materials and biological media play an important role in various fields, from physics to chemistry, biology and ecology [START_REF] Rice | Diffusion-limited reactions[END_REF][START_REF] Calef | Diffusion-controlled reactions[END_REF][START_REF] Weiss | Overview of theoretical models for reaction rates[END_REF][START_REF] Zhou | Rate theories for biologists[END_REF].

Email addresses: denis.grebenkov@polytechnique.edu (Denis S. Grebenkov), sergtray@mail.ru (Sergey D. Traytak) Preprint submitted to Journal of Computational Physics February 12, 2018 The geometric structure of these media is often modeled by packs of spheres. Some spheres can be just inert reflecting obstacles to diffusing particles, the others can fully or partially absorb the particles, while the third ones allow for diffusive exchange between interior and exterior compartments. In the stationary regime, the local concentration of diffusing particles, n(x), obeys the Laplace equation, ∇ 2 n(x) = 0, subject to appropriate boundary conditions. Similar boundary value problems arise in various sciences such as heat transfer [START_REF] Carslaw | Conduction of Heat in Solids[END_REF][START_REF] Samarskii | The Finite Difference Methodology[END_REF], electrostatics [START_REF] Jackson | Classical Electrodynamics[END_REF], hydrodynamics [START_REF] Milne-Thomson | Theoretical Hydrodynamics[END_REF], geophysics [START_REF] James | Transformation of spherical harmonics under change of reference frame[END_REF], and probability theory [START_REF] Gardiner | Handbook of stochastic methods[END_REF][START_REF] Keilson | Green's function methods in probability theory[END_REF]. Although the Laplace equation is probably the most well studied partial differential equation (PDE), its explicit analytical solutions are available only for a very limited number of three-dimensional domains [START_REF] Carslaw | Conduction of Heat in Solids[END_REF][START_REF] Crank | The Mathematics of Diffusion[END_REF]. Among them one usually distinguishes solutions obtained by separation of variables in separable curvilinear coordinate systems that are determined by the Euclidean symmetry group of the Laplace equation [START_REF] Miller | Symmetry and separation of variables[END_REF]. The most common examples are a sphere and a circular cylinder. In more complicated but also more practically relevant cases, one has to resort to numerical methods. Except for Monte Carlo simulations, essentially all numerical methods aim to reduce the PDE to an infinite system of linear algebraic equations (ISLAE) that is then solved numerically. The efficiency of a numerical method depends thus on the chosen reduction scheme. For instance, in a finite element method (FEM), the PDE is projected onto basic functions which are piecewise polynomials on each element of a meshed computational domain. The unknown coefficients in front of these functions are then determined by solving a system of linear equations. Without relying on specific geometrical properties of the domain, the FEM is a powerful tool to solve general PDEs in general bounded domains [START_REF] Reddy | An Introduction to the Finite Element Method[END_REF][START_REF] Eun | Influence of neighboring reactive particles on diffusion-limited reactions[END_REF]. In turn, a (much) higher computational efficiency is expected for a method that is specifically adapted to the geometrical structure of the domain. In particular, when the domain has disconnected spherical boundaries, one can profit from the underlying local spherical symmetries to build up more efficient but less generic methods.

In this paper, we revisit the so-called generalized method of separation of variables (GMSV) that goes back to Rayleigh's seminal paper on the conductivity of heat and electricity in a medium with cylindrical or spherical obstacles arranged in a rectangular array [START_REF] Rayleigh | On the influence of obstacles arranged in rectangular order upon the properties of a medium[END_REF]. In its modern form, the GMSV was thoroughly developed for studying diffraction of electromagnetic waves on surfaces of several bodies [START_REF] Ivanov | Diffraction of electromagnetic waves on two bodies[END_REF] and then applied in various fields. It is striking how many names were given to the method under consideration by different authors: "the method of addition theorems" [START_REF] Erofeenko | Addition theorems and solutions to the boundary value problems of mathematical physics[END_REF], "the method of reduction to the ISLEA" [START_REF] Alexandrov | Problems of the mechanics of a continuous medium with mixed boundary conditions[END_REF], "the method of irreducible Cartesian tensors" [START_REF] Traytak | Theory of recondensation of N drops[END_REF], "the method based on the theory of multipole expansions" [START_REF] Clercx | Many-body electrostatic interactions in electrorheological fluids[END_REF], "the generalized Fourier method" [START_REF] Nikolaev | The generalized Fourier method in spacial problems of the elasticity theory for canonical multi-connected bodies[END_REF], "the Rayleigh multipole method" [START_REF] Yardley | Addition formulas and the Rayleigh identity for arrays of elliptical cylinders[END_REF], "the method of twin multipole expansions" [START_REF] Tsao | Competitive diffusion into two reactive spheres of different reactivity and size[END_REF], "the direct method of re-expansion" [START_REF] Traytak | Methods for solution of the boundary value problems in domains with disconnected boundary[END_REF], "a twin spherical expansions method" (or just "twin expansions technique") [START_REF] Mcdonald | Diffusion and reaction for a spherical source and sink[END_REF], "the method of a bispherical expansion" [START_REF] Mcdonald | Competitive interaction between two different spherical sinks[END_REF], "the multipole re-expansion method" [START_REF] Gumerov | Computation of scattering from clusters of spheres using the fast multipole method[END_REF], "the multipole expansion method" [START_REF] Kushch | Multipole expansion method in micromechanics of composites[END_REF], and a particular case of "the method of series" [START_REF] Guz | Diffraction of the elastic waves in multi-connected bodies[END_REF]. The GMSV has been successfully applied in the elasticity theory [START_REF] Miyamoto | On the problem of the theory of elasticity for a region containing more than two spherical cavities[END_REF], heat transfer [START_REF] Jefferey | Conduction through a random suspension of spheres[END_REF], diffraction theory [START_REF] Guz | Diffraction of the elastic waves in multi-connected bodies[END_REF][START_REF] Linton | Lattice sums for the Helmholtz equation[END_REF] and other branches of mathematical physics [START_REF] Erofeenko | Addition theorems and solutions to the boundary value problems of mathematical physics[END_REF]. In the diffusion context, Mitra and later Goodrich were the first to apply the GMSV to find the steady-state diffusion field in the vicinity of two identical ideal spherical sinks (drops) [START_REF] Mitra | A new method of solution of the boundary value problems of Laplace's equation relating to two spheres -Part-I[END_REF][START_REF] Goodrich | On the diffusion field in the neighborhood of two identical spheres[END_REF]. Felderhof investigated diffusion-controlled reactions in regular and random arrays of static ideal spherical sinks [START_REF] Felderhof | Wigner solids and diffusion controlled reactions in a regular array of spheres[END_REF][START_REF] Mattern | Rate of diffusion-controlled reactions in a random array of spherical sinks[END_REF]. Later Venema and Bedeaux applied the GMSV to study similar problems for a periodic array of penetrable spherical sinks [START_REF] Venema | Reaction-diffusion on a periodic array of penetrable spherical sinks[END_REF]. Traytak employed the GMSV in the form of expansions with respect to irreducible Cartesian tensors [START_REF] Traytak | Theory of recondensation of N drops[END_REF][START_REF] Traytak | Methods for solution of the boundary value problems in domains with disconnected boundary[END_REF][START_REF] Traytak | The diffusive interaction in diffusion-limited reactions: the steady-state case[END_REF][START_REF] Traytak | On the irreducible tensors method in the theory of diffusive interaction between particles[END_REF]. Traytak and Tachiya investigated diffusive interaction between two spherical sinks in an electric field by means of the GMSV [START_REF] Traytak | Diffusion-controlled reactions in an electric field: Effects of an external boundary and competition between sinks[END_REF]. Tsao, Strieder et al. and Traytak et al. used the GMSV to calculate rigorously the electric field effects and to study reactions on two different spherical sinks and on spherical source and sink [START_REF] Tsao | Competitive diffusion into two reactive spheres of different reactivity and size[END_REF][START_REF] Mcdonald | Diffusion and reaction for a spherical source and sink[END_REF][START_REF] Mcdonald | Competitive interaction between two different spherical sinks[END_REF][START_REF] Seki | Rigorous calculation of electric field effects on the free energy change of the electron transfer reaction[END_REF][START_REF] Traytak | Competition effects in diffusion-controlled bulk reactions between ions[END_REF][START_REF] Traytak | On the time-dependent diffusive interaction between stationary sinks[END_REF][START_REF] Strieder | Series reactions A → B → C on successive spheres[END_REF]. A more general form of the GMSV was elaborated to compute the steady-state reaction rate for an irreversible bulk diffusion-influenced chemical reaction between a mobile point-like particle and static finite three-dimensional configurations of spherical active particles [START_REF] Traytak | Methods for solution of the boundary value problems in domains with disconnected boundary[END_REF][START_REF] Traytak | Boundary-value problems for the diffusion equation in domains with disconnected boundary[END_REF][START_REF] Galanti | Theory of diffusion-influenced reactions in complex geometries[END_REF][START_REF] Galanti | Reaction rate of a composite core-shell nanoreactor with multiple nanocatalysts[END_REF]. Since diffusion-reaction processes among spherical sinks is a long-standing problem, many other theoretical methods (such as variational estimates and perturbative analysis) have been employed [START_REF] Samson | Exact solution for the diffusion controlled rate into a pair of reacting sinks[END_REF][START_REF] Kayser | Diffusion in a Medium with a Random Distribution of Static Traps[END_REF][START_REF] Kayser | Reaction diffusion in a medium containing a random distribution of nonoverlapping traps[END_REF][START_REF] Berezhkovskii | Mutual influence of traps on the death of a Brownian Particle[END_REF][START_REF] Berezhkovskii | Diffusion-limited reactions with correlated traps[END_REF][START_REF] Berezhkovskii | Trap correlation influence on diffusion-limited process rate[END_REF][START_REF] Berezhkovskii | Brownian-particle trapping by clusters of traps[END_REF][START_REF] Yu | Role of trap clustering in the trapping kinetics[END_REF][START_REF] Yu | Effect of polydispersity on Brownianparticle trapping by clusters of traps[END_REF][START_REF] Oshanin | Kinetic description of diffusion-limited reactions in random catalytic media[END_REF][START_REF] Yang | Kinetic theory of bimolecular reactions in liquid. I. Steady-state fluorescence quenching kinetics[END_REF][START_REF] Kipriyanov | Many-particle derivation of binarykinetic equation of reaction A + B → B in liquid solutions[END_REF][START_REF] Barzykin | Kinetics of diffusion-assisted reactions in microheterogeneous systems[END_REF][START_REF] Torquato | Random Heterogeneous Materials: Microstructure and Macroscopic Properties[END_REF][START_REF] Nguyen | A Spectral Approach to Survival Probability in Porous Media[END_REF].

In a nutshell, the GMSV for many spheres consists in representing the solution of a boundary value problem as a linear combination of partial solutions written in local spherical coordinates of each sphere. The coefficients in front of the underlying solid harmonics are fixed to respect boundary conditions by using so-called addition theorems for solid harmonics to switch between local spherical coordinates. As in other numerical methods, the original PDE is reduced to an ISLAE that in general has to be solved numerically. However, the natural choice of the solid harmonics in local spherical coordinates as basis functions preserves the intrinsic symmetries of the domain and provides superior computational efficiency. In particular, the resulting ISLEA can be truncated at smaller sizes, yielding faster and more accurate solutions [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF]. Moreover, there is no need for meshing the domain: once the coefficients in front of partial solutions are computed, the concentration field can be easily and very rapidly evaluated at any point due to the explicit analytical dependence on the coordinates. For exterior problems, this method does not require imposing a distant artificial outer boundary that is needed for many other numerical methods (such as FEM) to deal with a bounded domain. It is important to note that the GMSV is not limited to spherical shapes and can be applied to other canonical domains such as spheroids, cylinders, cones, etc. [START_REF] Traytak | Diffusion-influenced reaction rates for active "sphere-prolate spheroid" pairs and Janus dimers[END_REF].

To our knowledge, the GMSV has not been applied to compute Green functions for Laplacian boundary value problems in three-dimensional domains with disconnected spherical boundaries (however, see [START_REF] Chen | Derivation of Green's function using addition theorem[END_REF] for the planar case). The Green function is harmonic everywhere in a given domain except for a fixed singularity point, and satisfies imposed homogeneous boundary conditions. The corresponding boundary value problems are known to be well-posed in simply-connected three-dimensional domains bounded by piecewise smooth boundaries (i.e., as in our setting) [START_REF] Courant | Methods of mathematical physics[END_REF][START_REF] Yu | Partial differential equation I: Foundations of the classical theory[END_REF][START_REF] Melnikov | Green's functions: Construction and applications[END_REF][START_REF] Bogolyubov | Green's function of the Laplace operator[END_REF]. We compute the Green functions for both exterior and interior domains and for all most common boundary conditions: Dirichlet, Neumann, Robin, and conjugate one (also known as the fourth boundary condition, transmission condition and exchange condition). We describe all the steps of the method, from analytical derivations to numerical implementations. We deduce the semi-analytical formula for the Green function, which is the key ingredient to solve general boundary value problems for Laplace and Poisson equations and to determine various characteristics of stationary diffusion such as reaction rate, escape probability, harmonic measure, residence time, and mean first passage time, to name but a few. An implementation of this method as a Matlab package is released and made freely accessible.

The paper is organized as follows. Section 2 presents the main results and their derivation, for both interior and exterior boundary value problems. With increasing complexity, we treat the exterior Dirichlet problem (Sec. 2.1), the exterior Robin problem (Sec. 2.2), the interior Robin problem (Sec. 2.3), and the conjugate (or exchange) problem (Sec. 2.4). To illustrate the general scheme, we summarize in Sec. [START_REF] Weiss | Overview of theoretical models for reaction rates[END_REF] several examples for which the solution is fully explicit. In Sec. 4, we discuss various applications of the derived semi-analytical formula for the Green function. Section 5 is devoted to a practical implementation of the proposed method and some numerical results. Section 6 concludes the paper, whereas some technical points are moved to Appendix A.

Semi-analytical solution

In this section, we present the detailed derivation of the Green function for the exterior Dirichlet problem (Sec. 2.1), the exterior Robin problem (Sec. 2.2), the interior Robin problem (Sec. 2.3), and the conjugate problem (Sec. 2.4).

For exterior problems, we consider an unbounded domain Ω -outside N non-overlapping balls

Ω i = {x ∈ R 3 : x -x i < R i } of radii R i , centered at x i , with i = 1, N (see Fig. 1(a)): Ω -:= R 3 \ N i=1 Ω i , Ω i := Ω i ∪ ∂Ω i , (1)
where ∂Ω i is the surface of the i-th ball, and • is the Euclidean distance. The nonoverlapping condition reads

Ω i ∩ Ω j = ∅ (i = j). (2)
The N -connected boundary of the domain Ω -is

∂Ω -= N i=1 ∂Ω i , (3)
i.e., partial surfaces ∂Ω i are the connected components of the boundary of the exterior domain Ω -. In the literature, domains like Ω -are called "periphractic domains" [START_REF] Maxwell | Treatise on electricity and magnetism[END_REF], "perforated domains" [START_REF] Marchenko | Homogenization of partial differential equations[END_REF][START_REF] Höfer | The Method of Reflections, Homogenization and Screening for Poisson and Stokes Equations in Perforated Domains[END_REF] and "domains with disconnected boundary" [START_REF] Traytak | Methods for solution of the boundary value problems in domains with disconnected boundary[END_REF].

For interior problems, we consider that N formerly introduced non-overlapping balls Ω i are englobed by a larger spherical domain Ω 0 = {x ∈ R 3 : xx 0 < R 0 } of radius R 0 , centered at x 0 . The interior boundary value problems are posed in a bounded interior domain Ω + (see Fig. 1(b))

Ω + := Ω 0 \ N i=1 Ω i , Ω i ⊂ Ω 0 . (4)
The (N + 1)-connected boundary of the domain Ω + is

∂Ω + = N i=0 ∂Ω i , (5)
that includes the outer boundary ∂Ω 0 . The Cartesian coordinates x 1 , x 2 , x 3 of the centers of these balls are given in some fixed global coordinate system. In turn, a local spherical coordinate system, (r i , θ i , φ i ), is associated with each ball. The Green function G(x, y) is computed at any pair of points x and y of Ω -. (b) Illustration of a bounded interior

x 1 R 1 1 x 2 R 2 2 x 3 R 3 3 x y - z x y L 12 L 13 L 23 (a) x 1 R 1 1 x 2 R 2 2 x 3 R 3 3 x y + z x y L 12 L 13 L 23 x 0 R 0 0 (b)
domain Ω + = Ω 0 \ N i=1
Ω i with three balls (N = 3).

Exterior Dirichlet problem

We first consider a general exterior Dirichlet boundary value problem for the Poisson equation in the unbounded domain Ω -⊂ R 3 :

-∇ 2 u = F (x ∈ Ω -), ( 6a
)
u| ∂Ωi = f i (i = 1, N ), ( 6b
)
u| x →∞ → 0, (6c)
where

F ∈ L 2 (Ω -) ∩ C 1 (Ω -
) is a given function of "sources" and f i ∈ C(∂Ω i ) are given continuous functions. The last relation (6c) is the regularity condition at infinity. This problem is well posed and has a unique classical solution [START_REF] Courant | Methods of mathematical physics[END_REF][START_REF] Yu | Partial differential equation I: Foundations of the classical theory[END_REF]. There are at least two standard ways to get the classical solution of this problem. (i) One can use the fundamental solution,

G(x, y) = 1 4π x -y , (7)
which is the Green function of the Laplace operator in R 3 \{y}:

-

∆ x G(x, y) = δ(x -y) (x ∈ R 3 \{y}), ( 8
)
where δ is the Dirac distribution, and y is a fixed point-like "source". Multiplying Eqs. (6a, 8) by G(x, y) and u(x) respectively, subtracting them, integrating over x ∈ Ω -and using the Green formula, one gets

u(y) = Ω - dx F (x) G(x, y) + ∂Ω - ds G(x, y) ∂u(x) ∂n x -u(x) ∂G(x, y) ∂n x x=s , (9)
where ∂/∂n x is the normal derivative at the surface point x, directed outward the domain Ω -. Combined with the boundary condition (6b), this representation yields a boundary integral equation on ∂u(x)/∂n x , whose solution then determines u(y) according to Eq. ( 9). The solution of this integral equation for domains with spherical and prolate spheroidal boundaries has been recently proposed by Chang et al. [START_REF] Chang | Green's Function Problem of Laplace Equation with Spherical and Prolate Spheroidal Boundaries by Using the Null-Field Boundary Integral Equation[END_REF]. Although the proposed method is conceptually close to the GMSV that we describe here, the use of translational addition theorems for solid harmonics significantly facilitates and speeds up the computation (see below). Most importantly, the dependence on the "source" point y will appear explicitly in our analysis.

(ii) The solution of the problem (6) can alternatively be written as

u(y) = Ω - dx F (x) G(x, y) + N i=1 ∂Ωi ds f i (s) - ∂G(x, y) ∂n x x=s , (10)
where G(x, y) is the Dirichlet Green function in Ω -satisfying for any y ∈ Ω -the boundary value problem

-∇ 2 x G(x, y) = δ (x -y) (x ∈ Ω -), (11a)
G| ∂Ωi = 0 (i = 1, N ), (11b)
G| x →∞ → 0. (11c)
The representation [START_REF] Gardiner | Handbook of stochastic methods[END_REF] is obtained by multiplying Eqs. (6a, 11a) by G(x, y) and u(x) respectively, subtracting them, integrating over x ∈ Ω -and using the Green formula. In spite of apparent similarity between Eqs. (9, 10), the major difference is that Eq. ( 10) is an explicit solution in terms of yet unknown Green function G(x, y), whereas Eq. ( 9) is an integral equation on ∂u(x)/∂n x involving the known fundamental solution G(x, y).

We recall that the Green function G(x, y) can be physically interpreted as the electric potential at x created by a charge at y with grounded balls [START_REF] Duffy | Green's functions with applications[END_REF].

To compute the Green function, one can represent it as

G(x, y) = G(x, y) -g(x; y), (12)
with an auxiliary function g(x; y) satisfying for any point y ∈ Ω -:

-∇ 2 x g(x; y) = 0 (x ∈ Ω -), (13a)
g(x; y)| x∈∂Ωi = G(x, y)| x∈∂Ωi (i = 1, N ), (13b)
g(x; y)| x →∞ → 0. ( 13c
)
In other words, one can separate the universal singular part G(x, y) (yielding the Dirac δ distribution) and the remaining regular part g(x, y) (ensuring the boundary conditions).

In spite of the known reciprocity,

G(x, y) = G(y, x),
we will treat the point y as a fixed parameter. In the remaining part of this subsection, we focus on the particular problem [START_REF] Miller | Symmetry and separation of variables[END_REF], bearing in mind that its solution gives access via Eq. ( 10) to a solution of any exterior Dirichlet problem [START_REF] Samarskii | The Finite Difference Methodology[END_REF]. 6

We search for the solution of Eqs. [START_REF] Miller | Symmetry and separation of variables[END_REF] in the form of superposition

g(x; y) = N i=1 g i (r i , θ i , φ i ; y), (14)
where g i is the partial solution in the local spherical coordinates of the ball Ω i , with (r i , θ i , φ i ) being the local spherical coordinates of point x, i.e., the spherical coordinates of xx i . The above expression follows immediately from the representation (9) of the function g(x; y) (with F ≡ 0) and the additivity of the Riemann integral over the disconnected boundary ∂Ω -. Each partial solution can be expanded onto a complete basis of functions {ψ - mn } outside the i-th ball:

g i (r i , θ i , φ i ; y) = ∞ n=0 n m=-n A i mn ψ - mn (r i , θ i , φ i ), (15)
where A i mn are the unknown coefficients (depending parametrically on y). Basis functions {ψ - mn } are the irregular (also called singular with respect to the origin) solid harmonics:

ψ - mn (r, θ, φ) := 1 r n+1 Y mn (θ, φ), (16)
where

Y mn (θ, φ) := P m n (cos θ) e imφ (17)
are the (non-normalized) spherical harmonics. The (non-normalized) associated Legendre functions P m n (z) of degree n (n = 0, 1, 2, . . .) and order m (m = -n, -n + 1, . . . , n -1, n) are

P m n (z) = (-1) m (1 -z 2 ) m/2 d m dz m P n (z) (m = 0, n), P -m n (z) = (-1) m (n -m)! (n + m)! P m n (z) (m = 1, n), (18)
where P n (z) are Legendre polynomials of degree n. Note that irregular solid harmonics ψ - mn and the coefficients A i mn are sometimes called the "multipoles" and the "moments of the expansion", respectively [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF]. The irregular solid harmonics are well defined in the exterior of a ball. In turn, the regular solid harmonics,

ψ + mn (r, θ, φ) := r n Y mn (θ, φ), (19)
are well defined in the interior of a ball. In order to satisfy boundary conditions, each partial solution g i , written in the local spherical coordinates of the ball Ω i , should be represented in local spherical coordinates of other balls. Such representations can be efficiently performed by so-called translational addition theorems (TATs) [START_REF] Miller | Symmetry and separation of variables[END_REF][START_REF] Epton | Multipole translation theory for the three-dimensional Laplace and Helmholtz equations[END_REF] which express a basic of solid harmonics {ψ ± mn (x i )} in local coordinates (O; x i ) via a new basis of solid harmonics {ψ ± mn (x j )} in translated coordinates x j = x i + L ij . There are three TATs that "translate" regular to regular (R→R), irregular to regular (I→R), and irregular to irregular (I→I) solid harmonics (noting that regular harmonics cannot be expanded onto irregular ones). We have thus R→R:

ψ + mn (r j , θ j , φ j ) = n l=0 l k=-l U (+j,+i) mnkl ψ + kl (r i , θ i , φ i ), (20a)
I→R:

ψ - mn (r j , θ j , φ j ) = ∞ l=0 l k=-l U (-j,+i) mnkl ψ + kl (r i , θ i , φ i ) (r i < L ij ), (20b)
I→I:

ψ - mn (r j , θ j , φ j ) = ∞ l=0 l k=-l U (-j,-i) mnkl ψ - kl (r i , θ i , φ i ) (r i > L ij ), ( 20c
)
where L ij = x j -x i is the vector connecting the centers of balls j and i, and (L ij , Θ ij , Φ ij ) are the spherical coordinates of the vector L ij :

x j = x i + L ij sin Θ ij cos Φ ij , y j = y i + L ij sin Θ ij sin Φ ij , z j = z i + L ij cos Θ ij . (21)
The coefficients U (±j,±i) mnkl are the matrix elements of the translation operator [START_REF] Epton | Multipole translation theory for the three-dimensional Laplace and Helmholtz equations[END_REF], which are also known as mixed-basis matrix elements [START_REF] Miller | Symmetry and separation of variables[END_REF]. For i = j, we have

U (+j,+i) mnkl = (n + m)! (n -l + m -k)! (k + l)! ψ + (m-k)(n-l) (L ij , Θ ij , Φ ij ), (22a)
U (-j,+i) mnkl = (-1) k+l (n + l -m + k)! (n -m)! (l + k)! ψ - (m-k)(n+l) (L ij , Θ ij , Φ ij ), (22b)
U (-j,-i) mnkl = (-1) m-n+l-k (l -k)! (n -m)! (m -n + l -k)! ψ + (m-k)(l-n) (L ij , Θ ij , Φ ij ). ( 22c
)
Note that we use the convention that ψ ± mn is zero for n < 0 or |m| > n. For instance, the elements of the matrix U (-j,-i) mnkl are zero when the inequalities l ≥ n and n + m -l ≤ k ≤ m -n + l are not fulfilled. Similarly, there is no contribution of terms with the index k such that |m -k| > n -l, i.e., the second sum in Eq. (20a) runs over k from -min(l, n -l -m) to min(l, n -l + m). In particular, the first sum in Eq. (20a) can be formally extended to +∞, as in other expansions. Since the relation (20a) involves polynomials on both sides, it is applicable for any r i . In turn, two other TATs are applicable for r i < L ij and r i > L ij respectively.

Here and throughout the text, we use the triple indices (i, m, n) and (j, k, l) to encode the elements of the involved vectors and matrices. These indices facilitate the visual interpretation of these elements, the superscript i (or j) always referring to the ball number, while the subscript mn (or kl) to the order m and the degree n of the solid harmonics (see also Sec. 5 for details on the numerical implementation). Note that the sign in front of i (or j) refers to regular (plus) and irregular (minus) solid harmonics. We will also employ the shortcut summation notations:

The unknown coefficients A i mn are fixed by the boundary condition (13b). To fulfill this condition at the surface ∂Ω i , one needs to represent g(x; y) in the local coordinates of the i-th ball. Combining Eqs. (14, 15, 20b), one gets any

i = 1, N g(x; y) = n,m A i mn r 2n+1 i + N j( =i)=1 l,k A j kl U (-j,+i) klmn ψ + mn (r i , θ i , φ i ), (24)
which is valid for x in a close vicinity of the i-th ball. We also use the Laplace expansion for the Newton's potential (see Appendix A.1) to expand the right-hand side of Eq. (13b) on the regular solid harmonics in the local coordinates of the i-th ball:

G(x, y) = n,m V i mn ψ + mn (r i , θ i , φ i ) (r i < L i ), (25)
with

V i mn = (-1) m 4π ψ - (-m)n (L i , Θ i , Φ i ), (26)
where L i = yx i and (L i , Θ i , Φ i ) are the spherical coordinates of L i . Equating Eqs. ( 24) and ( 25) at r i = R i and using (13b), one gets the equality that must be fulfilled for all points on ∂Ω i (i.e., all θ i and φ i ), implying that the coefficients in front of ψ + mn (R i , θ i , φ i ) must be identical:

A i mn R 2n+1 i + N j( =i)=1 l,k A j kl U (-j,+i) klmn = V i mn . (27)
Multiplying by R 2n+1 i and denoting

Û ij mnkl = R 2n+1 i U (-j,+i) klmn (i = j), δ mk δ nl (i = j), (28)
V i mn = R 2n+1 i V i mn , (29)
one rewrites the above relation as an ISLAE:

N j=1 l,k Û ij mnkl A j kl = V i mn (i = 1, N , n = 0, ∞, m = -n, n). (30)
Writing this system in a matrix form, one gets the vector A of coefficients A i mn by inverting the matrix Û:

A = W V, W = Û-1 . (31)
The Green function follows with the aid of Eqs. [START_REF] Crank | The Mathematics of Diffusion[END_REF][START_REF] Reddy | An Introduction to the Finite Element Method[END_REF][START_REF] Eun | Influence of neighboring reactive particles on diffusion-limited reactions[END_REF]:

G(x, y) = G(x, y) - N i=1 n,m A i mn ψ - mn (r i , θ i , φ i ). ( 32
)
In the second term, the dependence on x is captured explicitly via (r i , θ i , φ i ), whereas the dependence on y is also explicit but hidden in the coefficients A i mn via V i mn . In turn, the mixed-basis matrix U (and thus Û) depends on the positions and radii of the balls but is independent of the point y. In practice, one first inverts numerically a truncated matrix Û and then, for each point y, rapidly computes the vector V, from which the coefficients A i mn are found (see Sec. 5). As said earlier, the Green function allows one to solve any exterior Dirichlet problem [START_REF] Samarskii | The Finite Difference Methodology[END_REF]. The general representation [START_REF] Gardiner | Handbook of stochastic methods[END_REF] includes both the Green function and its normal derivative at spheres ∂Ω i , which is also known as the harmonic measure density [START_REF] Garnett | Harmonic Measure[END_REF]

ω y (s) := - ∂G(x, y) ∂n x x=s (s ∈ ∂Ω -). (33)
In Appendix A.2, we deduce the decomposition of this density onto irregular solid harmonics:

ω i y (s) := ω y (s)| ∂Ωi = 1 R i n,m (2n + 1)A i mn ψ - mn (R i , θ i , φ i ). ( 34
)

Exterior Robin problem

The above technique can be extended to finding the Green function with Robin boundary conditions

-∇ 2 x G(x, y) = δ (x -y) (x ∈ Ω -), (35a)
a i G + b i R i ∂G ∂n x x∈∂Ωi = 0 (i = 1, N ), (35b)
G| x →∞ → 0, (35c)
with a fixed source point y ∈ Ω -and nonnegative constants a i and b i such that a i +b i > 0 for each i. Using again Eq. ( 12), one gets the Robin boundary conditions for another auxiliary function g(x; y):

a i g + b i R i ∂g ∂n x x∈∂Ωi = a i G + b i R i ∂G ∂n x x∈∂Ωi . (36)
Re-writing Eqs. (24, A.6) as

g(x; y)| x∈∂Ωi = n,m ÛA i mn ψ - mn (R i , θ i , φ i ), (37a) ∂g(x; y) ∂n x x∈∂Ωi = 1 R i n,m (2n + 1)A i mn -n( ÛA) i mn ψ - mn (R i , θ i , φ i ), ( 37b
)
we represent the left-hand side of Eq. ( 36) as

n,m (a i -nb i )( ÛA) i mn + (2n + 1)b i A i mn ψ - mn (R i , θ i , φ i ). ( 38
)
Using Eq. (25, A.5), we represent the right-hand side of Eq. ( 36) as

n,m a i -nb i V i mn ψ - mn (R i , θ i , φ i ). ( 39
)
Equating Eqs. [START_REF] Venema | Reaction-diffusion on a periodic array of penetrable spherical sinks[END_REF][START_REF] Traytak | The diffusive interaction in diffusion-limited reactions: the steady-state case[END_REF], one finally gets the equalities on the coefficients A i mn in the Robin case for i = 1, N , n = 0, ∞, m = -n, n:

(2n + 1)b i A i mn + (a i -nb i ) ÛA) i mn = (a i -nb i ) V i mn . (40)
This ISLAE generalizes Eq. ( 30) to the Robin boundary condition. Representing the multiplication by a i , b i , and n in a matrix form by diagonal matrices â, b, and n, the ISLAE can be written in a matrix form as:

(2n + I) b + (â -nb ) Û A = (â -nb ) V, (41)
where I stands for the identity matrix. Inverting the matrix in front of A, one represents the vector of coefficients A i mn as

A = W V, W = (2n + I) b + (â -nb ) Û -1 (â -nb ). ( 42
)
In the Dirichlet case (b i = 0 and a i = 1), this expression is reduced to Eq. ( 31). The coefficients A i mn fully determine the Robin Green function:

G(x, y) = G(x, y) - N i=1 n,m A i mn ψ - mn (r i , θ i , φ i ). (43)
With this Green function, the solution of a general exterior Robin boundary value problem,

-∇ 2 u = F (x ∈ Ω -), (44a)
a i u + b i R i ∂u ∂n x ∂Ωi = f i (i = 1, N ), ( 44b
)
u| x →∞ → 0, (44c)
can be represented as

u(y) = Ω - dx F (x) G(x, y) + N i=1 ∂Ωi ds f i (s) ω i y (s), (45)
where

ω i y (s) = G(x, y)| x∈∂Ωi b i R i = - ∂G(x, y) a i ∂n x x∈∂Ωi (46)
is the spread harmonic measure density on the sphere ∂Ω i [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 to partially absorbing boundaries with Robin boundary condition. When both a i and b i are nonzero, two representations in Eq. ( 46) are equivalent due to Eq. (35b). In turn, one uses the first representation for the Neumann case (a i = 0) and the second representation for the Dirichlet case (b i = 0). For b i = 0, we use Eqs. [START_REF] Traytak | Methods for solution of the boundary value problems in domains with disconnected boundary[END_REF][START_REF] Jefferey | Conduction through a random suspension of spheres[END_REF] to get

ω i y (s) = 1 b i R i n,m V i mn -ÛA i mn ψ - mn (R i , θ i , φ i ), (47)
whereas Eq. ( 34) is used for the Dirichlet case (b i = 0, a i = 1).

Interior Robin problem

In many applications, a domain is limited by an outer boundary which can significantly affect the diffusion characteristics. A prominent example is the mean first passage time which is infinite for unbounded domains. In order to deal with such problems, one needs to incorporate an outer boundary, transforming the exterior problem to the interior problem in a bounded domain Ω + from Eq. ( 4), with N non-overlapping balls Ω i (i = 1, N ), englobed by a larger ball Ω 0 of radius R 0 and centered at x 0 . The Robin Green function for the interior problem in Ω + satisfies for any y ∈ Ω + :

-∇ 2 x G(x, y) = δ (x -y) (x ∈ Ω + ), (48a)
a i G + b i R i ∂G ∂n x x∈∂Ωi = 0 (i = 0, N ), (48b)
with nonnegative parameters a i and b i such that a i +b i > 0 for each i, and

a 0 +• • •+a N > 0.
The last inequality excludes the case with Neumann conditions at all boundaries, for which the Green function of an interior problem does not exist. Since the Green function is now defined in a bounded domain, there is no regularity condition (35c) at infinity. As previously, one represents the Green function as in Eq. ( 12) and then searches for an auxiliary function g(x; y) in the form

g(x; y) = N i=0 g i (r i , θ i , φ i ; y), (49)
to which a new function g 0 is added

g 0 (r 0 , θ 0 , φ 0 ; y) = n,m A 0 mn ψ + mn (r 0 , θ 0 , φ 0 ), (50)
where A 0 mn are the unknown coefficients, and (r 0 , θ 0 , φ 0 ) are the spherical coordinates of xx 0 . As this function describes the behavior inside the ball Ω 0 , one uses regular harmonics ψ + mn instead of irregular ones for other functions g i . The remaining derivation is similar to the exterior case, i.e., one needs to find the coefficients A i mn from the boundary conditions.

At the boundary of an inner ball Ω i , one re-expand ψ - mn (r j , θ j , φ j ) for j = 1, N (with j = i) as previously. In turn, one needs the R→R addition theorem (20a) to re-expand the function g 0 in the local coordinates of the i-th ball:

g 0 (r 0 , θ 0 , φ 0 ; y) = n,m A 0 mn l,k U (+0,+i) mnkl ψ + kl (r i , θ i , φ i ). ( 51
)
At the boundary ∂Ω i , one finds then

g 0 | ∂Ωi = n,m ψ + mn (R i , θ i , φ i ) l,k U (+0,+i) klmn A 0 kl , (52a)
R i ∂g 0 ∂n x ∂Ωi = - n,m n ψ + mn (R i , θ i , φ i ) l,k U (+0,+i) klmn A 0 kl . (52b)
Combining these contributions with other g i , we retrieve Eqs. [START_REF] Mattern | Rate of diffusion-controlled reactions in a random array of spherical sinks[END_REF], in which the matrix Û from Eq. ( 28) is modified by adding a new column j = 0 (with i > 0):

Û i0 mnkl = R 2n+1 i U (+0,+i) klmn . ( 53
)
As a consequence, the Robin boundary condition (35b) at each ∂Ω i implies the ISLAE [START_REF] Traytak | On the irreducible tensors method in the theory of diffusive interaction between particles[END_REF], as for the exterior problem. Here, the effect of the outer boundary is captured through the additional elements of the matrix Û in Eq. ( 53).

Moreover, one has to fulfill the Robin boundary condition (48b) at the outer boundary ∂Ω 0 . For this purpose, each g i is re-expanded by using the I→I addition theorem (20c) as

g i (r i , θ i , φ i ; y) = n,m A i mn l,k U (-i,-0) mnkl ψ - kl (r 0 , θ 0 , φ 0 ) (r 0 > L 0i ), (54)
from which

(g i )| ∂Ω0 = n,m ψ - mn (R 0 , θ 0 , φ 0 ) l,k U (-i,-0) klmn A i kl , (55a)
R 0 ∂g i ∂n x ∂Ω0 = - n,m (n + 1)ψ - mn (R 0 , θ 0 , φ 0 ) l,k U (-i,-0) klmn A i kl , (55b)
where

L 0i = x i -x 0 .
In addition, Eq. ( 50) implies

(g 0 )| ∂Ω0 = n,m R 2n+1 0 A 0 mn ψ - mn (R 0 , θ 0 , φ 0 ), (56a)
R 0 ∂g 0 ∂n x ∂Ω0 = n,m nR 2n+1 0 A 0 mn ψ - mn (R 0 , θ 0 , φ 0 ). ( 56b
)
Combining these relations, one finds

g| ∂Ω0 = n,m ψ + mn (R 0 , θ 0 , φ 0 ) ÛA 0 mn , (57a)
R 0 ∂g ∂n x ∂Ω0 = n,m ψ + mn (R 0 , θ 0 , φ 0 ) (2n + 1)A 0 mn -(n + 1)( ÛA) 0 mn , (57b)
where the matrix Û is modified by adding a new row at i = 0 (with j > 0) as

Û 0j mnkl = R -(2n+1) 0 U (-j,-0) klmn . ( 58
)
On the other hand, using again the Laplace expansion for the Newton's potential (see Appendix A.1), one can write the fundamental solution G(x, y) as

G(x, y) = n,m Ṽ 0 mn ψ - mn (r 0 , θ 0 , φ 0 ) (r 0 > L 0 ), ( 59
) with Ṽ 0 mn = (-1) m 4π ψ + (-m)n (L 0 , Θ 0 , Φ 0 ), (60) 13
where L 0 = yx 0 . One also gets

R 0 ∂G ∂n x ∂Ω0 = - n,m (n + 1) Ṽ 0 mn ψ - mn (R 0 , θ 0 , φ 0 ). ( 61
)
Combining the above expressions for g and G and their normal derivatives according to the Robin boundary condition (48b) at the outer boundary ∂Ω 0 , one gets the ISLAE for all n = 0, ∞ and m = -n, n:

(2n + 1)b 0 A 0 mn + (a 0 -(n + 1)b 0 )( ÛA) 0 mn = (a 0 -(n + 1)b 0 ) V 0 mn , (62)
where the vector V is modified at i = 0 as

V 0 mn = Ṽ 0 mn R -2n-1 0 . ( 63
)
Combining Eqs. [START_REF] Traytak | On the irreducible tensors method in the theory of diffusive interaction between particles[END_REF][START_REF] Torquato | Random Heterogeneous Materials: Microstructure and Macroscopic Properties[END_REF], one gets a complete ISLEA that fully determines all the coefficients A i mn . As previously, the solution can be written in a matrix form as

A = W V, W = (2n + I) b + (â -n b) Û -1 (â -n b), (64)
where the new matrix n includes the change of n to n + 1 in front of b 0 in Eq. ( 62) for the outer boundary:

n ij mnkl = δ ij δ mk δ nl (n + δ i0 ). (65)
The Green function reads

G(x, y) = G(x, y) - N i=1 n,m A i mn ψ - mn (r i , θ i , φ i ) - n,m A 0 mn ψ + mn (r 0 , θ 0 , φ 0 ). (66)
Note that when R 0 goes to infinity, the elements V 0 mn , as well as the nondiagonal elements of the matrix Û corresponding to A 0 mn , vanish, so that A 0 mn = 0 and one retrieves the solution for the exterior problem.

With this Green function, the solution of a general interior Robin boundary value problem,

-∇ 2 u = F (x ∈ Ω + ), (67a)
a i u + b i R i ∂u ∂n x ∂Ωi = f i (i = 0, N ), (67b)
can be represented as

u(y) = Ω + dx F (x) G(x, y) + N i=0 ∂Ωi ds f i (s) ω i y (s), (68)
where the spread harmonic measure density ω i y (s) is expressed through G(x, y) by Eq. ( 46). For i > 0, ω i y (s) is given by Eq. ( 47) for b i > 0, and by Eq. ( 34) for b i = 0. In turn, for i = 0, one finds

ω 0 y (s) =          1 b 0 R 0 n,m V 0 mn -ÛA 0 mn ψ + mn (R 0 , θ 0 , φ 0 ) (b 0 > 0), 1 a 0 R 0 n,m (2n + 1)A 0 mn ψ + mn (R 0 , θ 0 , φ 0 ) (b 0 = 0). (69)

Conjugate problems

In many biological applications, the diffusive transport occurs in heterogeneous media, with distinct diffusion coefficients in different regions. A pack of balls is a basic model of a tissue that is formed by individual cells located in the extracellular space [START_REF] Grebenkov | Pulsed-gradient spin-echo monitoring of restricted diffusion in multilayered structures[END_REF]. A diffusing molecule can cross cell membranes and move from a cell to the extracellular space and back. Such diffusion processes are often described with conjugate boundary conditions on the surface between any two adjacent "compartments" of the medium (also known as the fourth boundary condition, transmission condition, and exchange condition). When the surface is fully permeable, the concentration u of diffusing molecules obeys two conditions: (i) the continuity of the concentration,

u| ∂Ω -= u| ∂Ω + ; (70)
and (ii) the continuity of the diffusive flux at the surface,

-D -∂u ∂ nx ∂Ω - = D + ∂u ∂ nx ∂Ω + , (71)
where D ± are diffusion coefficients on both sides of the surface (denoted by ∂Ω ± ). Note that the normal derivatives ∂/∂ nx on both sides are directed outwards the corresponding compartment and thus opposite. When the membrane presents some "resistance" to exchange between compartments, the first condition is replaced by

-D -∂u ∂ nx ∂Ω - = κ u| ∂Ω --u| ∂Ω + , (72)
which states that the diffusive flux is proportional to the drop of concentrations on both sides. Here κ ≥ 0 is the permeability of the surface quantifying how difficult is to cross it: the limit κ = 0 corresponds to a fully impermeable boundary (in which case one recovers two uncoupled Neumann conditions on both sides), whereas the limit κ → ∞ describes the former situation of a fully permeable surface (in which case Eq. ( 72) is reduced to Eq. ( 70)). Only in the case of a fully impermeable surface, one can treat the boundary value problem separately in two compartments, in particular, one can use solutions from previous subsections to describe separately intracellular and extracellular diffusions. In contrast, whenever κ > 0, the two problems are coupled and should thus be treated simultaneously. As a consequence, the conjugate problems are more difficult to solve. Moreover, even a general formulation of conjugate problems is more challenging because one can imagine a large compartment (e.g., the extracellular space) filled with smaller compartments, each of them is filled with even smaller compartments, and so on (like Russian nested dolls). Although the GMSV can still be applied to such complicated cases (when all compartments are spherical), we do not consider this general setting. For illustration purposes, we limit ourselves to the practically relevant situation of N non-overlapping balls Ω i and the extracellular space Ω -= R 3 \ ∪ N i=1 Ω i . We search for the Green function G(x, y) that satisfies general conjugate boundary conditions

a i G + b i R i ∂G ∂n x x∈∂Ω - i = āi G + bi R i ∂G ∂n x x∈∂Ω + i , (73a)
c i G + d i R i ∂G ∂n x x∈∂Ω - i = ci G + di R i ∂G ∂n x x∈∂Ω + i , (73b)
where the parameters a i , b i , c i , d i characterize the exterior compartment Ω -(near the surface ∂Ω i ), while the parameters āi , bi , ci , di characterize the interior spherical compartment Ω i . As one needs to relate the Green function in the exterior compartment to that in the interior compartment, there are two conjugate relations at ∂Ω i , in contrast to former Robin boundary conditions with a single relation. These two relations should be linearly independent, i.e. one relation should not be reduced to the other (e.g., as Eqs. [START_REF] Yu | Partial differential equation I: Foundations of the classical theory[END_REF][START_REF] Melnikov | Green's functions: Construction and applications[END_REF]).

For convenience, we denote the restrictions of G(x, y) to Ω -and to Ω i as G -and G +i , respectively. We consider separately two cases: y ∈ Ω -and y ∈ Ω i .

(i) When y ∈ Ω -, each function G +i (x, y) satisfies the Laplace equation in Ω i , ∇ 2

x G +i (x, y) = 0, and it is thus naturally decomposed onto the regular solid harmonics in the local spherical coordinates of the ball Ω i :

G +i (x, y) = n,m Āi mn ψ + mn (r i , θ i , φ i ), (74)
with unknown coefficients Āi mn . In turn, the function G -(x, y) can be represented as G(x, y) -g(x; y), with an auxiliary function g satisfying the Laplace equation in Ω -, and G given by Eq. ( 7). The conjugate boundary conditions [START_REF] Bogolyubov | Green's function of the Laplace operator[END_REF] read

a i g + b i R i ∂g ∂n x + āi G +i + bi R i ∂G +i ∂n x x∈∂Ωi = a i G + b i R i ∂G ∂n x x∈∂Ωi , (75a)
c i g + d i R i ∂g ∂n x + ci G +i + di R i ∂G +i ∂n x x∈∂Ωi = c i G + d i R i ∂G ∂n x x∈∂Ωi , (75b)
where ∂Ω ± i were replaced by ∂Ω i , as the appropriate side of the surface is now clear from notations. The function g(x; y) is represented again as the sum (14) of partial solutions. With the aid of addition theorems, one can express g in the local coordinates of the ball Ω i , whereas the left-hand side of Eqs. ( 75) is an explicit function, which can be decomposed over the regular solid harmonics. Repeating the steps of Sec. 2.2, we get the ISLAE for i = 1, N , n = 0, ∞, m = -n, n:

(2n + 1)b i A i mn + (a i -nb i ) ÛA) i mn + R 2n+1 i (ā i + n bi ) Āi mn = (a i -nb i ) V i mn , (76a) (2n + 1)d i A i mn + (c i -nd i ) ÛA) i mn + R 2n+1 i (c i + n di ) Āi mn = (c i -nd i ) V i mn , (76b)
where the elements of Û and V were defined by Eqs. [START_REF] Gumerov | Computation of scattering from clusters of spheres using the fast multipole method[END_REF][START_REF] Kushch | Multipole expansion method in micromechanics of composites[END_REF]. These relations generalize Eqs. ( 40) by the inclusion of the terms with Āi mn that account for coupling between exterior and interior compartments. Although the number of unknowns is doubled (A i mn and Āi mn ), the number of equations is also doubled. Writing these equations in a matrix form, one can solve the truncated system to determine the unknown coefficients and thus the Green function.

(ii) When y ∈ Ω k for some k, each function G +i (x, y) with i = k satisfies the Laplace equation in Ω i and can thus be searched in the form [START_REF] Chang | Green's Function Problem of Laplace Equation with Spherical and Prolate Spheroidal Boundaries by Using the Null-Field Boundary Integral Equation[END_REF]. Moreover, G -satisfies the Laplace equation in Ω -so that one can set G -= g and search it as the sum (14) of partial solutions. In turn, the function G +k (x, y) can be represented as G(x, y)-g +k (x; y), with

g +k (x; y) = n,m Āk mn ψ + mn (r k , θ k , φ k ). ( 77
)
The conjugate boundary conditions [START_REF] Bogolyubov | Green's function of the Laplace operator[END_REF] read

a k g + b k R k ∂g ∂n x + āk g +k + bk R k ∂g +k ∂n x x∈∂Ω k = a k G + b k R k ∂G ∂n x x∈∂Ω k , (78a)
c k g + d k R k ∂g ∂n x + ck g +k + dk R k ∂g +k ∂n x x∈∂Ω k = c k G + d k R k ∂G ∂n x x∈∂Ω k (78b)
for i = k, and

a i g + b i R i ∂g ∂n x -āi G +i -bi R i ∂G +i ∂n x x∈∂Ωi = 0, ( 79a
)
c i g + d i R i ∂g ∂n x -ci G +i -di R i ∂G +i ∂n x x∈∂Ωi = 0 ( 79b
)
for i = k. One can repeat again the steps of Sec. 2.2 to derive linear equations on the unknown coefficients. The only difference is that one needs to employ another representation of the fundamental solution inside Ω k (similar to Eq. ( 59) for the interior problem):

G(x, y) = n,m Ṽ k mn ψ - mn (r k , θ k , φ k ) (r k > L k ), (80)
where

Ṽ k mn = (-1) m 4π ψ + (-m)n (L k , Θ k , Φ k ), (81)
with L k = yx k . We get thus

(2n + 1)b k A k mn + (a k -nb k ) ÛA) k mn + R 2n+1 k (ā k + n bk ) Āk mn = (a k -(n + 1)b k ) V k mn , (82a)
(2n + 1)d k A k mn + (c k -nd k ) ÛA) k mn + R 2n+1 k (c k + n dk ) Āk mn = (c k -(n + 1)d k ) V k mn (82b)
for i = k, and

(2n + 1)b i A i mn + (a i -nb i ) ÛA) i mn -R 2n+1 i (ā i + n bi ) Āi mn = 0, (83a)
(2n + 1)d i A i mn + (c i -nd i ) ÛA) i mn -R 2n+1 i (c i + n di ) Āi mn = 0 (83b) for i = k, with n = 0, ∞ and m = -n, n.
The above computation can be straightforwardly extended to the case when the extracellular space is bounded by a large ball Ω 0 . In this case, the analysis of the exterior part (i.e., the evaluation of the function g) should follow Sec. 2.3 instead of Sec. 2.2. Finally, one can consider more general problems, in which some balls are partially absorbing sinks or impermeable obstacles (with Robin or Neumann boundary conditions), whereas the other ball allow for interior diffusion (with conjugate boundary conditions). One just combines the corresponding conjugate conditions with Robin boundary conditions. Moreover, it is worth noting that the conjugate problem naturally includes the Robin boundary value problem as a particular case. In fact, setting G +i ≡ 0 in Eq. (73a) and removing Eq. (73b), one recovers the Robin boundary condition (35b).

Basic examples

Interior problem for two concentric spheres

As an example of an interior Robin problem, we determine the Green function in a bounded domain Ω + between two concentric spheres, centered at x 0 = x 1 = 0 and of radii R 0 > R 1 . Although this problem could be solived via the spectral decomposition over the known Laplacian eigenfunctions, our derivation yields a more explicit formula and serves as an illustration for the proposed method.

One can check that

Û 10 mnkl = R 2n+1 1 δ ln δ km , Û 01 mnkl = R -2n-1 0 δ ln δ km , (84)
i.e., the matrix Û is formed by four diagonal matrices. This structure is preserved by multiplication by diagonal matrices â, b, n, and n in the matrix relation [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF]. The resulting matrix W has the block three-diagonal structure:

W 11 mnkl = δ nl δ km w 11 n , W 10 mnkl = -δ nl δ km w 10 n R 2n+1 1 , W 01 mnkl = -δ nl δ km w 01 n R -2n-1 0 , W 00 mnkl = δ nl δ km w 00 n , (85)
where

w 11 n = w n (a 1 -nb 1 )(a 0 + nb 0 ), w 10 n = w n (a 1 -nb 1 ) 2 , w 01 n = w n (a 0 -(n + 1)b 0 ) 2 , w 00 n = w n (a 1 + (n + 1)b 1 )(a 0 -(n + 1)b 0 ), (86)
and

w n = 1 (a 1 + (n + 1)b 1 )(a 0 + nb 0 ) -(a 1 -nb 1 )(a 0 -(n + 1)b 0 )(R 1 /R 0 ) 2n+1 . ( 87
)
Since

L 1 = y -x 1 = y = y -x 0 = L 0 , one has V 1 mn = R 2n+1 1 (-1) m 4π ψ - (-m)n (L 1 , Θ 1 , Φ 1 ) = (-1) m 4π R 2n+1 1 r n+1 y Y (-m)n (θ y , φ y ), (88a)
V 2 mn = R -2n-1 0 (-1) m 4π ψ + (-m)n (L 0 , Θ 0 , Φ 0 ) = (-1) m 4π r n y R 2n+1 0 Y (-m)n (θ y , φ y ), (88b)
where (r y , θ y , φ y ) are the spherical coordinates of y. One gets thus

A 1 mn = (-1) m 4π Y (-m)n (θ y , φ y ) R 2n+1 1 r n+1 y w 11 n -w 10 n (r y /R 0 ) 2n+1 , (89a)
A 2 mn = (-1) m 4π Y (-m)n (θ y , φ y ) r n y R 2n+1 0 w 00 n -w 01 n (R 1 /r y ) 2n+1 , (89b)
from which

G(x, y) = G(x, y) - 1 4π ∞ n=0 P n (x • y) x y R 2n+1 1 ( x y ) n+1 w 11 n -w 10 n ( y /R 0 ) 2n+1 + ( x y ) n R 2n+1 0 w 00 n -w 01 n (R 1 / y ) 2n+1 , (90)
where the sum over m was calculated by using the classical addition theorem for two unit vectors e and e :

n m=-n (-1) m Y mn (θ, φ)Y (-m)n (θ , φ ) = P n (e • e ). ( 91
)
The spread harmonic measure density ω y follows from its definition [START_REF] Traytak | Boundary-value problems for the diffusion equation in domains with disconnected boundary[END_REF]:

ω 1 y (s) = - ∂G a 1 ∂n x x=s∈∂Ω1 = 1 4πa 1 R 2 1 ∞ n=0 P n (s • y) R 1 y (92) × R n+1 1 y n+1 n + (n + 1)w 11 n + n(R 1 /R 0 ) 2n+1 w 01 n - y n R n+1 1 R 2n+1 0 nw 00 n + (n + 1)w 10 n ,
from which the absorption probability (see also Sec. 4) reads

p 1 (y) = a 1 ∂Ω1 ds ω 1 y (s) = w 11 0 R 1 / y -w 10 0 R 1 /R 0 . (93)
Similarly, we get

ω 0 y (s) = - ∂G a 0 ∂n x x=s∈∂Ω0 = 1 4πa 0 R 2 0 ∞ n=0 P n (s • y) R 0 y 1 R n 0 × y n (n + 1) + (n + 1)w 10 n (R 1 /R 0 ) 2n+1 + nw 00 n - R 2n+1 1 y n+1 (n + 1)w 11 n + nw 01 n , (94)
from which p 0 (y) = a 0

∂Ω0 ds ω 0 y (s) = 1 + w 10 0 R 1 /R 0 -w 11 0 R 1 / y . (95)
Note that p 1 + p 0 = 1 as expected.

We are not aware of earlier derivations of the Robin Green function and the spread harmonic measure for two concentric spheres in such simple forms. The Green functions in four limiting cases (i.e., Dirichlet-Dirichlet, Dirichlet-Neumann, Neumann-Dirichlet, and Neumann-Neumann conditions on the inner and outer spheres) was provided in [START_REF] Chang | Green's Function Problem of Laplace Equation with Spherical and Prolate Spheroidal Boundaries by Using the Null-Field Boundary Integral Equation[END_REF]. These solutions can be easily deduced from our general formula [START_REF] Grebenkov | Residence times and other functionals of reflected Brownian motion[END_REF]. For instance, in the Dirichlet case (b 1 = b 0 = 0), one has w ij n = w n , and the formula (90) reads

G(x, y) = G(x, y) - ∞ n=0 P n (x • y) x y (96) × R 2n+1 1 (R 2n+1 0 -x 2n+1 ) + y 2n+1 ( x 2n+1 -R 2n+1 1 ) 4π x n+1 y n+1 (R 2n+1 0 -R 2n+1 1 )
.

Note also that an explicit form of the Dirichlet Green function for two nonconcentric spheres was derived by using bispherical coordinates in [START_REF] Chen | A Green's Function for the Domain Bounded by Nonconcentric Spheres[END_REF].

3.2.

Interior Robin problem for one sphere

In the limit R 1 → 0, Eq. ( 90) is reduced to

G(x, y) = G(x, y) - 1 4π ∞ n=0 P n (x • y) x y x n y n R 2n+1 0 a 0 -(n + 1)b 0 a 0 + nb 0 , (97)
i.e., we get the Green function for the interior Robin problem in a ball of radius R 0 . In the Dirichlet case, setting b 2 = 0 and using the identity

1 1 -2qz + z 2 = ∞ n=0 P n (q) z n , (98)
one retrieves the classical result

G(x, y) = 1 4π x -y - R 0 / y 4π x -yR 2 0 / y 2 , (99)
which is usually deduced by the image method [START_REF] Courant | Methods of mathematical physics[END_REF]. From the identity (98), one also gets

qz -z 2 (1 -2qz + z 2 ) 3/2 = ∞ n=0 nP n (q) z n , (100)
that helps to deduce the harmonic measure density

ω y (s) = 1 4πR 0 R 2 0 -y 2 y -s 3 . ( 101
)
Note that the Green function does not exist for the interior Neumann problem. This can be seen directly from Eq. (97) which diverges as a 0 → 0.

Exterior Robin problem for one sphere

In the limit R 0 → ∞, Eq. ( 90) is reduced to

G(x, y) = G(x, y) - 1 4π ∞ n=0 P n (x • y) x y R 2n+1 1 x n+1 y n+1 a 1 -nb 1 a 1 + (n + 1)b 1 , (102)
i.e., we get the Green function for the exterior Robin problem outside the ball of radius R 1 . The integral of the spread harmonic measure over the sphere ∂Ω 1 yields the absorption probability on the partially absorbing sink of radius R 1 :

p 1 (y) = a 1 a 1 + b 1 R 1 y . ( 103
)
In the Dirichlet case, the sum in Eq. ( 102) is again reduced to (99) (in which R 0 is replaced by R 1 ), whereas the (non-normalized) harmonic measure density becomes

ω y (s) = 1 4πR 1 y 2 -R 2 1 y -s 3 . ( 104
)
Its integral over the sphere yields the classical result

p 1 (y) = R 1 y . (105)
In the Neumann case, the sum in Eq. (102) becomes

G(x, y) = G(x, y) + 1 4π ∞ n=0 P n (x • y) x y R 2n+1 1 x n+1 y n+1 1 - 1 n + 1 . ( 106
)
The first sum was already computed in the Dirichlet case, whereas the second sum can be evaluated by taking the integral of Eq. ( 98) from 0 to t, yielding

G(x, y) = 2 4π x -y - R 1 / y 4π x -yR 2 1 / y 2 (107) + 1 4πR 1 ln R 2 1 -(x • y) + y x -yR 2 1 / y 2 x y -(x • y)
.

Note that the Dirichlet Green function for exterior and interior problems for a prolate spheroid was recently analyzed in [START_REF] Xue | Green's function and image system for the Laplace operator in the prolate spheroidal geometry[END_REF].

Applications

The knowledge of the Green function G(x, y) provides the solution of any boundary value problem associated to the Laplace or Poisson equation. In this section, we just mention several quantities that often appear in various applications and can be directly deduced by using our solution.

Hitting and splitting probabilities

As mentioned earlier, the normal derivative of the Green function yields the harmonic measure density ω i y (s), which characterizes the likelihood for Brownian motion started from y to arrive at the absorbing boundary for the first time in a vicinity of the boundary point s ∈ Ω i [START_REF] Garnett | Harmonic Measure[END_REF]. Integrating Eq. (34) over the sphere ∂Ω i , one gets the probability of the first arrival onto the ball Ω i :

p i (y) = ∂Ωi ds ω i y (s) ∂Ωi = 4πA i 00 , (108)
where only the rotation-invariant term with n = m = 0 survived. This is also known as the hitting probability or the splitting probability, i.e., the probability of arrival at the ball Ω i before arriving onto other balls or escaping at infinity. From this relation, the escape probability P ∞ (y) from a fixed starting point y to infinity is

P ∞ (y) = 1 - N i=1 p i (y) = 1 -4π N i=1 A i 00 . ( 109
)
In other words, P ∞ is the probability that the particle does not hit any absorbing sink. This is a nontrivial quantity in three dimensions because of the transient character of Brownian motion (in contrast, P ∞ is always zero in two dimensions). We recall that the dependence of P ∞ on y enters through the coefficients A i 00 that are expressed as linear combinations of V i mn . When the balls are only partially absorbing with Robin boundary conditions, the first arrival onto the ball does not necessarily imply absorption or chemical reaction, as the particle can be reflected. The absorption can thus be realized after numerous returns to the ball. The probability density of such absorption events is called the spread harmonic measure density ω i y (s) and given by Eq. ( 47). Integrating again this density over the boundary ∂Ω i , one gets the probability of absorption on the partially absorbing sphere ∂Ω i as

p i (y) := ∂Ωi ds - ∂G(x, y) ∂n x x=s∈∂Ωi = a i ∂Ωi ds ω i y (s) = 4πA i 00 (a i > 0), ( 110
)
the last equality coming from Eq. [START_REF] Traytak | On the irreducible tensors method in the theory of diffusive interaction between particles[END_REF]. Note that p i (y) = 0 if a i = 0 that corresponds to the Neumann boundary condition. The escape probability is still given by Eq. ( 109). Finally, when the balls Ω i are englobed by a larger ball Ω 0 , the diffusing particle cannot escape to infinity but can be absorbed by the outer boundary ∂Ω 0 . The corresponding spread harmonic measure density ω 0 y (s) is given by Eq. ( 69). Integrating this density over the sphere ∂Ω 0 , one gets

p 0 (y) := a 0 ∂Ω0 ds ω 0 y (s) = 4πR 0 × V 0 00 -ÛA 0 00 a 0 /b 0 (b 0 > 0), A 0 00 (b 0 = 0). ( 111
)
Note also that if at least one a i is nonzero, then the probabilities p i in a bounded domain satisfies

N i=0 p i (y) = 1. ( 112
)
This relation can be obtained by integrating Eq. (11a) over x ∈ Ω + , applying the Green formula and using Eq. ( 46). In probabilistic terms, it simply means that a particle released at y unavoidably arrives at some sink in a bounded domain.

Diffusive flux and reaction rate

In chemical kinetics, the escape probability P ∞ (y) from Eq. ( 109) can be interpreted as a concentration n(y) of species B diffusing from infinity towards partially absorbing sinks (species A) [START_REF] Tachiya | General method for calculating the escape probability in diffusion-controlled reactions[END_REF]. Although this particular problem was thoroughly investigated by using the GMSV in [START_REF] Galanti | Theory of diffusion-influenced reactions in complex geometries[END_REF], one can easily re-derive the former results from our more general semi-analytical solution for the Green function. In fact, the latter problem is conventionally formulated as an exterior boundary value problem

-∇ 2 n B = 0 (x ∈ Ω -), ( 113a
)
a i n B + b i R i ∂n B ∂n x ∂Ωi = 0 (i = 1, N ), ( 113b
)
n B | x →∞ → n 0 , (113c)
i.e., the field of concentration with a constant n 0 at infinity and partially absorbing sinks.

Setting n(x) = n 0 [1 -u(x)], one easily shows that this is a specific case of the general problem [START_REF] Traytak | On the time-dependent diffusive interaction between stationary sinks[END_REF] so that the solution is

u(y) = N i=1 ∂Ωi ds ω i y (s) = N i=1 p i (y). ( 114
)
As a consequence,

n B (y) = n 0 P ∞ (y), (115)
i.e., the concentration field is proportional to the escape probability.

The flux onto the sink Ω i can be computed as (see Appendix A.3)

J i := ∂Ωi ds -D ∂n B ∂n y y=s = πn 0 DR i N j=1 W ji 0000 , ( 116
)
where D is the diffusion coefficient and the matrix W is defined by Eq. ( 42). The total flux is just the sum of J i :

J := N i=1 J i = 4πn 0 D N i,j=1 W ji 0000 R i . (117)
In the case of a single spherical sink, this formula yields the classical Collins-Kimball relation [START_REF] Rice | Diffusion-limited reactions[END_REF][START_REF] Collins | Diffusion-controlled reaction rates[END_REF]

J = 4πn 0 DR 1 1 + b 1 /a 1 , (118)
which for b 1 = 0 is reduced to the famous Smoluchowski formula.

In some applications, the source of particles cannot be treated as infinitely distant. To account for a finite distance to the source, one assumes that the particles are constantly released from an outer spherical boundary ∂Ω 0 , in which case the concentration of particles, n B (x), satisfies

-∇ 2 n B = 0 (x ∈ Ω + ), ( 119a
)
a i n B + b i R i ∂n B ∂n x ∂Ωi = 0 (i = 1, N ), ( 119b
)
n B | ∂Ω0 = n 0 . ( 119c
)
The solution of this interior problem is simply

n B (x) = n 0 p 0 (x), ( 120
)
where p 0 is given by Eq. ( 111). This problem has found numerous applications in physics, electrochemistry, and biology [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 | Mathematical Basis for a General Theory of Laplacian Transport towards Irregular Interfaces[END_REF][START_REF] Gill | Modeling oxygen transport in human placental terminal villi[END_REF]. In particular, the diffusive flux can be expressed by using the Dirichlet-to-Neumann operator [START_REF] Grebenkov | Mathematical Basis for a General Theory of Laplacian Transport towards Irregular Interfaces[END_REF]. Our general solution allows one to investigate the spectral properties of this pseudo-differential operator in various configurations of spherical sinks.

Residence time and other functionals of Brownian motion

The Dirichlet Green function is related to the expectation of functionals of Brownian motion B t started from a point y, according to the formula [START_REF] Mörters | Brownian motion, Cambridge Series in Statistical and Probabilistic Mathematics[END_REF]

E y    τ 0 dt f (B t )    = Ω ± dx f (x) G(x, y), (121)
where E y is the expectation, f is a measurable function, and τ is the first passage time to the boundary of Ω ± : τ = inf{t > 0 : B t ∈ ∂Ω ± } (this is valid for both exterior and interior cases). In particular, if

f (x) = 1 D I C (x),
where D is the diffusion coefficient and I C (x) is the indicator function of a subset C ⊂ Ω ± , i.e.

I C (x) := 1 if x ∈ C, 0 if x / ∈ C,
then the functional (121) is the residence (or occupation) time in C, i.e., the time that Brownian motion spends in C until the first arrival onto the boundary ∂Ω ± , or escape at infinity [START_REF] Mörters | Brownian motion, Cambridge Series in Statistical and Probabilistic Mathematics[END_REF][START_REF] Grebenkov | Residence times and other functionals of reflected Brownian motion[END_REF]. When C is a ball Ω I of radius R I that is centered at x I and does not intersect the sinks, the addition theorem (20b) allows one to compute the residence time as (see Appendix A.4)

T (y) = 1 D Ω I dx G(x, y) = 4πR 3 I 3D 1 4πL I - N j=1 n,m A j mn ψ - mn (L Ij , Θ Ij , Φ Ij ) , (122)
where L Ij = x j -x I , (L Ij , Θ Ij , Φ Ij ) are the spherical coordinates of L Ij , and L I = yx I . Note that this result can also be extended to the case when C is an arbitrary union of non-overlapping balls.

Mean first passage time

For the interior problem, an immediate application of the semi-analytical form of the Green function is related to the mean first passage time (MFPT), T i (y), to the sink Ω i when a particle is started from y and reflected from all other sinks. The MFPT satisfies

-∇ 2 T i = 1/D (y ∈ Ω + ), (123a)
T i | ∂Ωi = 0, ( 123b
) ∂T i ∂n y ∂Ωj = 0 (j = i). (123c)
By definition of the Green function, one has

T i (y) = 1 D Ω + dx G(x, y), (124)
where the Green function satisfies the boundary conditions (48b), with a j = δ ij and b j = 1-δ ij . Note that the integral in Eq. ( 124) can be computed explicitly (see Appendix A.5). Moreover, Eq. ( 124) for the MFPT resembles Eq. ( 122) for the residence time, the main difference between two quantities lying in the boundary conditions and thus in the coefficients A i mn . More generally, one can find the MFPT to any combination of absorbing/reflecting sinks or with more general partial reflections. In addition, one can consider the spacedependent diffusion coefficient, in which case the factor 1/D(x) would remain under the integral.

Numerical aspects

Implementation

Our semi-analytical solution for the Green function in both exterior and interior problems is exact and valid for any configuration of non-overlapping balls (with or without an outer spherical boundary). An approximation is only involved at the implementation step of this solution that requires truncation of infinite-dimensional matrices, vectors and series. By setting the maximal degree n max of solid harmonics, one truncates all the series for n > n max or l > n max . We replace the triple index (i, m, n) of A i mn by a single index of a vector A of size M = N (n max + 1) 2 for exterior problems, and of size M = (N + 1)(n max + 1) 2 for interior problems (note that the size is doubled for conjugate problems). Similarly, the truncated matrix Û and the truncated vector V are of sizes M × M and M , respectively. For implementing Robin boundary conditions, one also constructs the truncated diagonal matrices â, b, n, and n as illustrated in Table 1. This table also shows one possible ordering of the coefficients A i mn as elements of the truncated vector A. For a given configuration of balls, the truncated matrix Û has to be computed only once. If the parameters a i and b i of Robin boundary conditions are fixed, the truncated matrix W in Eq. [START_REF] Seki | Rigorous calculation of electric field effects on the free energy change of the electron transfer reaction[END_REF] or Eq. ( 64) has to be computed only once by a numerical inversion. When M is large, this is the most time-consuming operation. Once the truncated matrix W is found, the coefficients A i mn and the resulting Green function are computed rapidly. In particular, the Green function G(x, y) can be evaluated at any spatial points x and y with a low computational cost. Note also that many other diffusion characteristics such as the escape probability, the residence time, and the mean first passage time (see Sec. 4) are immediately accessible from the computed A i mn and Û.

We implemented the computation of the Green function in three-dimensional domains with disconnected spherical boundaries as a Matlab package "GreenBallsL" that can be freely downloaded at https://pmc.polytechnique.fr/pagesperso/dg/GBL/gbl.html

In this package, one needs to specify the radii, positions and surface properties (coefficients a i and b i ) of the non-overlapping balls, as well as sets of points x and y, at which the Green function should be calculated. In spite of the mathematical condition (2) needed to formally prove the convergence of the solution, the package allows one to consider touching balls as well. . . .

Table 1:

Ordering the coefficients A i mn as elements of the truncated vector A of size M = N (nmax + 1) 2 (for an exterior problem), where N is the number of balls and nmax is the maximal degree of spherical harmonics. The diagonal elements of the matrices n, â and b, involved in the boundary conditions [START_REF] Traytak | Diffusion-controlled reactions in an electric field: Effects of an external boundary and competition between sinks[END_REF], are also shown.

Monopole approximation

When the absorbing sinks are small, one can resort to the monopole approximation (MOA) which consists in truncating all expansions to the zeroth degree: n max = 0. For diffusion problems, this approximation was first proposed by Borzilov and Stepanov to study the growth of N drops immersed in an unbounded gas medium [START_REF] Borzilov | On derivation of the equation for condensation of an array of drops[END_REF] and by Deutch et al. to get approximate solutions of the trap problem in regular arrays with N ideal sinks [START_REF] Deutch | Competitive effects in diffusion-controlled reactions[END_REF]. Later on, this approximation was often employed by many authors (e.g., see [START_REF] Traytak | The diffusive interaction in diffusion-limited reactions: the steady-state case[END_REF][START_REF] Galanti | Theory of diffusion-influenced reactions in complex geometries[END_REF][START_REF] Biello | Competitive effects between stationary chemical reaction centres: A theory based on off-center monopoles[END_REF] and references therein). For the exterior problem, one only needs the elements

Û ij 0000 = R i L ij (i = j), V i 00 = R i 4πL i , (125)
while Eq. ( 40) is reduced to the set of N linear equations on A i 00 :

(a i + b i ) A i 00 R i + a i N j( =i)=1 A j 00 L ij = a i 4πL i . (126)
The monopole approximation accounts for the inter-sink distances L ij but fully ignores the angular part. The monopole approximation for the interior problem is summarized in Appendix B.

Numerical analysis and validation

In order to illustrate the efficiency of the proposed method, we consider two basic configurations of sinks.

Two concentric spheres

We start with the case of two concentric spheres for which an explicit solution in Eq. ( 90) was derived. First, we evaluate how the contribution of the n-th term in Eq. ( 90) decreases with n. This analysis assesses the accuracy of the truncated explicit solution that will serve as a reference point to check the accuracy of our numerical implementation of the GMSV. Figure 2(a) shows that the contribution of the n-th term decreases exponentially fast with n, regardless of the type of boundary condition used. In particular, the truncation size of n = 40 provides the accuracy of the order of 10 -7 which is enough for our illustrative purposes. Now, we use the analytical solution as a reference to check the accuracy of our implementation of the GMSV for the same geometric configuration. For this purpose, we compute the Dirichlet Green function G(x, y) inside the domain Ω + between two concentric spheres of radii R 1 = 1 and R 0 = 2 in two ways: analytically via Eq. ( 90) and numerically according to Eq. (66) truncated at n max . The starting point y is fixed at (0, 0, 1.5). The Green function is computed on a set of 10 000 points x uniformly distributed in the domain. The maximal error, i.e., the L ∞ -norm of the difference between analytical and numerical solutions, is then evaluated. Figure 2(b) shows the maximal error as a function of the truncation degree n max . One can see that the error decreases exponentially fast.

Co-axial configurations of spheres

Now we switch to a co-axial configuration of balls that are englobed by a larger ball. This type of configurations is particularly suitable because the axial symmetry facilitates both a visualization of the obtained results and a numerical solution by a finite element method that we use as an independent verification scheme. In fact, one can use cylindrical coordinates, (z, ρ, φ), to reduce the original three-dimensional problem to an effectively two-dimensional problem if the boundary conditions do not depend on the angular coordinate φ. For illustrative purposes, we investigate the stationary concentration of particles, n B (y), that are constantly released from a source at the outer sphere ∂Ω 0 and diffuse towards partially reactive inner sinks. According to Eq. (120), this concentration is proportional to the absorption probability p 0 (y). Setting n 0 = 1, we focus on the latter quantity. On one hand, we solve the boundary value problem (119) by a FEM implemented in Matlab PDE toolbox. As a numerical method, the FEM provides an approximate solution whose accuracy depends on the maximal mesh size h max used. To control the accuracy of the FEM, we compute p 0 (y) with two values of h max : 0.05 and 0.02. On the other hand, we calculate p 0 (y) from Eq. ( 111) by using the GMSV and computing the underlying matrices. In the following, we analyze how the accuracy of the GMSV depends on the truncation degree n max and on the reactivity of the spherical sinks.

Figure 3 shows the absorption probability p 0 (y) for the domain Ω + composed of two inner spherical sinks of radii R 1 = R 2 = 1 centered at (0, 0, ±2), englobed by the outer source of radius R 0 = 5 centered at (0, 0, 0). Setting the Dirichlet boundary condition Figure 4 shows the difference between the solutions obtained by the FEM and by the GMSV. In the top panel, the solution by the GMSV is compared to the coarser FEM solution with h max = 0.05. Increasing the truncation degree n max from 1 to 7, one progressively improves the accuracy of the GMSV solution. The maximal absolute difference (i.e., the L ∞ -norm) is reported in Table 2. One can see that this difference stops to decrease for n max ≥ 5. This reflects the fact that the difference is further determined by the limited accuracy of the FEM solution. In the bottom panel of Fig. 4, a more accurate FEM solution with the maximal mesh size 0.02 is used for comparison. In this case, the maximal absolute difference progressively decreases for all considered n max up to 7. Similar behavior is observed for mixed Dirichlet-Neumann boundary conditions set on two inner sinks (Fig. 5).

Qualitatively, the accuracy of the FEM solution with h max = 0.05 is comparable to that of the GMSV with n max = 5. However, this FEM solution involves the triangulation of the planar computational domain with 99 494 triangles and K = 50 190 vertices and thus requires to solve numerically the system of K linear algebraic equations. In turn, finding the GMSV solution relies on solving the system of 3(5+1) 2 = 108 linear algebraic equations. In addition to a 500-fold reduction in the number of equations, the GMSV provides the solution in an analytic form that can be easily manipulated. Moreover, we chose here the co-axial configuration of sinks just to facilitate the use the FEM solution in a planar computational domain. In general, a three-dimensional computational domain has to be discretized that would drastically increase the number of linear equations in the FEM (to keep the same h max and thus the same accuracy). In contrast, the computational cost of the GMSV does not depend on whether the configuration is co-axial or not. Finally, solving an exterior problem for two absorbing sinks without an outer boundary by the GMSV involves even a smaller system of linear algebraic equations, whereas the addition of an artificial outer boundary is mandatory for the FEM. We conclude that the GMSV significantly outperforms the FEM for three-dimensional domains with disconnected spherical boundaries, and is particularly valuable for exterior problems.

Conclusion

Using the classical translational addition theorems for solid harmonics, we elaborated a general semi-analytical solution for boundary value problems associated to the Laplace 2.

operator in arbitrary configurations of non-overlapping balls in three dimensions. We considered both exterior and interior problems with the most common Dirichlet, Neumann, and Robin boundary conditions. We also treated the conjugate boundary value problems with diffusive exchange between interior and exterior compartments. In all cases, the solution is based on the derived semi-analytical formula for the Green function G(x, y), in which the dependence on points x and y enters analytically through solid harmonics ψ ± mn while the associated coefficients are obtained numerically by truncating and solving the established system of linear algebraic equations. In other words, although the solution is exact, its practical implementation requires matrix truncation and inversion. The desired accuracy of the solution is achieved by varying the truncation degree. The natural choice of solid harmonics as basis functions that respect intrinsic symmetries of the domain, implies a very rapid convergence of the numerical solution, as confirmed with several examples. Even the truncation to the zeroth degree, n max = 0, known as the monopole approximation, can yield accurate results, especially when the balls are small as compared to the inter-ball distances. Moreover, the computation does not involve meshing of the domain that is often a limiting factor, especially in three dimensions. Once the coefficients in front of solid harmonics are found, one can easily and rapidly evaluate the Green function at any point of the domain. Since irregular solid harmonics decay at infinity, there is also no need for imposing an artificial outer boundary to transform an exterior problem to an interior problem that is needed for most other methods. The long range character of the fundamental solution G(x, y) implies that the impact of such an artificial boundary onto the solution can be significant even for distant boundaries. To reduce this impact in conventional methods, one would need to put the outer boundary far away from balls that would greatly increase the number of discrete elements and thus the number of linear equations. In contrast, the GMSV is even simpler for exterior domains and provides superior computational efficiency. To summarize, the major advantages of the GMSV are: semi-analytical form of the solution, mesh-free computation, very rapid convergence, and no need for imposing artificial outer boundary to treat exterior problems.

The Green function is also the key ingredient to access various characteristics of stationary diffusion among partially reactive sinks such as reaction rates, escape probability, harmonic measure, residence time and mean first passage time, for which we provided semi-analytical formulas. Although our main focus was on applications to diffusioninfluenced chemical reactions, the proposed method is also valuable in other fields in which the Laplace and Poisson equations are relevant. For instance, one can describe molecular motion in biological tissues and heat transfer in heterogeneous media, both modeled as non-overlapping balls (e.g., cells or tumors) immersed in an exterior space. The exchange between these compartments is accounted via conjugate boundary conditions. In electrostatics, the Dirichlet Green function G(x, y) can be interpreted as the electric potential created by a point charge at y in the presence of grounded balls. In fluid dynamics, one can compute the velocity potential of an incompressible flow in a pack of non-overlapping spheres which is often used as a basic model of heterogeneous porous media.

In spite of our focus on domains with disconnected spherical boundaries, the GMSV is applicable to other canonical domains and their combinations [START_REF] Traytak | Diffusion-influenced reaction rates for active "sphere-prolate spheroid" pairs and Janus dimers[END_REF]. For instance, one can consider spherical sinks englobed by a parallelepiped or by a cylindrical tube; moreover, spherical sinks can be replaced by spheroids, parallelepipeds, cylinders, or their combinations. As a consequence, such combinations of canonical domains provide a very flexible and versatile tool for modeling structured or disordered media and the related diffusionreaction processes. Future developments of the GMSV for other canonical domains is a promising perspective of the present work. Recalling Eq. ( 31), one gets a simpler form ∂g(x, y) ∂n x ∂Ωi = 1 R i n,m (2n + 1)A i mn -n V i mn ψ - mn (R i , θ i , φ i ).

(A.7)

Combining these results, we get Eq. ( 34) for the harmonic measure density.

Appendix A.5. Integrals over balls

One can compute the integral of ψ - mn (r j , θ j , φ j ) over any ball Ω I (of radius R I and centered at x I ), which is not overlapping with the ball Ω j . In fact, denoting the local spherical coordinates associated to Ω I as (r I , θ I , φ I ), one can use the I→R addition theorem (20b) for r I < L Ij to write We only consider the case when each ball Ω i can be either included into Ω I (i.e., Ω i ⊂ Ω I ), or lie outside Ω I (i.e., Ω i ∩ Ω I = ∅). In other words, we do not allow the ball Ω I to cut any ball Ω i . In this case, the integral over ΩI is simply the integral over Ω I minus the integrals over each Ω i . First, we have Ωj dx ψ - mn (r j , θ j , φ j ) = δ n0 δ m0 2πR 2 J (A.19) (although ψ - mn is singular at r j = 0, this singularity is integrable for n = 0 due to the radial weight r 2 , whereas the symmetry of the integration domain Ω j cancels the contribution from other harmonics with n > 0). Second, the integral of ψ - mn (r j , θ j , φ j ) over Ω i (with i = j) is given by Eq. (A.13). Combining all these results, we get ΩI dx ψ - mn (r j , θ j , φ j ) = 4π δ n0 δ m0

R 2 I -L 2 Ij -R 2 J 2 + U (-j,+I) mn00 L 3 Ij 3 - i R 3 i 3 U (-j,+i) mn00 , (A.20)
where the last sum is taken over the balls Ω i (except Ω j ) which are included in Ω I . This formula allows one to integrate the solution over any ball Ω I that does not cut balls Ω i . Using the addition theorem (20c), one can compute an integral over a large sphere ∂Ω I that englobes a ball Ω j . In fact, since R I > L Ij because Ω j ⊂ Ω I , one has = 1. Note that this result depends neither on the location, nor on the radius of the ball Ω j .

Appendix B. Monopole approximation for interior problems

The monopole approximation for the interior problem of finding chemical reaction rates was discussed in [START_REF] Galanti | Theory of diffusion-influenced reactions in complex geometries[END_REF][START_REF] Galanti | Reaction rate of a composite core-shell nanoreactor with multiple nanocatalysts[END_REF][START_REF] Galanti | Diffusion-influenced reactions in a spherical cavity[END_REF]. Here, we briefly present its extension for computing the Green function.

For the interior problem, one needs to modify the elements of Û and V corresponding to the outer boundary ∂Ω 0 :

Û i0 0000 = R i (i > 0), Û 0j 0000 = 1 R 0 (j > 0), V 0 00 = 1 4πR 0 . (B.1)
With this modification, the boundary conditions read

(a i + b i )A i 00 + a i R i N j( =i)=1 L -1 ij A j 00 + a i R i A 0 00 = a i R i 4πL i0 (i =
A i 00 + c i N j( =i)=1 L -1 ij A j 00 - c i C N j=1 A j 00 N k( =j)=1 c k L -1 kj = c i 4πL i0 - c i C N k=1 c k 4πL k0 - 1 4π , (B.

8) or

A i 00 1 -

c i i + c i N j( =i)=1 A j 00 L -1 ij - c i j = c i 4π 1 L i0 + 1 - 1 C N k=1 c k L k0 , (B.9)
where we denoted

Figure 1 :

 1 Figure 1: (a) Illustration of an unbounded exterior domain Ω -= R 3 \ N i=1

  (-n)n (-n+1)n ... nn 2n+1 . . .

1 a 1 a 1 a 1 3 . 1 a N a N a N 3 . 1 b N b N b N 3 .

 131313

Figure 2 :

 2 Figure 2: (a) The maximal contribution of the n-th term in Eq. (90) computed numerically for two concentric spheres of radii R 1 = 1 and R 0 = 2, with the starting point y = (0, 0, 1.5) and three combinations of boundary conditions at inner/outer spheres: Dirichlet-Dirichlet, Dirichlet-Neumann, and Neumann-Dirichlet. (b) The maximal error (L∞-norm) of the Dirichlet Green function obtained via our numerical implementation of the GMSV for the same configuration, as a function of the truncation degree nmax. The numerical solution is compared to the exact formula (90) truncated at n = 40.

Figure 3 :

 3 Figure 3: The absorption probability p 0 (y) for the domain composed of two inner spherical sinks of radii R 1 = R 2 = 1 centered at (0, 0, ±2), englobed by the outer source of radius R 0 = 5 centered at (0, 0, 0), with a 0 = 1 and b 0 = 0. On two inner sinks, we set either Dirichlet conditions (a 1 = a 2 = 1, b 1 = b 2 = 0, left), or Robin conditions (a 1 = a 2 = 1, b 1 = 0.5, b 2 = 2, middle), or Dirichlet-Neumann conditions (a 1 = b 2 = 0, a 2 = b 1 = 0, right).

  (a 0 = 1 and b 0 = 0) at the outer source, we compare several combinations of boundary conditions at the inner sinks: two fully absorbing sinks (a 1 = a 2 = 1 and b 1 = b 2 = 0), two partially reflecting sinks (a 1 = a 2 = 1, b 1 = 0.5, b 2 = 2), and one absorbing sink with one reflecting obstacle (a 1 = b 2 = 1 and a 2 = b 1 = 0). These solutions are obtained by the GMSV with the truncation degree n max = 7.

Figure 4 :

 4 Figure 4: Difference between the absorption probabilities p 0 (y) obtained by the GMSV and by the FEM of Matlab PDE toolbox, for the domain composed of two inner balls of radii R 1 = R 2 = 1 centered at (0, 0, ±2), englobed by the outer source of radius R 0 = 5 centered at (0, 0, 0). We set Dirichlet boundary conditions: a 1 = a 2 = a 0 = 1 and b 1 = b 2 = b 0 = 0. Top/bottom rows correspond to two maximal mesh sizes hmax of the FEM: 0.05 (coarser) and 0.02 (finer). Plots from left to right correspond to different truncation degrees of the GMSV: nmax = 1, 3, 5, 7. The maximal absolute errors are reported in Table2.

Figure 5 :

 5 Figure 5: Difference between the absorption probabilities p 0 (y) obtained by the GMSV and by the FEM of Matlab PDE toolbox, for the domain composed of two inner balls of radii R 1 = R 2 = 1 centered at (0, 0, ±2), englobed by the outer source of radius R 0 = 5 centered at (0, 0, 0). We set Dirichlet-Neumann boundary conditions at inner balls: a 1 = b 2 = a 0 = 1 and b 1 = a 2 = b 0 = 0. Top/bottom rows correspond to two maximal mesh sizes hmax of the FEM: 0.05 (coarser) and 0.02 (finer). Plots from left to right correspond to different truncation degrees of the GMSV: nmax = 1, 3, 5, 7. The maximal absolute errors are reported in Table2.

(- 1 ) 1 R

 11 from which Eq.[START_REF] Traytak | Methods for solution of the boundary value problems in domains with disconnected boundary[END_REF] follows. If x i = 0, then this formula is reduced toG(x, y) x , y } n max{ x , y } n+1 . (A.3)In the opposite case r i > L i , one has r > = r i and r < = L i so that1 m ψ + (-m)n (L i , Θ i , Φ i ) ψ - mn (r i , θ i , φ i ), (A.4)from which Eq. (59) follows.Appendix A.2. Derivation of the harmonic measure densityTaking the derivative of Eq. (25) with respect to r i , one finds∂G(x; y) ∂n x x∈∂Ωi =n,m nV i mn ψ - mn (1, θ i , φ i ). (A.5)Similarly, the derivative of Eq. (24) with respect to r i yields ∂g(x, y)∂n x ∂Ωi = ∂g i (r i , θ i , φ i ; y) j (r j , θ j , φ j ; y) mn (r i , θ i , φ i ) + mn (r i , θ i , φ i ) i mn ψ - mn (R i , θ i , φ i ) -n mn (R i , θ i , φ i ) mn (R i , θ i , φ i ),that can also be written as ∂g(x, y) ∂n x ∂Ωi = i n,m (2n + 1)A i mn -n ÛA i mn ψ - mn (R i , θ i , φ i ).(A.6)

mnkl 2πδ l0 δ k0 (R 2 I -L 2

 22 mn (r j , θ j , φ j ) = Ij , Θ Ij , Φ Ij ), (A.[START_REF] Miller | Symmetry and separation of variables[END_REF] where L Ij = x j -x I , (L Ij , Θ Ij , Φ Ij ) are the spherical coordinates of L Ij , and the mixed-basis elements are given by Eq. (22b). Similarly, the integral over the sphere ∂Ω I reads∂Ω I ds ψ - mn (r j , θ j , φ j ) = 4πR 2 I ψ - mn (L Ij , Θ Ij , Φ Ij ). (A.14)Now we consider a more complicated situation when Ω j ⊂ Ω I . We split the integration domain Ω I into two subsets, Ω < I and Ω > I , such thatΩ < I = {x ∈ Ω I : xx I < L Ij }, Ω > I = {x ∈ Ω I : xx I > L Ij }. (A.15)In each subset, we can use the appropriate addition theorem to compute the integral. Using Eq. (20b) for r I < L Ij and Eq. (20c) for r I > L Ij , we haveΩ < I dx ψ - mn (r j , θ j , φ j ) = mn (r j , θ j , φ j ) = kl (r I , θ I , φ I ) Ij ) = δ n0 δ m0 2π(R2Imay also need to compute the integral of ψ - mn (r j , θ j , φ j ) over Ω I without any ball Ω i : ΩI = Ω I \

  mn (r j , θ j , φ j ) = kl (r I , θ I , φ I ) = 4πR I δ n0 , (A.21) the last equality coming from the rotation symmetry of spherical harmonics Y kl and from the identity U (-I,-i) 0000

1 ,

 1

Table 2 .

 2

		h max \ n max	1	3	5	7
	DD	0.05 0.02	0.0965 0.0166 0.0099 0.0099 0.0965 0.0166 0.0028 0.0005
	DN	0.05 0.02	0.0703 0.0137 0.0137 0.0137 0.0704 0.0117 0.0019 0.0004

Table 2 :

 2 Maximal absolute errors between the absorption probabilities p 0 (y) obtained by the GMSV with the truncation degree nmax and by the FEM of Matlab PDE toolbox with the maximal mesh size hmax (see Figs.4, 5).

Acknowledgments

DG acknowledges the support under Grant No. ANR-13-JSV5-0006-01 of the French National Research Agency.

Appendix A. Technical derivations

Appendix A.1. Newton's potential

We use the Laplace expansion for the Newton's potential [START_REF] Epton | Multipole translation theory for the three-dimensional Laplace and Helmholtz equations[END_REF],

where where we used Eqs. (109, 115). According to Eq. ( 42), the derivative of A j 00 can be expressed as a linear combination of the derivatives of V k mn . We show that

from which Eq. ( 116) follows. Indeed, for j = i, the integral is

For j = i, we use the addition theorem (20b) to get

. Residence time

We use Eqs. (32, 25, 20b) to write the residence time T in a ball Ω I of radius R I centered at x I as

where L Ij = x j -x I , (L Ij , Θ Ij , Φ Ij ) are the spherical coordinates of L Ij , L I = y -x I , and V I mn is given by Eq. ( 26) which is modified for the ball Ω I .

If a 0 = 0, one can express A 0 00 from the last equation and substitute it into the former ones that yields a closed system of linear equations on A i 00 for i = 1, N :

Finally, if a 0 = 0 (i.e., the Neumann boundary condition at the outer boundary), the last relation in Eq. (B.2) is reduced to

In this case, Eqs. (B.2) can be written as

where C = c 1 + . . . + c N . Expressing A 0 00 from this relation, one gets a closed system of linear equations on A i 00 for i = 1, N :