The generalized method of separation of variables for diffusion-influenced reactions: Irreducible Cartesian tensors technique
Sergey D. Traytak


AIP/123-QED



 








sergtray@mail.ru
Semenov Federal Research Center for Chemical Physics, Russian Academy of Sciences, 4 Kosygina St., 119991 Moscow, Russian Federation








             


Motivated through various applications of the trapping diffusion-influenced reactions theory in physics, chemistry and biology, this paper deals with irreducible Cartesian tensors (ICT) technique within the scope of the generalized method of separation of variables (GMSV). Presenting a survey from the basic concepts of the theory, we spotlight the distinctive features of the above approach against known in literature similar techniques. The classical solution to the stationary diffusion equation under appropriate boundary conditions is represented as a series in the ICT. By means of proved translation addition theorem we straightforwardly reduce the general boundary value diffusion problem for N spherical sinks to the corresponding resolving infinite system of linear algebraic equations with respect to unknown tensor coefficients. These coefficients comprise explicit dependence on the arbitrary three-dimensional configurations of N sinks with different radii and surface reactivities. Specific application of the ICT technique is illustrated by some numerical calculations for the reactions, occurring in the arrays of two and three spherical sinks. For the first time we treat the statement and solution to the above reaction-diffusion problem, applying the smooth manifold concept.


The generalized method of separation of variables for diffusion-influenced reactions: Irreducible Cartesian tensors technique
    Sergey D. Traytak
    April 29, 2024
=============================================================================================================================

































§ MOTIVATION




It is presently accepted that for the most part, chemical reactions in micro-heterogeneous liquid media are contact diffusion-influenced <cit.>. The latter means that the rate of these reactions is determined significantly by the rate of encounter of reactants due to diffusion; that is, the reaction rate is not controlled by only the chemical requirement of overcoming an activation energy barrier. 

Due to their high abundance, diffusion-influenced reactions play a decisive role for a wide and diversified range of applications often occurring in physics, chemistry, biology and nanotechnology. Examples include excitation quenching of donors by acceptors, heterogeneous catalysis, crystal defect annealing, crystal growth, Ostwald ripening, evaporating or burning of liquid droplets, nutrient consumption by living cells and cell metabolism to name just a few.

To describe microphysics of diffusion-influenced reactions theoretically often so-called Smoluchowski's trapping model is used, i.e., assuming that particles diffuse in a continuum medium containing immobile uniformly distributed absorbing sinks. Smoluchowski's theory is essentially a one-sink theory, which does not account for the influence of neighboring sinks. Quite apparently the presence of other neighboring absorbing sinks decreases the local concentration of diffusing particles, and so each sink captures less particles than it would a single isolated. This, obviously, implies that many-sink effects should be incorporated into theoretical analysis. According to accepted terminology, many-sink effects are called the diffusive interaction (for an extended discussion see Sec. <ref>). Corresponding mathematical problems are well-defined, however, the theoretical description of the diffusive interaction is mostly so complicated that solving of such problems analytically at least approximately is very difficult if it is even possible. 

We shall investigate here micro-heterogeneous medium in which spherical obstacles including absorbing sinks, of different radii are distributed in a second distinct continuous phase. The dispersed obstacles assumed to be non-overlapping and impenetrable.

It has been shown explicitly that the diffusive interaction effects become the most profound in the steady state <cit.> and it is this regime that will be studied in this work. A number of different analytical approaches proposed earlier to solve the posed diffusion problems for the steady state were reviewed in Ref. <cit.>. The most complete list of the commonly used analytical and numerical methods to solve the posed problems along with some comments was given later in Refs. <cit.>.

Here, we consider a modification of the GMSV that goes back to Rayleigh's seminal paper on the conductivity of heat and electricity in a medium comprising cylindrical or spherical sink systems <cit.> (see Sec. <ref>). Although the GMSV by means of ICT seems to be fairly convenient at least for numerical calculations, however it still did not draw due proper attention to researchers. The work extends our previous studies via the GMSV <cit.> to the general case of the steady-state diffusion-influenced processes. We focus on the development a detailed theory of the diffusive interaction effects (see Sec. <ref>) by means of the ICT technique within the framework of microscopic Smoluchowski's trapping model. It turns out that  solution in terms of the ICT elucidates rather subtle mathematical facets of the multipole method and, what is more important, it significantly simplifies the calculation of the diffusion field and corresponding reaction rates. 

This research contains all the pertinent mathematical details including terminology, definitions, geometrical structure, convergence and step by step description of the GMSV algorithm with the ICT technique to solve the general diffusion boundary value problem within the scope of Smoluchowski's trapping model.

Therefore, the main purpose of this research is twofold. First, we intend to pose rigorously and in full details the boundary value problem for the most general case of boundary conditions, describing the small particles diffusion and absorption in a medium comprising both static inert obstacles and reactive sinks. Secondly, providing the overall view onto the GMSV, we solve the posed problem by means of the ICT technique, keeping in mind its practical aspects. Besides, since reviews about the problem are lacking in the existing literature, so it is our intention to provide a rather comprehensive overview of the current research activities in the field. Immediately, note here that, it seems expedient to give brief literature surveys on the subjects under consideration in the appropriate places of the text.

Thus, we hope that, first of all, this paper will repay the attention of mathematically inclined theoreticians working in a wide range of the applied scientific fields dealing with many-body effects that occur in diffusion-influenced reactions and other related topics.






§ INTRODUCTION




To provide a self-contained presentation, we start with a brief survey of the microscopic Smoluchowski trapping theory of the diffusion-influenced reactions. Without going into details, we particular emphasis on the main mathematical ideas in order to generalize easily its results to the most general case studied in this paper.

The standard notation: ℕ, ℤ, ℝ and ℝ_+ the sets of natural, integer, real and strictly positive real numbers will be used in the sequel. For the three-dimensional (3D) Euclidean space we use common ℝ^3. 
As is customary, let ∂Ω denotes the closed boundary of a 3D domain (connected open subset) Ω⊂ℝ^3 such that its closure is Ω= Ω∪∂Ω. 

We shall study real-valued functions f: Ω→ℝ_+. Let m ∈ℕ∪{ 0 } (particularly C(Ω):= C^0(Ω)) then C^m (Ω) stands for the class of m times continuously differentiable functions on  Ω and C^m(Ω) denotes a space of functions having m continuous derivatives in Ω⊂ℝ^3 and continuously extendable to the closure Ω of Ω. A function of class C^∞ (Ω) is an infinitely differentiable functions. In the following, the term smooth function  on  Ω refers to a
C^∞ (Ω)-function.

Let 𝐚 =
(a_1,a_2,a_3) ∈ℝ^3 be a fixed translation vector (bold symbols are used for vectors or tensors in ℝ^3). For any point 𝐱∈ℝ^3 define the translation operator (and its inverse one) T_𝐚 (T^-1_𝐚): ℝ^3 →ℝ^3 by the relations

    T_𝐚𝐱 = 𝐱 + 𝐚 ,    T^-1_𝐚𝐱 = 𝐱 - 𝐚 .

It is well-known that translations form a subgroup of the 3D Euclidean motion group  E(3) such that: corresponding composition operator  T_𝐜 = T_𝐚 T_𝐛 with translation vector 𝐜 = 𝐚 + 𝐛; and identity operator I:=T_0 = T^-1_𝐚 T_𝐚 <cit.>.





 §.§ A model description



To study diffusion-influenced reactions theoretically so-called trapping model is often used <cit.>. It has focused on studying diffusion of very small (point-like) reactants B in heterogeneous media, in which the disperse phase behaves as a collection of immobile obstacles immersed in an inert continuous medium commonly known as host medium. This is the case, e.g., when obstacles diffuse much more slowly than reactants (they are of much larger size). Provided arbitrary but finite number of obstacles (sinks) are regularly or randomly distributed within a given region of surrounding medium, for this system of obstacles terms "ensemble" or "array" are usually used.

Standard treatments of the theory of diffusion-influenced reactions generally assume that the host medium is unbounded (open system <cit.>), quiescent, isotropic and homogeneous, therefore mathematically it is well modelled by the 3D Euclidean space ℝ^3 (called physical space below). In turn reactants assumed to be identical non-interacting the point-like Brownian particles (from now on called B-particles or B-reactants) diffuse among 3D non-overlapping immobile spherical obstacles embedded in some domains from the physical space.

In considering microscopic kinetics of diffusion-influenced reactions, it is necessary to distinguish between obstacles. Each obstacle is characterized by its geometrical and physicochemical parameters, which drastically influence on the reaction kinetics. 
Any particle of kind A refers to the absorbing obstacle are usually called trap or sink  <cit.>.
Further, for definiteness sake, we use the term "obstacles" to refer to the fully reflecting ones only, otherwise we shall use the most common term "sinks" solely. Despite its relative simplicity, the trapping problem is fundamental for the understanding of kinetics of diffusion-influenced reactions <cit.>.

In this paper we deal with irreversible bulk contact diffusion-influenced reactions described by the simplest reaction
scheme <cit.>

    A+Bk⟶  A · B k_in⟶ A + Products  .

Here A· B is so-called encounter pair;  k and k_in stand for the reaction rate constant and  the intrinsic reaction rate constant for the reaction of the encounter pair A· B to form inert products, respectively <cit.>. Clearly, the rate k_in describes surface reactivity of sinks A <cit.>, which we assume to be spherically symmetric.

Within the scope of scheme (<ref>) the macroscopic action mass law equations for the bulk concentrations of reactants c_A(t ) and c_B (t ) are

    t_r ċ_B( t ) = - c_B ( t )  ,   ċ_A( t ) =0  ;

where t_r := (k c_A )^-1 is the characteristic time scale of the macroscopic kinetics. 


And, finally, it should be emphasized that one can assume infinite capacity approximation for sinks A provided the initial concentration of B-particles appreciably exceeds that of sinks, i.e., when condition c_A(0 ) ≪ c_B (0 ) holds <cit.>. To put it otherwise, reactions (<ref>) are the catalytic type <cit.>, when activity of the catalyst A remains unchanged. In this connection, reactions (<ref>) are called sometimes catalytic degradation reactions. 





 §.§ A microscopic Smoluchowski's theory




The absorbing rate coefficient of B-particles on a given sink surface is of primary importance for the microscopic theory of  diffusion-influenced reactions (<ref>). In his pioneering work of 1917 Smoluchowski, basing on Fick's laws of diffusion, proposed a rather simple calculation of the rate coefficient to describe absorption of the B-particles by the surface of a given fully absorbing spherical sink often referred to as a test sink <cit.>. 

Let us first treat the simplest case of the steady-state diffusion among randomly distributed static spherical sinks A of the same radii R. So, for the sink volume fraction ϕ_A in dilute enough sink ensembles the following condition holds ϕ_A:=  4 π c_A R^3/3 ≪ 1 <cit.>.




In other words diffusive interaction accepted to be weak and, therefore, it is commonly assumed that sinks mostly interact with the unbounded host medium <cit.>.

Smoluchowski's theory corresponds to the limiting case of infinite dilution (ϕ_A → 0) <cit.> and, therefore, only one-sink system should be treated. 

Let Ω_1 ⊂ℝ^3 be the test sink domain and consider its compliment
Ω_1^- := ℝ^3 \Ω_1.  It is expedient to use the standard polar spherical coordinates  { O; r, θ, ϕ} (θ, ϕ are the polar and azimuthal angles, respectively) coordinate systems attached to the origin O (0) coinciding with the test sink center. Consider  a current point P ∈Ω_1^-  with position vector 𝐫 such that  r = ‖𝐫‖ its radial distance, where ‖·‖ stands for the common Euclidean norm. Plainly, for the spherical sink domain and relevant reaction surface  we have: Ω_1:= { r < R } and  ∂Ω_1 = { r = R }. 

	For structureless point-like B particles the configuration space is a 3D domain composed of all B particle positions 𝐫.

Since diffusion of B particles occur in the exterior of the sink, domain Ω_1^- is the configuration space 

and clearly B particles local concentration is spherically-symmetric n_1 (r). Often function n_1 (r) is scaled with its bulk concentration c_B so that n_1 (r)/c_B treated <cit.>.  
However, we are going to use another function namely complementary normalized local concentration, which seems to be more appropriate to describe subsequently diffusive interaction effects

    u_1: Ω_1^- →(0, 1) ,     	u_1( r) = 1- n_1( r )/c_B .

For brevity we will still call this function local concentration, if there is no confusion to be appeared.

in three dimensions there are 13 coordinate systems in which Laplace's equation is separable

According to Smoluchowski's theory, desired concentration u_1( r) is governed by the classical diffusion equation under appropriate boundary conditions. Thus, relevant exterior boundary value diffusion problem for a partially reflecting test sink reads <cit.>

    ∇^2_r u_1 = 0     Ω_1^-  ,
    	
    . [ 4 π R^2   D∂_ru_1  - k_in (1 - u_1  )  ] | _r=R = 0  , 
    . u_1 |_r →∞→ 0  ,

where D is the translation diffusion coefficient and ∇_r^2 : = 2r^-1∂_r  + ∂_r^2 is the radial part of Laplacian. For brevity, henceforth ∂_ς stands for the partial derivative ∂/∂ς with respect to the independent variable ς. Note in passing that in physical literature condition (<ref>) is commonly referred to as the outer boundary condition <cit.>.

It has been well-known that solution to this problem is <cit.>

    u_1( r) = 1- R/rJ_0^1  .

In this formula we introdused dimensionless rate J_0^1 := k_CK k_S =κ_1(1 + κ_1) < 1,   
where κ_1 := k_in/k_S. Hereinafter J_0^1 is the Collins-Kimball rate constant k_CK( k_in) normalized by the Smoluchowski rate constant  k_S given by

    k_S = lim_k_in→∞k_CK = 4 π R D = Φ_1 ( R ) c_B  ,

where  Φ_1( R ) is the total flux of B-particles on the test sink of radius R <cit.>.

We face up to another important aspect. Lee and Karplus highlighted the fact that in terms of non-normalized concentration n_1 (r) the outer boundary condition is time dependent <cit.>, i.e.

    . n_1 |_r →∞→ c_B ( t )  .
 
To clarify the matter we note that the steady state reaction regime occurs under the following condition t_D ≪ t ≲ t_r <cit.>, where t_D = R^2/3D being the characteristic encounter time <cit.> and the frequency of reaction (<ref>) takes the form  t_r^-1= 4π DR J_1 c_A. 

Strictly speaking, one should consider the above boundary value problem (<ref>)-(<ref>) as a quasi-steady state one <cit.>. That is consistent with the fact that the characteristic time scale of the change in c_B ( t ) (see Eq. (<ref>)) is t_r ≫ t_D = t_r J_1 ϕ_A. 	
Really, one can see that "fast" (microscopic) time evolution relaxation towards a quasi-steady-state occurs for times t = O(t_D) while changing of bulk concentration c_B ( t ) takes place at "slow" (macroscopic) times  t = O(t_r). 

Concluding the section note that we can treat recalled here one-sink case as a reference problem for further generalization to the case of N sinks when diffusive interaction cannot be neglected.




§ MATHEMATICAL METHODS TO STUDY DIFFUSIVE INTERACTION



Subsequently, one-sink diffusion theory was generalized in many aspects. This section focuses on previous theoretical works developing above theory to account for the diffusive interaction between sinks. 

Smoluchowski's theory, in view of Eq. (<ref>), leads to the total flux of B's into the whole ensemble of N identical uncharged sinks as the additive sum Φ_N =  N Φ_1  ( R ) <cit.>. Nevertheless, over the years it has been realized that this theory faces challenges under attempts to understand various reaction-diffusion phenomena in systems comprising many sinks with small sink-sink separations. Indeed, each sink, due to chemical reaction (<ref>) occurring upon its surface, exerts influence on the B's concentration field around other  sinks. Simply speaking an influence arises between the sinks as a result of the fact that any sink "feels" the self-consistent diffusion field of B-particles, determined by the entire array of sinks. Obviously this fact is caused by the law of conservation of  B-particles in a host medium.


It must be recognized that, up to now, there has been no established terminology concerning above effects. The reason is that since these effects have wide applications in various fields, various technical languages, including definitions and notations were used. Really, note that a variety of terminology is used in the literature even to define such kind interaction. Besides the following terms are used: (1) competition <cit.>; (2)  diffusive (diffusional) interaction <cit.>; (3) effect of multiparticle interaction <cit.>; (4) Laplace's interaction <cit.>; (5) concentration effects <cit.>; (6) collective effects <cit.>;	(7) reactive interference <cit.>; (8) shielding effects <cit.>; (9) cooperative effects <cit.> and even (10) chemical interaction <cit.>.
That is why throughout the paper we shall clarify some of relevant terms which are used in the literature.

Regarding above-mentioned note that the use of the term "collective effects" seems an unfortunate one because of its frequent utilization in particle accelerators physics. At the same time we believe that the among pointed terms, term "Laplace interaction" (or, generally speaking,  "Laplacian transport" <cit.>) best describes the essence of the interaction under consideration for the steady-state diffusion, heat transfer and Stokes hydrodynamics <cit.>. Nonetheless, for definiteness sake we, as before, use only term diffusive interaction, taking in mind the fact, that obtained results may be applied to the another counterparts <cit.>.

As far as we know  Frisch and Collins pioneered in studying "competition among sinks for the diffusing molecules" to the diffusion-influence processes, particularly growth of aerosol particles by condensation <cit.>. Among early approaches it is also worth mentioning Ham's theory based on the Wigner-Seitz cell model <cit.>. Applying a variational procedure Reck and Prager had first established rigorous upper and lower bounds on the rate of a diffusion-controlled reactions <cit.>. Detailed exposition of this approach with applications can be found in the book Ref. <cit.>.

Later, during many years, a serious effort has been mounted toward calculation of the microscopic reaction rates, taking into account multisink effects. The monopole and dipole approximations for the steady state diffusion and reactions of B-particles in dense ensembles of fully absorbing spherical sinks has been extensively investigated by many authors <cit.>. However, it was noted long ago: "... the interactions between sinks for the competitive consumption of the solutes are not exactly accounted for in these approximations" <cit.>. Using the moment scattering, Bonnecaze and Brady derived analytically the reaction rate for cubic arrays of sinks up to the quadrupole level <cit.>. Incidentally, they have emphasized there: "To accurately compute the effective reaction rate at high volume fractions, higher order many-body multipole interactions are required."   

Most attention has been concentrated on an important particular case of two sinks. While the above illustration relied on an example with two spherical sinks that could alternatively be solved analytically by using bispherical coordinates <cit.> (see also Sec. <ref>).  

For N>2 the problem was solved by the images and reflections methods analytically <cit.>. It is noteworthy that the connection between the method of images and the method of reflections was discovered in Ref. <cit.>. An important point is that within the framework of the classical  method of images, the corresponding compensating solutions are based on point charge potentials inside sinks and, thereby, allow to obviate the need for the use of addition theorems for solid harmonics <cit.>.

Concerning theoretical methods to attack the multi-particle diffusion problem Ratke and Voorhees (Ref. <cit.>) noted: "We shall develop a solution to the diffusion problem that is consistent with the particles having a spherical morphology. This can be done in a self-consistent manner using multiple scattering theory <cit.>, irreducible Cartesian tensor <cit.> or using boundary integrals and multipole expansions <cit.>.". On this head we note that multiple scattering theory <cit.> based on the microscopic monopole approximation and  multipole expansions method was performed in spherical coordinates <cit.>. In its turn the irreducible Cartesian tensor approach <cit.> is nothing more but a scalar version of well-known induced force method (see Sec. <ref>).

Later on influence of many neighboring sinks on diffusion-controlled reactions was theoretically investigated in Ref. <cit.>. Particular focus has been given there to the fact that: "... in most analytical works, only two-sphere cases were mainly considered in the calculations of physical quantities such as rate constant. This is because of difficulty in solving differential equations with more than three spheres (or with many local boundary conditions)." On this basis authors, "to avoid the mathematical difficulty encountered in analytical approach", have been solved the corresponding boundary value problem for the 3D diffusion equation by means of the finite element method <cit.>. However, it turns out that the diffusive interaction is harder to describe by this method especially in the case of fully reflecting obstacles. 

We emphasize that interest to the study of the diffusive interaction was rekindled last decades due to problems on crowded systems <cit.> and active transport. More recently, there has been a growing interest to theoretical and experimental investigations on the diffusive interaction effects of neighboring droplets during their evaporation on a surface (substrate). Particularly in recent  experiments, there was evidence that  "...neglecting the diffusive interactions can lead to severe inaccuracies in the measurement of droplet concentration..." <cit.>. This circumstance becomes really important bearing in mind the fact that "droplet evaporation on surfaces is ubiquitous in nature and plays a key role in a wide range of industrial and scientific applications..." <cit.>. Note in passing that quasi-steady-state diffusion to an assembly of slowly growing truncated sphere on a substrate may be also treated with the help of spherical multipoles <cit.>.



At the same time we especially stress that the diffusive interaction effects  under consideration one should not be confused with so-called  diffusion interaction parameter widely studied experimentally in recent years (see Ref. <cit.> and the references therein). The latter effects are the connected with influence of concentration effects to the diffusion coefficient, which significantly is manifested in highly concentrated protein solutions.

Early works on the diffusive interaction were comprehensively reviewed in Ref. <cit.>, whereas the current status of research on this subject one can find in the recent survey <cit.>.

Thus, accurate theoretical description of the problems involving diffusive interaction in multi-particle ensembles of  spherical obstacles and sinks is a long-standing challenge due to their many-body nature. 






§ STATEMENT OF THE GENERAL PROBLEM



As we pointed out above the diffusive interaction is the most manifested in the steady state regime. Accordingly, when studying the diffusive interaction this circumstance allows us to simplify significantly the mathematical problem, ignoring the time-dependent effects. 

A rigorous formulation of problems on the microscopic theory of diffusion-influenced reactions comprises: (a) specifying the geometry of the configuration space; (b) account of the reactants properties and the host medium involved; (c) using adequate diffusion system of the continuity equations and constitutive relation; (d) prescribing appropriate boundary conditions. It is important to note from the outset that by solution to the diffusion-reaction problem, we mean solely its classical solution.

Basicly our mathematical statement of the problem is similar to that given in Ref. <cit.> (attention is drawn to the fact that sinks and diffusing particles denoted there as B and A, respectively). 

The diffusive interaction substantially depends on a given configuration of reactive and inert boundaries <cit.>, so first of all we specify the geometrical part of the problem.




 §.§ Domains definition



Consider a collection of arbitrary but finitely many N ≥ 1 microscopic spherical sinks of different radii R_i, immersed within an unbounded host fluid medium ℝ^3. Hereafter the sinks are labelled sequentially by the index i = 1, N, where symbol l, m means that all integers (including {∞}) from l to m are taken their values successively, i.e.
l, m:=l, l+1, …, m, where l (< m ≤∞ ), m ∈ℤ.

Geometrically we treat the i-th sink as a 3D body occupying a spherical domain (open ball) Ω_i ⊂ℝ^3.  Denote by Ω_i the corresponding closed bounded spherical domain, supposing that all Ω_i are non-intersecting or non-touching, i.e. Ω_i∩Ω_j= ∅  for i(≠ j)=1,N. So the total domain occupying by array of N sinks is Ω^+ := ⋃_i=1^N Ω_i.
Clearly this domain has multiply connected boundary ∂Ω^+ = ⋃_i=1^N ∂Ω_i, where corresponding boundary ∂Ω _i (i-th sink reaction surface <cit.>) is the i-th connected component of the total boundary ∂Ω^+.

On the other hand, it is apparent that the diffusion problem should be posed on the whole 3D physical space ℝ^3 excepting the total domain occupying by sinks, i.e., in sink-free region <cit.>. Hence the required domain outside all sinks is the exterior domain with respect to sinks: Ω^- := ℝ^3 \Ω^+ ⊂ℝ^3 and the whole physical space ℝ^3 can be naturally partitioned into two complementary subdomains  Ω ^+ and Ω ^-. 


Plainly, Ω^- is the set of accessible configurations for B's so, according to above definition <ref>, we term this domain as a configuration space. It is evident that diffusive interaction is highly affected by geometry of the configuration space Ω^-.

	For instance consider two sinks labeled  i and j with radii R_i and R_j, as depicted in Fig. <ref>. Here we placed a particle B is into some current point P of the configuration space  Ω^- := ℝ^3 \( Ω_i ∪Ω_j ).



Let 𝐫:= OP be a vector that represents the position of a current point P in Ω^- with respect to some fixed point O. Introduce a global Cartesian coordinate system, taking point O as the origin: { O; 𝐫} and attach the standard orthonormal Cartesian basis {𝐞_α}_α = 1^3 to this origin. The relaited components of the position vector 𝐫 we denote by r_α (α=1, 2, 3) or in ordered triple form by (r_1, r_2, r_3 ). 


Using the global coordinates, we designate by 𝐫^i_0 the position of the center of i-th sink and by 𝐋_ij := O_j O_i = 𝐫^i_0  - 𝐫^j_0 the separation vector between the two sinks centers, pointing from j-th sink toward i-th sink (see Fig. <ref>).

Clearly, any i-th sink spherical domain is:
Ω_i= {𝐫∈ℝ^3: ‖𝐫 - 𝐫^i_0 ‖ < R_i } with the boundary ∂Ω _i = {𝐫∈ℝ^3: ‖𝐫 - 𝐫^i_0 ‖ = R_i }. 



We shall consider diffusion of B's in the configuration manifold Ω ^- formed by an arbitrary but finite ensemble of N spherical sinks with fixed radii {R_i}_i=1^N centered at positions {𝐫^i_0 }_i=1^N immersed in the 3D host medium. Hence the geometry of this manifold is completely determined by N-sink configuration <cit.> or  briefly microstructure <cit.>:

    X^(N):={( 𝐫^i_0,  R_i ) }_i=1^N  .
 

In Ref. <cit.> we have shown that the 3D smooth manifold is a relevant mathematical concept to describe diffusion-influenced reactions in domains with multiply connected boundaries. 
Bearing in mind further development of the GMSV to sinks of different types of canonical domains (see below definition <ref>) we consider here a general formalism based on the manifold conception <cit.>. We shall show here that the manifold conception arises simply for the above microstructure (<ref>).
Nevertheless, in our view, the much more surprising thing is that researchers in various applied sciences still have no use in boundary value problems similar to treated here this natural and convenient language.  Especially since, paraphrasing Sawyer <cit.>, one may say, "We have all worked with manifolds but, like Molière's character who had been speaking prose all his life, we have not been aware of it." 





 §.§ Manifold structure on Ω^-




Let us introduce the manifold structure on the configuration space Ω ^-. First let us introduce open complements of the sink domains Ω_i in the physical space  ℝ^3, which play an important role in our study: Ω ^-_i=ℝ^3\Ω_i. It is clear that for all arrays (N ≥ 1) we have: Ω^- = ⋂_i=1^N Ω_i^- ⊆Ω_i^-. 

Configuration space Ω ^- is a 3D smooth manifold (strictly speaking sub-manifold in ℝ^3) <cit.> since:
(a) there is its coverage ⋃_i=1^N Ω_i^- ⊃Ω ^-; (b) Ω ^- is a metric space as a subset of (ℝ^3, ‖·‖); (c) there exists smooth homomorphisms  φ_i: Ω_i^- →ℝ^3; (d) For intersection Ω ^-_ij:= Ω ^-_i ∩Ω ^-_j ≠∅ the overlap map:  φ_ij:=φ_j ∘φ_i^-1: φ_i(Ω ^-_ij) →φ_j(Ω ^-_ij) is smooth for each pair of indices i(≠ j)=1,N.

A coordinate chart (or just a chart) on Ω ^- is a pair ( Ω ^-_i,φ _i).  A collection of all charts ( Ω ^-_i,φ _i) (i=1,N) is called a smooth atlas on Ω ^-, which defines the manifold structure on Ω^-. For short in the following  we shall simply term Ω^-_i the i-th chart. It is clear that boundary ∂Ω ^- is a 2D manifold with empty boundary.


Consider the position vector 𝐫 of point P ∈Ω^-_i, then φ_i ( 𝐫) = ( r_1^i, r_2^i, r_3^i) ∈ℝ^3. These 3 real numbers are called local coordinates of point P in the chart Ω^-_i. The corresponding local Cartesian coordinate system for any i-th sink it is convenient to designate by { O_i; r^i_1, r^i_2, r^i_3 }≡{ O_i;𝐫_i} with the same orientation as the global one and the origin O_i (𝐫^i_0) at the center of the i-th sink (see Fig. <ref>). So, for any fixed i we can choose standard local orthonormal Cartesian basis {𝐞_α^i }_α = 1^3 attached to the origin O_i (𝐫^i_0), such that r^i_α are coordinates in this basis. Introdused local coordinates are useful for computation of the diffusion field in arrays wth microstructure (<ref>). 


Let us take above position vector 𝐫 for the point P and consider it with respect to two local Cartesian coordinate systems { O_i; 𝐫_i}  and { O_j; 𝐫_j}. Corresponding point P position vectors are 

    φ_i ( 𝐫) = 𝐫_i = 𝐫 - 𝐫^i_0     φ_j ( 𝐫) = 𝐫_j=𝐫-𝐫^j_0  .

It is also evident that φ _i^-1( 𝐫_i) = 𝐫_i + 𝐫_0^i = 𝐫.
Perform the transition between these local coordinates given by the overlap map φ_ij=φ_j ∘φ_i^-1:𝐫_i  →𝐫_j with the rule (see Fig. <ref>):


    𝐫_j= φ_ij(𝐫_i) = φ_j((φ_i^-1(𝐫_i))) = 𝐫
    	_i+𝐋_ij ,     𝐋_ij=𝐫^i_0-𝐫^j_0 .


Hereafter 𝐋_ij being the vector, connecting centers of sinks i and j. In this work we also use the corresponding unit vector 𝐋_ij := 𝐋_ij/  L_ij, where L_ij is the distance between above centers L_ij= ‖𝐋_ij‖ = ‖𝐫^i_0-𝐫^j_0 ‖ = 𝐫_i-𝐫_j. 
Transformations (<ref>) and (<ref>) are clearly illustrated by the following triangle commutative diagram: 


	
		(m) [matrix of math nodes,row sep=3em,column sep=4em,minimum width=2em]
		
			{O_i;𝐫_i }    {O_j;𝐫_j }

			{O;𝐫 }    
;
		[-stealth]
		(m-1-1) edge node [left] φ_i^-1 (m-2-1)
		edge [double] node [below] φ_ij (m-1-2)
		
		
		(m-2-1) edge node [right] φ_j (m-1-2);
		
		
	

 

	Note that simple linear connection between local coordinates (<ref>) arose due to the use of Cartesian coordinates. Really, the application of polar spherical coordinate systems leads to the corresponding nonlinear dependence.  	 

In view of the foregoing, the configuration space Ω ^- can be called the configuration manifold.

An additional point to emphasize is that theoretical calculations are often used the distance between surfaces of domains Ω_i and  Ω_j, i.e. h_ij := L_ij - ( R_i + R_j) (see Fig. <ref>). Hence the condition describing configurations of N sinks, which are not intersecting or touching may be written explicitly as  

    h_ij > 0      i(≠ j)=1,N .






 §.§ Governing equation




To describe the steady state microscopic diffusion-reaction model introduce a local concentration (number density) of B-reactants, n:Ω^- →ℝ_+ depends on a given microstructure X^(N) (<ref>) and also different reactivities of sinks. To describe the latter it is convenient to introduce the set of sinks reactivities k_in^(N):={k^i_in}_i=1^N, where rate k^i_in corresponds to the i-th sink reactivity. Thus, in a fixed global Cartesian coordinate system  { O; 𝐫}, extending notations of Refs. <cit.>, one can write this function as n ( P; X^(N),   k_in^(N)). However, for short we denote it just by  n_N (P). Moreover, to investigate the effects of diffusive interaction as in the case of one sink Eq. (<ref>) we introduce  normalized concentration of B particles related to their rescaled local concentration 

    u_N: Ω^- →(0, 1) ,      u_N( P) = 1- n_N( P)/c_B .


Assuming that rescaled Fick's local diffusion flux constitutive relation holds <cit.> 

    𝐣 = D ∇ u_N( P)     Ω ^-

the rescaled local concentration u_N( P ) is governed by the B-particles
conservation law which leads to the Laplace equation

    ∇ ^2 u_N( P) = 0     Ω ^-  .

Hereinafter as usual ∇ := ∑_α=1^3𝐞_α∂_α and  ∇^2 stand for the  gradient operator and the Laplacian on ℝ^3, respectively. 

To describe diffusion-reaction phenomena, Eq. (<ref>) should be subjected to the appropriate boundary conditions on the sink surfaces ∂Ω_i and condition at infinity. 





 §.§ General boundary conditions



The present work is concerned with the whole possible range of the intrinsic rate constants 0 ≤k^i_in < + ∞ and in general case we can impose the local form of inhomogeneous Robin boundary conditions 

    . [ 4 π R_i^2 D  𝐧̂_i ·∇u_N - k^i_in (f_i( P_i)-  u_N  )]  | _∂Ω_i = g_i( P_i)  ,     i=1,N ;

where 𝐧̂_i  being the outer-pointing unit normal vector with respect to ∂Ω_i (see Fig. <ref>) and point P_i ∈Ω_i. We assume for simplicity that f_i , g_i ∈ C(∂Ω_i). 


For applications it is expedient to distinguish two particular cases of the general boundary conditions (<ref>).

(i) In the case when k_in^i → 0 conditions (<ref>) give inhomogeneous Neumann boundary conditions 

    4 π R_i^2 D. 𝐧̂_i ·∇u_N| _∂Ω_i = g_i( P_i)  .


(ii) Provided k_in^i →∞, conditions (<ref>) are reduced to inhomogeneous Dirichlet boundary conditions

    .  u_N  | _∂Ω_i=f_i( P_i).


For uniqueness of a solution to Laplace's equation (<ref>) under conditions (<ref>) one should require satisfying the regularity conditions at infinity: u_N( P) = O(‖𝐫‖^-1) and ∇u_N( P) = O(‖𝐫‖^-2) as ‖𝐫‖→∞ <cit.>. These conditions are equivalent to the following regularity condition at infinity

    .  u_N  ( P) |_‖𝐫‖→∞⇉ 0

and it uniquely determines the function u_N  ( P), which is harmonic in the infinite domain Ω ^- <cit.>. Hereafter the sign ⇉ denotes the uniform limit. 


Thus, the exterior inhomogeneous Robin boundary value problem (<ref>), (<ref>) and (<ref>) is formulated. The existence and uniqueness of the classical solution to the posed problem u_N  ( P) are a matter of common knowledge <cit.>. Furthermore, it is well-known that the smoothness class of solutions u_N  ( P) depends on a boundary operator: (a) for the Robin and Neumann problems u_N (P)∈ C^2(Ω^-) ∩ C^1(Ω^-); (b) for the Dirichlet  problem  u_N (P)∈ C^2(Ω^-) ∩ C(Ω^-) <cit.>.





 §.§ Boundary conditions for diffusion-influenced reactions



It must be emphasized that the general boundary value problem (<ref>), (<ref>) and (<ref>) describes a vast amount of diffusion-reaction phenomena, occurring in nature and industry (see Sec. <ref>).
Although the approach under consideration may be applied to that problem, for the present paper we confine ourselves to study an important special case: f_i( P_i) ≡ 1 and g_i( P_i) ≡ 0, describing the microscopic Collins-Kimball diffusion-influenced reactions <cit.>.

Specific values of k^i_in should be known experimentally or by means of more detailed kinetic theory, e.g., based on the Fokker-Planck-Klein-Kramers equation <cit.>. Provided magnitudes of k^i_in are positive and finite, we deal with diffusion-influenced reactions on so-called partially reflecting (absorbing) sinks <cit.>.  In the case of very large values of k^i_in (i.e., as k^i_in→ + ∞), the sinks are called fully (perfectly) absorbing, and the corresponding reactions, occurring on them, are known as fully diffusion-controlled <cit.>. Physicall that means that the activation barrier appeared to be very small <cit.>.  Particular attention should be given to the case of fully reflecting obstacles (when k^i_in=0), which becomes also important for dense arrays of sinks and obstacles. Formally speaking, we can treat a fully reflecting obstacle as a sink with the zero intrinsic rate. 

Of special interest in many applications is the case when sinks possess different surface reactivity for a  microstructure X^(N). To put it another way when intrinsic rate constants k_in^i are different for at least some numbers from the set { 1 < i ≤ N  } <cit.>. At first glance it would seem that this case is described by boundary conditions (<ref>) automatically. However, solution of the boundary value problem (<ref>), (<ref>) and (<ref>) may be faced subtle mathematical difficulties <cit.>.

In this connection we believe that the mixed boundary value problems should be distinguished by two types of boundary conditions. We adopt the following definitions.

	If the boundary conditions differ on different parts of a given connected
	component of a boundary ∂Ω_i these conditions are referred to as proper mixed boundary conditions.

Physically the proper mixed boundary condition means that a given connected component ∂Ω _i is characterized by a prescribed heterogeneous reactivity. So from analytical viewpoint one should define a piecewise constant function k^i_in: ∂Ω_i →  [0  + ∞). 

	Provided different boundary conditions are posed on different connected components they are termed improper mixed boundary conditions.

Therefore, for the improper mixed boundary conditions we have: .  k^i_in|_∂Ω _i = const.

Figure <ref> shows a geometric sketch of a simple array consisting of three chemically distinct reaction surfaces described by the improper mixed boundary conditions.

It is significant that, generally speaking, proper mixed boundary value problems are seldom amenable to analytical treatment because of their complexity, therefore, henceforth we shall treat here only improper mixed problems, omitting word "improper". For more details regarding derivations of improper mixed boundary conditions, the reader is referred to Refs. <cit.>.


In concluding this section, we note, that crowding and geometrical constraints effects on diffusion-influenced reactions <cit.> may be treated if one considers microstructures X^(N) containing 1 ≤ N_0 ≤ N inert obstacles.

The presence of fully reflecting obstacles might strongly influence the reaction rate when the number of these obstacles is large enough <cit.> (see also numerical calculations of the reaction rates for the mixed Dirichlet-Robin problem at N=3 in Sec. <ref>).





 §.§ The dimensionless formulation of the problem



Before subsequent mathematical study it is expedient to reduce the original boundary value problem to a non-dimensional form. 
In its turn, for spatial dimensionless independent variables we shall utilize corresponding local Cartesian coordinates ( O_i; 𝐫_i ) at the i-th sink of radius R_i as a characteristic unit:

    ξ_i=𝐫_i/R_i ,   ξ^i_0=𝐫^i_0/R_i   ξ^i_ν=r^i_ν/R_i ,   ξ_i=‖ξ_i ‖ .
 
Besides these dimensionless local coordinates for any pair of sinks i and j (≠ i) the following important parameters <cit.> 

    ε_ji=R_i/L_ij < 1 ,   ε_ij=R_j/L_ij < 1

are naturally arisen. Clearly, parameters (<ref>) are totally determined by a given microstructure X^(N) (<ref>) and in general case they obey the condition ε_ij≠ε_ji. Moreover, the unit sphere centered at a point ξ^i_0 we henceforth designate by

    ∂Ω_i^1 ≡∂Ω_i (ξ^i_0; 1 ) := {ξ_i ∈Ω_i^-: ξ_i=1 } .
 



Here it should be especially noted that performed normalization (<ref>) holds only locally in charts Ω_i^- and does not hold in the whole manifold Ω^-.  Note also that in the interest of readability for all charts Ω_i^-, transformed according to Eq. (<ref>), we retain the same designations. 

It is expedient now to write the diffusion problem immediately in the local coordinates { O_i;ξ_i}. However, sometimes to emphasize the fact that we deal with the field on a manifold coordinateless form u_N (P ) will be used. Clearly, the basic equation governing exterior steady-state Eq. (<ref>) under partially reflecting  boundary conditions (<ref>) and regularity condition at infinity (<ref>)  with respect to the local concentration u_N (ξ_i ) (i= 1,N) and dimensionless variables (<ref>) reads

    ∇_ξ_i^2  u_N (ξ_i ) =0     Ω_i^- ,
    	
    
    	-. ∂ _ξ _i u_N (ξ_i ) | _∂Ω_i^1=κ_i[1 -.
    	u_N (ξ_i ) | _∂Ω_i^1] > 0   ,  
    . u_N (ξ_i ) | _ξ_i →∞⇉ 0  .
	
Hereafter  ∇_ξ_i^2 stands for the corresponding dimensionless Laplacian written in local coordinates { O_i;ξ_i}. 


Besides we introduce the corresponding dimensionless set of sinks reactivities:  κ^(N):={κ_i }_i=1^N with	κ_i=k_in^i/k_S^i ∈[0, ∞),
where k_S^i=4π R_iD is the Smoluchowski rate constant (<ref>) for the fully absorbing i-th sink.

Thus, from a mathematical viewpoint, we deal with the exterior Robin boundary value problem (<ref>)-(<ref>) in an unbounded 3D manifold Ω ^- with N-connected boundary ∂Ω ^-.  An additional point to emphasize is that the boundary conditions (<ref>), being posed on all connected components ∂Ω _i of the boundary ∂Ω^-, are integral, which reflects the important fact that diffusive interaction between sinks are not pairwise additive. 





 §.§ The microscopic trapping rate



For the theory of diffusion-influenced reactions the most important value is the microscopic trapping rate for the i-th sink defined as k_i := Φ_i/c_B, where Φ_i stands for  the total flux of B's on the i-th sink surface. 
In turn with the help of known solution n_N ( P ), the total microscopic trapping rate on the reaction surface ∂Ω_ i can be calculated straightforwardly

    k_i =∮_∂Ω_i. 𝐧̂_i ·𝐣| _∂Ω_idS_i   ,

where dS_i is the i-th surface element. 
Then, recasting (<ref>) with respect to u_N (ξ_i ) in an exterior neighborhood of the i-th sink surface, one can easily obtain the desired microscopic trapping rate (<ref>) by the surface integral over the unit sphere 

    k_i (X^(N);κ^(N)) = - ∮_∂Ω_i^1 . ∂ _ξ _i u_N | _∂Ω_i^1dξ̂_i .

Hence, for dense enough sink arrays when diffusive interaction effects become important it is necessary to determine the appropriate correction factors to the absorbsion rates on all N sinks.
Therefore the above rate (<ref>) is sought in the form

    k_i (X^(N);κ^(N)) = k_S^i  J_0^i (κ_i) J_i(X^(N);κ^(N)) ,

where the i-th dimensionless Collins-Kimball rate given by 

    J_0^i (κ_i)= κ_i(1 + κ_i)      i= 1,N .
 
generalizes one-sink expression. Below, to simplify notations, we omitted corresponding parameters.
Evidently,  J_0^i is the unperturbed by the diffusive interaction  Collins-Kimball trapping rate for the i-th sink normalized by k_S^i. 

In representation (<ref>) the magnitudes J_i are commonly called the screening coefficients <cit.> (rate correction factors <cit.>). It is absolutely clear that due to dependence on configuration of sinks X^(N) and their reactivities κ^(N) screening coefficient J_i (usually, for simplicity, we shall omit arguments in the function notation) expresses the effect of the diffusive interaction between a given i-th sink and other N-1 sinks. It emerges from the physical standpoint that the diffusive interaction reduces the reaction rate of either sink compared with that of a single one. Hence, the stronger the diffusive interaction the smaller J_i is and the following relations hold true 

    0 < J_i  <1 ,   lim_ε_ij→ 0 J_i  = 1      i(≠ j)=1,N .

Limit (<ref>) simply means that the diffusive interaction disappears at large separations between sinks.

It is significant that, generally speaking, the microscopic trapping rate (<ref>) for N ≥ 2 does not coincide with the reaction rate constant k defined by Eq. (<ref>), which is a macroscopic value <cit.>. 

So, investigating diffusive interaction effects in systems with N ≥ 2 sinks we cannot use  the term "rate constant" for k_i given by Eq. (<ref>) any more <cit.>.





§ METHOD OF SOLUTION




The present section is devoted to an overall description of the generalized method of separation of variables, which is a powerful tool for solution of the Robin boundary value problem (<ref>)-(<ref>). 

However, before starting to discuss this topic, we briefly highlight the application of the boundary integral equations method. <cit.> One can reduce the above diffusion problem (<ref>)-(<ref>) to solution of a system of Fredholm boundary integral equations of the II kind. It turned out that this system of integral equations may be reduced to the corresponding infinite system of linear algebraic equations. <cit.> Despite the fact that this method is also applicable to the mixed boundary conditions <cit.> its implementation in the diffusion problem at issue  involves a number of disadvantages compared to the approach that would be considered here. 






 §.§ On the method of separation of variables



First recall some well-known mathematical facts on the standard method of separation of variables (MSV) <cit.>. An additional point to emphasize is that present study deals with multiplicatively separable solutions of the diffusion equation only.  

Let us fix a global Cartesian coordinate system  { O; 𝐫} in a bounded or unbounded domain Ω⊂ℝ^3.

	A domain Ω is said to be canonical if any solution u (𝐫) to a boundary value problem in Ω may be derived by the standard MSV.

Thus, in other words, there exists an orthogonal functional basis of solutions {φ_n (𝐫) }_n=0^∞ (called basis solutions) to a boundary value problem in a canonical domain  Ω⊂ℝ^3 and a set of real coefficients  { C_n }_n=0^∞ such that 

    u (𝐫) = ∑_n=0^∞ C_n φ_n (𝐫)     Ω ,

where series is absolutely and uniformly convergent. For 3D Laplace' equation, e.g., canonical domains are those bounded by surfaces: spheroids, cylinders, cones, paraboloids, and their combinations for interior and corresponding complements in ℝ^3 for exterior boundary value problems, respectively.

The specific values of the coefficients C_n to be determined from the appropriate boundary conditions.
Plainly, both introduced above domains Ω_i and their complements Ω_i^- (i = 1, N) are internal and exterior canonical ones for the steady-state diffusion problems, respectively. The relevant functional basises consist of regular and irregular solid harmonics.

The application of the MSV to the well-posed diffusion problems involves:

	
  * application of the local linear superposition principle;
	
  * determination of the basis solutions to the equation in a given canonical 
	domain;	
	
  * determination of unknown coefficients from the posed boundary conditions.






 §.§ Generalized method of separation of variables




The generalized method of separation of variables (GMSV) <cit.> naturally stems from the standard method of separation of variables and its incipience go back to the classic works by Maxwell <cit.> and  Basset. <cit.> However, the most clearly this idea was expressed in Rayleigh's seminal paper on the conductivity of heat and electricity in a medium with cylindrical or spherical sinks arranged in a rectangular array <cit.>. Nowadays, this approach is referred to as the Rayleigh multipole method.

Since then, the GMSV had been intensively developed and found numerous applications in electrostatics, hydrodynamics, mechanics, heat transfer, diffraction theory and many other fields. Surprisingly, these advances typically remained "hidden" within each discipline and unknown to the researchers, worked in other fields. Such a parallel development led to multiple "rediscoveries" of the same results in different fields. Wherein, various names were given to the GMSV by different authors.  It is also widely known as the "generalized Fourier method" <cit.>. Guz' and Golovchan treated the method as essentially a particular case of the "method of series" <cit.>. Finally, in micromechanics of heterogeneous materials mainly term "multipole expansion approach" is used <cit.>. 

Note that all above names capture only some features of the method under consideration. As this method relies on the separation of variables in local curvilinear coordinates, following Ivanov, it can be call "generalized method of separation of variables" <cit.>. We have to emphasize that here a particular case of the GMSV for the sinks of the same form (spherical) will be utilized. So we shall not apply the
re-expansion theorems connecting different basis functions corresponding to the different shapes of sinks.

Among the numerous works devoted to the development of GMSV, the Ivanov's book holds a special place <cit.>.
For the first time Ivanov presented the scheme of the GMSV in full details and, moreover, he derived a number of new and reviewed known then addition theorems for the sets of basis solutions to the Helmholtz equation, written in different curvilinear coordinate systems <cit.>.  Although Ivanov used the method to solve various problems concerning diffraction theory for two particles, but it is clear that the same approach is also valid for finite number of particles. 



The basic idea of the GMSV consists in reducing the original boundary value problem posed in the configuration manifold Ω^- to the corresponding N boundary value problems in charts Ω^-_i. These N problems are coupled through the boundary conditions, prescribed on the whole disconnected boundary ∂Ω^-. Thus, the GMSV is actually the method of separation of variables on the 3D manifold <cit.>.





 §.§ Outline of approach



The application of the GMSV to the boundary value problem (<ref>)-(<ref>) may be formulated as the GMSV algorithm which in turn can be divided up into following main steps.

	
  * Proof of the well-posedness of the solution to the posed boundary value problem;
	
  * Deconposition of the desired solution in the form of the sum of partial solutions by means of the general linear superposition principle;
	
  * Determination of the appropriate basis solutions to the diffusion boundary value problem posed in the given canonical domains with disconnected boundary;	
	
  * Determination of relevant addition theorems to obtain local regular expansions of partial solutions in a neighborhood of corresponding connected components of bounday;	
	
  * Application of the addition theorems in order to satisfy the boundary conditions; 
	
  * Reduction of the problem to a resolving self-consistent infinite system of linear algebraic equations;
	
  * Solution of the resolving infinite system of linear algebraic equations;
	
  * Calculation of the microscopic reaction rate.

In Fig. <ref> we provide a full scheme of the GMSV application to a typical  diffusion boundary value problem. One can see that the linearity of the boundary value problem is one of the primary requirements for the GMSV algorithm as well.

Simple inspection reveals that the fundamental differences between standard MSV and GMSV (see Fig. <ref>) are in the blocks depicted by long-dashed lines. Hence, contrary to the MSV the GMSV algorithm requires two supplementary rather non-trivial mathematical tools: the appropriate addition theorems for the basis solutions along with methods for of the resolving infinite system of linear algebraic equations.





§ THE GMSV WITH IRREDUCIBLE CARTESIAN TENSORS



Previously, by means of GMSV for 3D domains with N disconnected spherical boundaries, we presented the detailed derivation of the Laplacian Green functions for both the exterior and interior Dirichlet and Robin problems along with the Green function for the relevant conjugate (transmission) boundary value problem <cit.>. Some facets of the GMSV applications to the theory of stationary bulk diffusion-influenced reactions on spherical sinks were also studied and discussed in Refs. <cit.>. However, it is important that all analytical and numerical calculations were performed there using the local spherical coordinate systems.

Plainly the general solution algorithm is feasible in two ways: by means of either solid spherical harmonics with respect to a polar spherical coordinate system or, alternatively, with the help of the ICT (for details, see below Appendix <ref>). The wide use of solid harmonics in the theory of diffusion-influenced reactions is a matter of common knowledge (see, e.g., Refs. <cit.> and references cited there). In the meantime, studies devoted to similar applications of ICT are only very few <cit.>. This seems surprisingly, taking into account the facts that ICT formalism has a number important advantages against standard solid harmonics approach (see discussion below in Sec. <ref>). 
In this connection, it is particularly remarkable that recent decades investigations have inspired a renewed interest for use of the ICT technique to study a number of problems in systems with N-spheres microstructure (<ref>) <cit.>. 

Against these recent trends, the lack of new works concerning the applications of the ICT in diffusion-influenced reactions theory is especially noticeable. 

To our knowledge Lebenhaft and Kapral were the first who successfully applied the ICT technique to the diffusion problem (<ref>)-(<ref>) <cit.>. Subsequently, an explicit form of the GMSV with the aid of ICT was carried out in our paper Ref. <cit.>, where we utilized this form of the GMSV to solve problems on Ostwald ripening, taking into account multipole corrections. It should be stressed that our approach stems from the method of Cartesian ansatz suggested in 1978 by Schmitz and Felderhof <cit.>. Although Ref. <cit.> studied a hydrodynamic problem, authors also treated the Robin boundary value problem for Laplace's equation in the unbounded domain exterior to a sphere. In that regard it is pertinent to cite here above paper: "We solve this problem by making an ansatz in cartesian coordinates, rather than following the usual method of introducing spherical coordinates." 
The Schmitz-Felderhof ansatz means that one is seeking a solution in the form of mixture dependence of spherical radial ξ_i and corresponding Cartesian local coordinates  ξ_γ _m^i in every i-th chart (i = 1, N) <cit.>

    u_N (ξ_i ) = ∑_n=0^∞ w_n (ξ_i  ) C_γ _1 …γ _n^i ⊙^n ξ_γ _1^i…ξ_γ _n^i^[-.04in].01in.05in
    		.47in.01in[-0.04in].01in .05in    Ω_i^-  .

Hereinafter ξ_γ _1^i…ξ_γ _n^i^[-.04in].01in.05in
	.47in.01in[-0.04in].01in .05in is the ICT defined by Eq. (<ref>). The explicit form of functions w_n (ξ_i  ) may be found by substituting in Laplace's equation (<ref>) while the unknown tensor coefficients C_γ _1 …γ _n^i to be determined from the boundary conditions (<ref>).

Note in passing that so-called self-consistent field method that is another version of the above method of Cartesian ansatz <cit.>. A method of reflections with the help of ICT used in Ref. <cit.> to solve the problem of the spherical particles motion due to a temperature gradient is also worth mentioning here. It is important to mention that some approaches vide used just different modifications of the GMSV in terms of the  above treated version with the help of the ICT <cit.>.

In Ref. <cit.> a scalar version of well-known induced force method has been applied to the theory of diffusion-controlled reactions for the first time. It turns out, however, that both induced forces method and the GMSV by means of the ICT lead to the same second kind ISLAE with respect to unknown tensor constants involved <cit.>. 





§ THE ADDITION THEOREM IN TERMS OF THE IRREDUCIBLE TENSORS





Due to the fact that the Laplace equation (<ref>) remains invariant under the action of the group of Euclidean motions <cit.> irregular solid harmonics (see definitions in Appendix <ref>) written in local coordinates of any chart Ω_i^- (Ω_j^-) can be recast in terms of regular solid harmonics written in local coordinates of another chart Ω_j^- (Ω_i^-). 

One of the key point of the GMSV algorithm is the transformation of the solid harmonics under action of subgroup of translations (<ref>) <cit.>.
Provided charts Ω_i^- and  Ω_j^- possess the same form (complements of spheres in the present paper), the appropriate transformation formulas of that kind we shall call irregular to regular translation addition theorems (I→R TAT). 

	On the contrary, in the case when forms of domains Ω_i^- and  Ω_j^- are different (like the exteriors of a sphere Ω_1^- and a prolate spheroid Ω_2^- see Ref. <cit.>) corresponding relations it is more expedient to call translation re-expansion theorems.


The goal of this section is to present a proof of the I→R TAT within the scope of the GMSV in terms of the ICT.
We note in passing that, strictly speaking, we deal here with scalar translation addition theorems only.





 §.§ The degenerate translation addition theorem



First let us focus our treatment on a particular case of the I→R TAT , which plays a key role in the subsequent proof of the general TAT.

Assume that r > ‖𝐫^i_0 ‖, where r = ‖𝐫‖ is the distance between the global origin O and a point P (see Fig. <ref>).

Thereby, we can write the known Taylor's expansion

    1/‖𝐫_i‖ =	1/‖𝐫 - 𝐫^i_0 ‖ = ∑_k=0^∞(-1)^k/k!( 𝐫^i_0 ·∇)^k ( 1/r) .

Clearly, series (<ref>) converges absolutely and uniformly for all r > 0 <cit.>. 

Function 1/r is often termed the generating function.
Traditionally, Eq. (<ref>) is referred to as the multipoles expansion of the (global) fundamental solution of the Laplacian at the pole point 𝐫^i_0 in powers of r^-1 (see, e.g., Ref. <cit.>). 

For the further consideration, it is convenient to recast expansion (<ref>) in terms of the ICT defined by Eq. (<ref>) <cit.>. This fact allows us to formulate so-called the degenerate translation addition theorem for solid harmonics 

	For all r > 0 the fundamental solution of the Laplacian by means of the ICT may be expanded in an absolutely and uniformly convergent series
	
    1/‖𝐫-𝐫^i_0‖	=∑_k=0^∞ω_k r^-(2k+1)r_γ _1… r_γ _k^
    			[-.04in].01in.05in.45in.01in[-0.04in].01in
    			.05in⊙^k r_γ _1^0i… r_γ _k^0i^
    			[-.04in].01in.05in.45in.01in[-0.04in].01in
    			.05in ,

	where ω_k:=  ( 2k-1)!!/k!,  r_γ _ν and r_γ _m^0i are Cartesian coordinates of vectors 𝐫 and 𝐫^i_0, respectively.


This form of multipole expansion is very broadly utilized across various scalar fields of the physical and chemical sciences. Moreover, (<ref>) plays a key role in derivation of the translation addition theorem and that is why it is also called the degenerate addition theorem.






 §.§ The translation addition theorem



It has been known that I→R TAT holds true for solid harmonics ψ^±_nm( 𝐫) defined in polar spherical coordinate systems (see, e.g., Ref. <cit.> and references therein). Here we shall present a derivation of the similar I→R TAT in terms of ICT X^±_n (see Appendix <ref>), which was given in Refs. <cit.> without a proof.

To calculate the local concentration u_N( P)  at
any given point P ∈Ω ^-⊂Ω ^-_i ∩Ω ^-_j with position vector 𝐫 we consider functions connected with every coordinate chart ( Ω ^-_i,φ _i): u^'( 𝐫_i) =u_N( P) for 𝐫_i=φ _i( 𝐫) (i=1, N) (see <ref>). 
In any other chart ( Ω_j ^-,φ _j) at the same
point P ∈Ω ^- we have alternatively

    u^''( 𝐫_j) =u^'( 𝐫
    	_i) =u_N( 𝐫) , 
    	
    	u^''( φ_ij(𝐫_i)) = u^'( φ_ji(𝐫_j)  )=u_N( 𝐫)  .


For our problem the overlap map φ_ij is a simple translation (<ref>) so let us start with an important general

	Any two given local curvilinear coordinate systems { O_i;ζ_i} and { O_j;ζ_j} are said to be consistent if there exists a translation (<ref>) with vector 𝐚 such that ζ_j =T_𝐚ζ_i = ζ_i + 𝐚. 

This definition immediately brings to

	Any two local Cartesian coordinate systems { O_i;𝐫_i} and { O_j;𝐫_j} are consistent if and only if: (a) they are positively oriented; and (b) all their axes are pairwise parallel with respect to the global Cartesian coordinates.

One can see that local Cartesian coordinates { O_i; 𝐫_i}  and { O_j; 𝐫_j}  depicted in Fig. <ref> are consistent.  
Note in passing that some authors prefer to use term "parallel coordinate systems" rather than "consistent coordinate systems" <cit.>.




	A translation addition theorem for a given smooth function u: Ω ^-→ℝ_+ is a formula expressing its value u(𝐫 + 𝐚) in terms of values u(𝐫) and u(𝐚),  where 𝐫, 𝐫 + 𝐚∈Ω ^-, and also their derivatives.



	The translation addition theorem I→R yields an explicit formula for converting an irregular ICT X_n^- given in one origin O_j into a local expansion with respect to the corresponding regular ICT X_n^+ about a shifted origin O_i.

The meaning of the addition theorems is to derive corresponding addition formulas. 

Now, we can state and prove our main result in this section.

	The irregular to regular translation addition theorem for the irreducible Cartesian tensors  holds true 
	
    r_j^-( 2n+1) r _γ _1^j … r _γ
    				_n^j^[-.04in].01in.05in.5in.01in
    			[-0.04in].01in .05in  
    		=∑_k=0^∞ U^ij_γ _1 …γ _nμ
    			_1  …μ _k( 𝐋_ij) ⊙^k r_μ _1^i … r _μ _k^i^[-.04in].01in.05in.5in.01in
    			[-0.04in].01in.05in ,
	
	where the transformation matrix elements read
	
    U^ij_γ _1 …γ _nμ
    			_1  …μ _k( 𝐋_ij) =σ
    		_kn  L_ij^-(n+k+1 )Λ _γ _1 …γ _nμ
    			_1  …μ _k( 𝐋_ij)  , 
    σ _kn=( -1) ^n[ 2( k+n) -1]
    			!!/k!( 2n-1) !! , 
    Λ _γ _1 …γ _nμ _1 …μ _k( 𝐋_ij) =L̂_γ _1^ij…
    			L̂_γ _n^ijL̂_μ _1^ij…L̂_μ _k^ij^
    			[-.04in].01in.05in1.04in.01in[-0.04in].01in
    			.05in ,   L̂_μ _m^ij= L_μ _m^ij / L_ij  ,    m ∈ℕ

	wherein series (<ref>) converges uniformly and absolutely if inequalities (<ref>) hold true.

 
	Rewriting Eq. (<ref>) in the relevant local Cartesian coordinates of j-th and i-th spheres we have 
	
    r_γ _ν^j=r_γ _ν^i+L_γ _ν^ij ,

	where L_γ _ν^ij are coordinates of vector 𝐋_ij (γ _ν =1,3). Linear dependence (<ref>) leads to the evident relations ∂
	_r_γ _ν^j=∂ _r_γ _ν^i. 	
	Using them in Eq. (<ref>), one readily obtains
	
    r_j^-( 2n+1) r_γ _1^j … r_γ _n^j
    		^[-.04in].01in.05in.45in.01in
    			[-0.04in].01in
    			.05in=( -1) ^n/( 2n-1) !!∂
    		_r^i_γ _1…∂ _r^i_γ _n( 1/𝐫
    			_i-𝐋_ij)
	
	or, taking into account that ∂ _r_γ _ν^i=-∂_L_γ _ν^ij,
	one has the relation
	
    r_j^-( 2n+1) r_γ _1^j… r_γ _n^j
    		^[-.04in].01in.05in.45in.01in
    			[-0.04in].01in.05in  
    		=1/( 2n-1) !!∂ _L_γ
    				_1^ij…∂ _L_γ _n^ij( 1/𝐫_i-
    			𝐋_ij)  .
	
	Theorem <ref> yields 
	
    1/‖𝐫_i-𝐋_ij‖	=∑_k=0^∞( -1) ^k/k!∂ _L_γ _1^ij…∂ _L_γ
    				_k^ij(1/L_ij) ⊙^k r_γ _1^i… r_γ _k^i^
    			[-.04in].01in.05in.5in.01in[-0.04in].01in
    			.05in .
	
	Substituting this expansion into Eq. (<ref>) and taking into consideration the uniform convergence of the series (<ref>) we arrive at the required translation addition theorem (<ref>). 

Note that, following Miller, the transformation matrix U^ij_γ _1 …γ _nμ
	_1  …μ _k is commonly termed the mixed-basis matrix <cit.>.

For further investigation of the problem (<ref>)-(<ref>) it is expedient to recast expression (<ref>) in a dimensionless form

    ξ _j^-( 2n+1) ξ_γ _1^j …ξ_γ
    			_n^j^[-.04in].01in.05in.5in.01in
    		[-0.04in].01in .05in  
    	=∑_k=0^∞U^ij_γ _1 …γ _nμ
    		_1  …μ _k⊙^k ξ_μ _1^i …ξ_μ _k^i^[-.04in].01in.05in.5in.01in
    		[-0.04in].01in.05in ,
	
where dimensionless mixed-basis matrix elements are

    U^ij_γ _1 …γ _nμ
    	_1  …μ _k:=σ
    _knε _ij^n+1ε _ji^kΛ _γ _1 …γ _nμ
    	_1  …μ _k(𝐋_ij)  .

Particularly dimensionless form of the degenerate translation addition theorem (<ref>) reads

    ξ _j^-1 =∑_k=0^∞U^ij_μ
    		_1  …μ _k⊙^k ξ_μ _1^i …ξ_μ _k^i^[-.04in].01in.05in.5in.01in
    		[-0.04in].01in
    		.05in .





The above presented proof of the I→R TAT in terms of ICT seems to be the simplest among previously known theorems for solid harmonics in terms of polar spherical coordinates <cit.>. It is important to stress that using connection between spherical solid harmonics and ICT (see Appendix <ref>) the proved addition theorem (<ref>) can be reduced  to that for the spherical solid harmonics written in a spherical coordinate system. 

Finally, note that to solve diffusion problems under proper mixed boundary conditions <cit.>, when intrinsic rates are functions on the angular local coordinates κ_i (θ_i, ϕ_i ), we need the rotational addition theorem for the ICT. The latter theorem may be proved in entirely similar way.





§ SOLUTION TO THE PROBLEM



In this section, we implement the GMSV algorithm described in Sec. <ref> using ICT technique. Moreover we shall develop here the version of the GMSV elaborated previously in Refs. <cit.>.




 §.§ The GMSV algorithm implementation



For clarity sake we organize our treatment in according to the method outline <ref> distinguishing   8 steps inherent the GMSV.


Step 1. The mixed boundary value problem (<ref>)-(<ref>) is well posed and has a unique classical solution <cit.>.


Step 2. First, we can apply superposition principle. Decompose the general solution
u_N (P ) of above diffusion problem (<ref>)-(<ref>) with respect to the i-th sink, considering the set of partial solutions to Eq. (<ref>), which are defined on the charts: { u_i (ξ_i ) }_i=1^N (u_i: Ω_i^- →(0,1]). Thus, in the local Cartesian coordinates { O_i; ξ_i} we can write

    u_N (ξ_i ) = u_i (ξ_i  ) +  ∑_j( ≠i) =1^N u_j (ξ_j  )  in  Ω^-_i  .

Here, by virtue of Eqs. (<ref>) and (<ref>), one has evident relations:  ε_ijξ_j  = ε_jiξ_i  + 𝐋_ij (see Fig. <ref>). From the physical standpoint representation (<ref>) means that to find the desired reaction rate on the i-th sink we need to know only local behavior of the field u_N (P ) in a vicinity of this sink.

Clearly, by definition functions u_i (ξ_i ) are harmonic in charts Ω ^-_i, i.e.

    ∇_ξ_i^2 u_i(ξ_i ) = 0  , 
    . u_i(ξ_i ) |_ξ _i →∞⇉ 0  .

Therefore, the original diffusion boundary value problem (<ref>)-(<ref>) in the manifold Ω^- (intersection of all charts Ω ^-_i) is reduced to N coupled problems for u_i(ξ_i ) in simpler charts Ω^-_i, with the advantages of the local Cartesian coordinates { O_i; ξ_i}.  

Alternatively, for inhomogeneous partially reflecting boundary conditions (<ref>) in an exterior ϵ-neighborhood of the boundary ∂Ω_i: Ω ^-_i(ϵ):= {ξ_i: 1 < ξ_i  < 1 +ϵ}⊂Ω ^-_i  it is convenent to recast solution u_N (ξ_i ) as

    u_N (ξ_i ) = u_( 0) ^i ( ξ_i )  + δ u_i (ξ_i  )   in  Ω ^-_i(ϵ)  ,

where u_( 0) ^i ( ξ_i ) is the unperturbed one-sink solution and δ u_i (ξ_i  ) is a perturbation due to influence of each from other N-1 sinks for j( ≠i) =1, N. Owing to Eq. (<ref>) we assume that these solutions obey inhomogeneous and homogeneous partially reflecting boundary conditions, respectively

    - . ( ∂ _ξ _i  u_( 0) ^i -κ_i u_( 0) ^i)| _∂Ω_i^1=κ_i , 
    . ( ∂ _ξ _iδ u_i -κ_iδ u_i)| _∂Ω_i^1=0 .

Therewith the unperturbed by diffusive interaction solution

    u_( 0) ^i ( ξ_i ) = J_0^i  ξ_i^-1

obeys inhomogeneous boundary condition (<ref>). Here J_0^i  is the i-th Collins-Kimball rate (<ref>).  

One important point to emphasis is that Eqs. (<ref>) and (<ref>) are global and local representations of the same function defined on the manifold:  u_N: Ω^- →(0,1].


Step 3. It can be shown that tensors fields introduced by Eqs. (<ref>) and (<ref>) are families of irregular  {𝐗_n^-(ξ_i  ) }_n=0^∞ and regular {𝐗_n^+(ξ_i  ) }_n=0^∞ solid Cartesian harmonics, which form desired basis functions (a complete and orthogonal sets) in domains Ω^-_i and Ω^+_i, respectively.
Thus, we can represent a partial solution u_i (ξ_i  ) with the help of the irregular ICT, satisfying regularity condition at infinity 

    u_i (ξ_i  ) = ∑_n=0^∞ξ _i^-( 2n+1) 
    	A_γ _1...γ _n^i ⊙^n  ξ
    		_γ _1^i…ξ_γ _n^i^[-.04in].01in.05in
    		.5in.01in[-0.04in].01in
    		.05in    Ω^-_i  .


On the other hand, in an exterior neighborhood of i-th sink boundary ∂Ω_i we can represent the perturbation function in Eq. (<ref>) as a series with respect to the regular ICT 

    ∑_j( ≠i) =1^N u_j (ξ_j  ) = ∑_n=0^∞ B_γ _1 …γ _n^i ⊙^n ξ_γ _1^i…ξ_γ _n^i^[-.04in].01in.05in
    		.47in.01in[-0.04in].01in .05in .

In Eqs. (<ref>) and (<ref>) A_γ _1 …γ _n^i and B_γ _1 …γ _n^i are tensor coefficients to be determined from the boundary conditions (<ref>).

Step 4. Theorem <ref> in the dimensionless form (<ref>) gives the required translation addition theorem to tackle the diffusion problem (<ref>)-(<ref>). 

Step 5. Now, by means of dimensionless translation addition theorem (<ref>), we satisfy the Robin boundary conditions (<ref>). Hence,  taking advantage of the ICT linear independence, we can express the coefficients A_γ _1 …γ _n^i in terms of B_γ_1 …γ _n^i 

    A_0^i=J_0^i (1-B_0^i)       n=0   , 
    
    	A_γ _1…γ _n^i=n-κ_i/1+n+κ_iB_γ _1 …γ _n^i
    	     n ∈ℕ  .


With the aid of these connections formula (<ref>) for the perturbed diffusion field in an exterior neighborhood of the surface ∂Ω_i (i = 1, N) yields

    δ u_i (ξ_i  )  = ∑_k=0^∞[ 1+( k-κ_i/1+k+κ_i) ξ
    	_i^-( 2k+1)
    	] B_γ
    		_1 …γ _k^i ⊙^k ξ_γ _1^i…ξ_γ _k^i^
    		[-.04in].01in.05in.47in.01in[-0.04in].01in
    		.05in  .

One can can easily verify that this solution obeys homogeneous partially reflecting boundary condition (<ref>).

And, in addition, using Eq. (<ref>) under relationships (<ref>) and (<ref>) for domaind Ω^-_j we immediately write down an auxiliary relation

    ∑_j( ≠i) =1^N u_j (ξ_j  ) = ∑_j ( ≠i) =1^N u_( 0) ^j ( ξ_j ) 
    	
    	+∑_j ( ≠i) =1^N 
    	∑_k=0^∞( k-κ_j/1+k+κ_j) ξ _j^-(
    		2k+1) B_γ
    		_1 …γ _k^j ⊙^kξ_γ _1^j…ξ_γ _k^j^
    		[-.04in].01in.05in.45in.01in[-0.04in].01in
    		.05in .

It is noteworthy that all terms here satisfy the regularity condition at infinity (<ref>) as opposed to the local expression (<ref>).

Step 6. Now in the local Cartesian coordinates { O_i; ξ_i} one can carry out a self-consistent procedure to find unknown tensor coefficients B_γ _1 …γ _k^i appearing in formulas (<ref>) and (<ref>).
First we apply the dimensionless translation addition theorem (<ref>) including its degenerate form (<ref>) to the right-hand side of Eq. (<ref>) in order to recast it in the local coordinates { O_i; ξ_i}. Then with allowance made for expression (<ref>) and orthogonality property of the ICT on the unit spheres ∂Ω_i^1 (<ref>) we derive the required self-consistent infinite system of linear algebraic equations (ISLAE) of the II kind with respect to unknown tensor coefficients B_γ _1 …γ _k^i. We should note that these infinite systems are commonly referred to as the resolving ISLAE <cit.>.

According to the statement of the diffusion problem given in Sec. <ref> we should distinguish three particular cases of the resolving ISLAE. 

(a) For partialy reflecting sinks (when 0 < κ_i < ∞) the resolving ISLAE is

    B_γ _1 …γ _k^i = B_γ _1 …γ _k^i(
    		0) + ∑_j ( ≠ i ) =1^N ∑_l=0^∞( k-κ_j/1+k+κ_j
    	) 
    ×U^ij_γ _1 …γ _kμ
    		_1  …μ _l⊙^l  B_μ _1 …μ _l^j  ,     k = 0, ∞ .

Hereinafter for short we denoted

    B_γ _1 …γ _k^i( 0) = ∑_j ( ≠ i ) =1^N J
    	_0^j U^ij_γ
    		_1  …γ _k .

Another two important ISLAE are special cases of (<ref>).

(b) In case of fully absorbing sinks (when κ_i→∞) the resolving ISLAE is simplified to

    B_γ _1 …γ _k^i = B_γ _1 …γ _k^i(
    		0) - ∑_j ( ≠ i ) =1^N ∑_l=0^∞U^ij_γ _1 …γ _kμ
    		_1  …μ _l⊙^l  B_μ _1 …μ _l^j  ,

where k = 0, ∞.

(c) For fully reflecting (when κ_i→ 0) the resolving ISLAE (<ref>) obviously yields

    B_γ _1 …γ _k^i =  ( k/1+k)∑_j ( ≠ i ) =1^N ∑_l=0^∞U^ij_γ _1 …γ _kμ
    		_1  …μ _l⊙^l  B_μ _1 …μ _l^j  ,

where k = 1, ∞.

Step 7. Consider here the most general case of arrays with partialy reflecting sinks. In other words to find B_γ _1 …γ _k^i the system (<ref>) requires the inversion of an corresponding infinite-dimensional matrix <cit.>. So we deal with a typical problem which may be solved by methods of functional analysis <cit.>. For instance, provided the ISLAE (<ref>) is of normal Poincare-Koch type it implies that the Fredholm-Hilbert alternative holds true <cit.>. Hence, a unique solution of (<ref>) exists and it may be found by the method of reduction to any degree of accuracy <cit.>. The latter means that the resolving ISLAE (<ref>) should be truncated and then inverted obtained finite system numerically to get approximations of the coefficients B_γ _1 …γ _k^i. Under some conditions the resolving ISLAE (<ref>) possesses solution which may be obtained by means of simple iterations <cit.> (see also discussion below in Sec. <ref>).  Obviously the zero and the first iterations for B_0^i are

    B_0^i( 0)= ∑_j ( ≠ i ) =1^N J
    	_0^j U^ij_0 ,   
    	
    	B_0^i(1) = B_0^i(
    		0)-∑_j ( ≠ i ) =1^N J
    	_0^j U^ij_0B_0^j( 0) - ∑_j ( ≠ i ) =1^N∑_l=1^∞ J
    	_0^j  
    	U^ij_μ
    		_1  …μ _l⊙^l  B_μ _1 …μ _l^j(
    		0) .

Here and below U^ij_0 is the monopole approximation to the  mixed-basis matrix.	
Interestingly that in the course of solving the ISLAE (<ref>) by simple iterations angular dependence arises starting from terms 𝒪( ε ^4 ).

Step 8. Once the local concentration u_N (P) is found one can deduce the total flux of diffusing particles B onto the i-th sink. At that one should utilize the local representation (<ref>) of function u_N (P). Due to orthogonality condition for the ICT (<ref>) explicit surface integration on the unit sphere in the general formula (<ref>) will yield the desired rates (<ref>).

	Quite apparently, functions ξ̂_γ _1^i…ξ̂_γ _n^i^
		[-.04in].01in.05in.47in.01in[-0.04in].01in
		.05in are harmonics polynomianls defined on the i-th unit sphere ∂Ω_i^1 and, in accordance with the Weierstrass approximation theorem, 
	any function w_i ( ξ̂_i)∈ C(∂Ω_i^1) may be expanded in the series 
	
    w_i (ξ̂_i)
    		= ∑_n=0^∞α _γ _1 …γ _n^i ⊙^n ξ̂_γ _1^i…ξ̂_γ _n^i^
    			[-.04in].01in.05in.45in.01in[-0.04in].01in
    			.05in ,

	where α _γ _1 …γ _n^i are some constant tensor coefficients to be determined by means of the orthogonality condition (<ref>).  One can see that expression (<ref>) makes it possible to apply the ICT method being considered to the diffusion problems under general Robin boundary conditions (<ref>).  	 






 §.§ Main results obtained by using of the ICT technique



Within this subsection we present our main results, which are conveniently to formulate in a form of two theorems. 

Thus, for an arbitrary sinks microstructure X^(N) and reactivities set κ^(N) we have proved the following assertions.

(a) On global representation of the solution u_N (P ).

	Using the local Cartesian coordinates { O_j; ξ_j} (j=1,N) the solution u_N (P ) (P ∈Ω ^-, see Fig. <ref>) to the exterior Robin boundary value problem (<ref>)-(<ref>) may be repesented as follows:
	
    u_N (P ) = ∑_j =1^N    [  J_0^j (1-B_0^j)ξ_j^-1 + 
    		∑_k=0^∞( k-κ_j/1+k+κ_j)  . 
       .   ×ξ _j^-(
    			2k+1) B_γ
    			_1 …γ _k^j ⊙^kξ_γ _1^j…ξ_γ _k^j^
    			[-.04in].01in.05in.45in.01in[-0.04in].01in
    			.05in]  ,

	where tensor coefficients B_γ _1 …γ _k^i (k = 0, ∞) are solution of the resolving ISLAE (<ref>).


(b) On local representation of the solution u_N (P ).

	Solution u_N (ξ_i ) to the exterior Robin boundary value problem (<ref>)-(<ref>) in a small exterior ϵ-neighborhood of the unit sphere ∂Ω_i^1 with respect to the local Cartesian coordinates { O_i; ξ_i} reads
	
    u_N (ξ_i ) = J_0^i ξ_i^-1 + ∑_k=0^∞[ 1+( k-κ_i/1+k+κ_i) ξ
    		_i^-( 2k+1)
    		]  
    ×  B_γ
    			_1 …γ _k^i ⊙^k ξ_γ _1^i…ξ_γ _k^i^
    			[-.04in].01in.05in.47in.01in[-0.04in].01in
    			.05in    Ω ^-_i(ϵ)   i = 1, N ;
 
	where tensor coefficients B_γ _1 …γ _k^i (k = 0, ∞) are solution of the resolving ISLAE (<ref>).

So, if we know the solution u_N (ξ_i ) in an exterior neighborhood adjacent to the i-th reaction surface (<ref>), we can calculate the B-particles flux on the i-th sink and, therefore, the corresponding absorbing rate and, in this way, screening coefficient (<ref>). Thus, taking advantage of Theorem <ref>, one can readily prove an important


	For the exterior Robin boundary value problem (<ref>)-(<ref>) the screening coefficient in Eq. (<ref>) has the form   
	
    J_i = 1 - B_0^i      0 ≤ B_0^i <1  ,

	where coefficients B_0^i  are determined by means of solution to the resolving ISLAE (<ref>).


It is worth noting here that in Ref. <cit.> Krasovitov thoroughly reproducing the notation, derivations and even misprints of our paper <cit.>. However, he has been the first who kept all terms comprising B_γ _1...γ _n^i (n ∈ℕ) in the series with respect to orthogonal system of ICT after integration over the unit sphere (<ref>) (see detailed discussion of these points in Ref. <cit.>).






 §.§ Monopole approximation



As we already noted above in Sec. <ref> the monopole approximation is the simplest tool to study the problem under consideration and, therefore, the most commonly used in various applications. Here we dwell on this approximation briefly.

It is commonly accepted the following

	Monopole approximation corresponds to truncations of series (<ref>) and (<ref>) for any partial solutions u_i (ξ_i  ) (i=1,N) such that: A_0^i,  B_0^i ≠ 0  and  A_γ _1 …γ _k^i, B_γ _1 …γ _k^i  ≡ 0 for all k ∈ℕ.


It is notable that one can see directly from the resolving ISLAE (<ref>) that B_0^i=0 and, therefore, the fully reflecting obstacles cannot be described in the monopole approximation. Besides,  for this case due to Eq. (<ref>) we have J_i = 1 as should be from the physical standpoint.

For definiteness sake  below we will consider the case of fully diffusion-controlled reactions (i.e., the case (b) in step 6.) only.

Utilizing system (<ref>) one infers that for the monopole approximation, values of the total flux 
J_i may be determined by the much simpler finite N × N system of equations:

    J_i=1-∑_j =1^N V^ij_0  J_j ,	   i=1,N ;  
    V^ij_0 := {[    U^ij_0 j ≠ i   ,;         0   i=j   . ].

Therefore, expressions (<ref>) and (<ref>) are simplified to 

    u_N=u_( 0) ^i+B_0^i( 1-1/ξ_i)  ,  
    
    	J _i = 1-B_0^i  .

In the important particular case at N=2 system (<ref>) gives 

    J_i=1- ε_ij J_j ,    (i ≠ j)=1, 2

with the first iteration

    J_i=1- ε_ij ,    (i ≠ j)=1, 2 .

It is evident that system (<ref>) has the exact solution

    J_i= 1-ε_ij/1-ε_jiε_ij , 	   (i ≠ j)=1, 2 .

Note in passing that sometimes in literature Eq. (<ref>) is called the monopole approximation, however, the monopole approximation is given by Eq. (<ref>). Besides Eq. (<ref>) may be also obtained by the reflection method for the monopole approximation <cit.>.
























 §.§ General remarks on the obtained solution



In Sec. <ref> we already stated that exact solution to the diffusion boundary value problem (<ref>)-(<ref>) is possible only for two fully absorbing sinks by means of bispherical coordinates <cit.>. Using the bispherical coordinates for the problem under partially reflecting boundary conditions (<ref>) leads to rather cumbersome recurrence relations and therefore does not provide any advantages compared to other methods. On the contrary, for our approach the kind of boundary conditions does not matter. Moreover, the calculation of the diffusion field u_N (P)  from the boundary value problem (<ref>)-(<ref>) in general case of N sinks microstructure (<ref>), when N >2 is a far too complicated task to be solved analytically in closed form. 


A comparison of the general resolving ISLAE (<ref>) with that (<ref>) for fully absorbing sinks shows the weakening of the diffusive interaction effects in arrays comprising the partialy reflecting sinks. Clearly the diffusive interaction vanishes for the fully reflecting obstacles. 

An N-sink calculation reveals that the pure diffusion field of B-particle around any i-th sink u_( 0) ^i ( ξ_i ) (<ref>) is influenced by a contribution mediated by all other N-1 sinks (<ref>), leading, in it turn, to the screening effects similar to that in an electrostatics <cit.>.


We emphasize again that discussed in Sec. <ref> mixed diffusion boundary value problems with three kinds of obstacles may be treated straightforwardly by the GMSV in terms of the ICT.
We will bring a pattern of numerical calculations of the screening coefficients J_i (<ref>) for the Dirichlet-Robin proper mixed problem at N=3 below in Sec. <ref>.


The GMSV in terms of the ICT may be applied to catalytically activated diffusion-influenced trapping reactions occuring in host media comprising arrays of both finitely many sinks and immobile catalytic domains <cit.>. 

Finally, it is worth noting that stationary diffusion equation under so-called conjugate boundary conditions <cit.> may be also tackled with above method quite similarly.





§ CONVERGENCE ANALYSIS OF THE GMSV



Thus, in principle, we can describe theoretically the trapping of B-reactants for any finite number N of spherical sinks. However, for some sink arrays one faces the difficulties with convergence of iterative solution to the resolving ISLAE (<ref>). Alternative numerical solution of the aforementioned truncated system of equations corresponding to (<ref>) also leads to serious difficulties. 
So, estimations on validity of the GMSV with the help of the ICT for solution of the diffusion problem at issue naturally arises. This question appeared to be rather complicated and, therefore, we discuss it very shortly. Throughout this section we shall deal with the problem (<ref>)-(<ref>) for the case of fully diffusion-controlled reactions.




 §.§ Relations with method of reflections




The point is that according to the weak maximum principle for harmonic functions we have inequality

    max_ P ∈Ω^-| u_N ( P ) |≤max_P ∈∂Ω^-| u_N ( P ) | .

Let a harmonic function v_N (P) be an approximation to the exact solution u_N (P) of the diffusion problem (<ref>)-(<ref>). It immediately follows from inequality (<ref>) that the maximum absolute error of the residual for any approximation | v_N (P) -  u_N (P) | is achieved on the boundary of the configuration manifold ∂Ω^-. So, it is a simple matter to find a priori the upper estimate of the error for the residual of a given approximation in the configuration  manifold Ω^-.

As will become clear below, our consideration concerning convergence of the GMSV procedure it is expedient to start with some results on the method of reflections (in mathematical literature it is also known as the generalized Schwarz’s method) <cit.>.  For two fully absorbing sinks the reflection method leads to a series which reveals slow down convergence as sinks become closer up to the contact, still being convergent (see discussion below). However, even in case of three sinks convergence of the method is not guaranteed at all intersink separations <cit.>. Indeed, studing similar problem on the hydrodynamic interaction in the Stokes flow past three spheres, Kim drew attention to the fact that even in three spheres case convergence of the method of reflections can be violated <cit.>. We do not intend to discuss this subject here and so we refer the interested reader to  Ref. <cit.>. 

According to the method of reflections one shold construct a functional sequence of relevant approximations { v_m( P) }_m=0^∞ so-called reflections (see p. 638 of Ref. <cit.> and  Ref. <cit.> for details) such that this sequence converges in Ω ^- uniformly to the exact solution 

    v_m( P) ⇉ u_N (P)    as    m →∞ .

In Ref. <cit.> we proved the following theorem on the conditions for uniformly convergence of the method of reflections (<ref>).

	For the relationship (<ref>) to hold the following sufficient (a) and necessary (b) conditions must be satisfied
	
    M_N :=  max_i=1,N∑_k( ≠
    			i) =1^N( ε_ik/1-ε_ki) <1  , 
       m_N :=   min_i=1,N∑_k( ≠
    			i) =1^N( ε_ik/1+ε_ki) <1  .


It is important to highlight that currently it is believed these inequalities are still the best ones to estimate convergence of the method of reflections <cit.>.

Plainly, Theorem <ref> implies the following two assertions.

	For sinks of identical radius R_i =R (i= 1,N) conditions (<ref>) and (<ref>)  are simplified to 
	
    M_N =   max_i=1,N∑_k( ≠
    			i) =1^N( L_ik/R -1 )^-1 < 1  , 
       m_N =  min_i=1,N∑_k( ≠
    			i) =1^N( L_ik/R + 1 )^-1 < 1  .



	For the monopole approximation conditions (<ref>) and (<ref>)  are reduced to 
	
    M_N = max_i=1,N∑_k( ≠
    			i) =1^Nε_ik < 1   ,        m_N =  min_i=1,N∑_k( ≠
    			i) =1^Nε_ik < 1 .


One can see that above condition (a) (<ref>)
coincides with well-known sufficient condition for the convergence of the iterative solution of the system (<ref>): 

    ‖V^ij_0 ‖ _∞:= max_i=1,N∑_k=1^NV^ik_0<1 ,

where V^ij_0 is the hollow matrix (<ref>).

Moreover, previously in Ref. <cit.> we have shown that this approach is tightly connected with the corresponding iterative solution within the scope of the GMSV and proved

	Solution to the resolving ISLAE (<ref>) by simple-iteration method is equivalent to the method of reflections.

In this way Theorem <ref> immediately leads to the following remarkable result

	Conditions (<ref>) and (<ref>) may be used to estimate the convergence of the simple-iteration method when applied to the resolving ISLAE (<ref>).


Solvability of the resolving ISLAE by simple-iteration method is very important because of its simplicity and, due to this fact, considerable current use for numerical implementation of the approach at issue <cit.>.  





 §.§ Some estimates of the convergence



Estimates of the convergence for symmetric configuration of sinks are the most simple. First consider  two examples of axisymmetric configurations.

Two sinks of different radii.

For two sink arrays we shall always denote L = L_12=L_21.




Consider first two sinks supposing, for definiteness, that R_1 >  R_2. 
In this particular case quite apparently the sufficient condition of convergence (<ref>) is satisfied and it yields inequality M_2 = R_1/ ( 1 - R_2) < 1. Hence we have 
h = L - ( R_1 + R_2) > 0, where h is the distance between sink surfaces.
In its turn, the latter inequality coincides with non-overlapping sinks condition (<ref>).

Then, by Theorem <ref> and Corollary <ref> we arrive at the analog of the known Goluzin theorem <cit.>.

	For any two disjoint fully absorbing sinks the resolving ISLAE (<ref>) may be solved by the method of simple iterations.


To show how one can directly apply the Corollary <ref>, as examples consider arrays of N fully absorbing sinks of identical radius R_i =R (i = 1,N).  

Sinks of equal radii are located along a straight line.

In addition let us, for simplicity sake, assume that all N sinks are equally spaced at spacing a along their line of centers (see Fig. <ref> (a)). We have performed calculations at N=100. It turned out that the simple-iteration method using to solve the resolving ISLAE (<ref>) is not valid when R/a < 4.9, and when R/a < 9.4 its applicability requires additional justification.



A triangular arrangement of equal sinks.

Diffusion-controlled reactions in the three-sink arrays, particularly equal sinks arranged in an equilateral triangle  have been the subject of numerous studies <cit.>. 

So, consider now the array containing three identical sinks arranged so their centers are at the apices of an equilateral triangle with sides of length a (see Fig. <ref> (b)). One can easy verify that sufficient condition for the convergence of the simple-iteration method (<ref>) fails when 1/3 < R/a < 1/2.
This is entirely consistent with the conclusion by Kim in Ref. <cit.>, where were found numeically that for three sphere system relevant reflections series will diverge when 0.46 ≲  R/a < 0.5 and it becomes convergent for R/a < 0.4. Moreover, Kim stated there: "From a rigorous viewpoint, it is not entirely clear at this point whether this is because of an inherent limitation of a method based on the expansion in R/a (method of reflections) or of because some error in the analysis."

For further results concerning the convergence of the simple-iteration method in the problem at issue, interested readers are encouraged to refer to the recent works Refs. <cit.>.





§ ILLUSTRATIVE EXAMPLES



In this section as illustrations of the above GMSV by means of the ICT we numerically evaluated the diffusive interaction for a few simple examples of many-sink arrays.




 §.§ Iso-concentration surfaces for two sinks



To be specific, assume that the sinks are fully absorbing and start from the simplest, and yet, important case of two sinks.




By means of Theorem <ref> the local concentration field may be provided in the entire manifold Ω^-. For this purpose let us introduce the concentration level surfaces: S ( h ) ⊂Ω^- such that  .  u_2 ( P ) | _S = h ∈(0, 1 ). As usually for discretization we perform uniform partition of the concentration range (0, 1 ) with a fixed step value ▵ h, such that 1 correspond to the reaction surfaces and 0 to infinity: h_l = l ▵ h for l = 0,l_0, where l_0=1/▵ h.

Therewith, for concreteness, we carry out numerical calculations assuming that R_1 > R_2. 
It is convenient to consider dependence of the local concentration  u_2 ( P) on two dimensionless geometrical parameters in units of the smaller sink: reduced sink 1 radius R^∗_1:=R_1/R_2 and the sink spacing L^∗=L/R_2.

Calculated cross sections through the appropriate level surfaces have been plotted in the order of decreasing of the hight h. In  Fig. <ref> and (a) of Fig. <ref> we used  decrement step ▵ h = 0.05, whereas calculations in plot (b) of Fig. <ref> were performed at ▵ h = 0.01.

Analysis of above numerical results for two fully absorbing sinks allows us to formulate general

	There exists a number h_e ∈(0, 1 ) such that one can treat the domain 
	Ω _e:={ P: u_N( P)
	< h_e} as an effective fully absorbing sink with respect to the scaled concentration u^∗_N( P):=  u_N( P)  h_e given on the configuration manifold Ω^- _e:={ P: u^∗_N( P)
	> 1}.


For instance, at N=2 Fig. <ref> shows that one can take h_e = 0.95 to get snowman-shaped effective fully absorbing sink with the configuration manifold Ω_e^-= { P: u_2( P) 
> 0.95}.





 §.§ Calculations of the screening coefficients



Within this subsection to calculate the screening coefficients due to siffusive interaction we shall use the Corollary <ref> of the Theorem <ref> on local representation of the solution.

For two fully absorbing sinks of different sizes numerical estimations show that the corresponding screening coefficients obey inequality J_1 ≥ J_2 at all sink-sink center separations L > R_1 + R_2, wherein J_1 → 1 as R_2 → 0.  Furthermore, there is no diffusive interaction between sinks and any number N_p (1 ≤ N_p ≤ N) of isolated points. Really for the isolated point sinks R_j = ε_ij=0 and, therefore, corresponding mixed-basis matrix elements in the resolving ISLAE (<ref>) are equal to zero.


Two equal fully absorbing sinks.

Again assume that both sinks are fully absorbing.
We chose this case taking into account the fact that iteration solution of the resolving ISLAE (<ref>) converges up to the contact of sinks at N=2 <cit.>.

For simplicity sake consider here two equal sinks R_1=R_2=R, so we have: ε_21  = ε_12 = ϵ∈ (0, 0.5) and for screening coefficients J_1(ϵ)=J_2(ϵ)=J(ϵ). Solving the ISLAE (<ref>) by simple iterations as ϵ→ 0+ we find:

    J(ϵ) = J_a(ϵ) +  O(ϵ^6)  , 
     
    	J_a(ϵ) = 1-ϵ + ϵ^2 - ϵ^3 + 2ϵ^4 - 3 ϵ^5  .
 
It is noteworthy that angular dependence in (<ref>) arises starting from the terms O(ϵ^4). 

To show that J_a(ϵ) very nearly approximates to the exact screening coefficient we need error analysis.  we quantify the difference between the solutions of these test problems. 
For this purpose let us estimate the above approximation relative error bound defined as

    δ_m := sup_ϵ∈ (0, 0.5)δ(ϵ) =  lim_ϵ→ 0.5-δ(ϵ) ,

where δ(ϵ):=  | J (ϵ) - J_a(ϵ) | /J (ϵ) is the relative error.

In the case under consideration with the aid of bispherical coordinate system Samson and Deutch obtained the exact reaction rate <cit.>

    J(ϵ)=∑_n=1^∞( -1) ^1+nsinh( μ_0 ) /sinh( nμ_0 ) ,

where 
μ_0 :=cosh ^-1(  1/2ϵ).
For the contact, i.e. when ϵ→ 0.5- (or μ_0 → 0+) we have fully absorbing dumbbell system and formula (<ref>) yields the exact limit

    lim_ϵ→ 0.5-J(ϵ) =ln 2
 and, therefore, we get desired upper bound for the relative error

    δ(ϵ) < δ_m = 1- 21/32 ln 2≈ 0.053 .

One can see that relative error decreases rather rapidly when ϵ decreases, e.g., δ(0.4 ) ≈ 0.015, δ(0.3 ) ≈ 0.003 and besides J (ϵ) > J_a(ϵ) for all ϵ∈ (0, 0.5).




Fully absorbing sink and two partially reflecting sinks of equal radii.

Finally, we also give a numerical example for screening coefficients (see Subsection <ref>) in the case of three sinks of the same radii (R_i = R, i = 1,3) equally spaced at spacing a along their line of centers. Since the sinks numbering system is arbitrary we number sinks from left to right (see Fig. <ref>). Moreover we suppose that sinks possess different reactivities, namely, sink 2 is fully absorbing (κ_2 →∞) while the other two sinks are of the same, but varied reactivity (κ_1 = κ_3 =κ∈[0, ∞)). Thus, the reactivities set is simplified to κ^(3):={κ, +∞, κ}. Introduce the notations
J(κ):= J_1(κ^(3))=J_3(κ^(3)) and J_2(κ):= J_2(κ^(3)) and assume that a/R = 3. 

Fig. <ref> shows that in both sinks 1 (3) and sink 2 cases influence of diffusive interaction slighty decreases with increasing intrinsic rate constant κ. Furthermore, screening coefficient J(κ) (J_2(κ)) is monotonically increasing  (decreasing)
approaches corresponding horizontal asyymptotes 

    lim_κ→ +∞ J(κ) ≈ 0.70   ,   lim_κ→ +∞ J_2(κ) ≈ 0.55  .

In Fig. <ref> short-dashed straight line corresponds to the horizontal asymptotes (<ref>) for sinks 1 (3) and sink 2, repectively.


For the fully reflecting obstacles limit we have (see corresponding short-dashed straight line)

    lim_κ→ 0+  J(κ) = 0   ,   lim_κ→ 0+  J_2(κ) ≈ 0.97  .

Clearly, since sink 2 is hidden from both sides by neighbors, even inert obstacles reduce the screening coefficient.

Calculations were carried out with by means of the reduction method <cit.>. Note in passing that, although for the microstructure under consideration when κ→∞ relation (<ref>) gives M_3 = 1,  the simple-iteration method also converges.





§ DISCUSSION



This section contains discussion of some distinctive characteristics of the method under consideration and some important conclusions of its application for description of the diffusive interaction between sinks. 

First point out that regardless of its technical implementation the GMSV belongs to the class of methods, when corresponding approximate solution satisfies the diffusion equation (<ref>) and regularity at infinity (<ref>), while the boundary conditions (<ref>) are obeyed only approximately. This fact is of primary concern to theoretical modelling of the diffusion-influenced reactions in an arbitrary ensembles of N sinks. Here, it is worth noting that, e.g., performing solution of the exterior elliptic boundary value problems in unbounded domains by grid methods one faces serious difficulties concerning condition at infinity.







 §.§ Advantages of the ICT technique



We have solved the problem on diffusive interaction between spherical sinks by means of local Cartesian coordinate systems, rather than by introducing polar spherical coordinates. As it might be appeared at first sight, Cartesian coordinates serve well to solve diffusion problems for rectangular geometries, and not for the problems containing spherical domains.
However, it turns out that above assertion generally is not a case. In connection with this it is significant that Hinsen and Felderhof claimed the following: "... the spherical harmonics have disadvantages in numerical calculations." and below they continued: "in both analytical and numerical calculations the Cartesian moments are often more convenient" <cit.>.

It appears that applying of irreducible Cartesian tensors technique leads to rather fast convergence of the local concentration expansion (<ref>). Besides, this approach, contrary to standard one based on the irregular solid spherical harmonics ψ^-_n ( r_i, θ_i , ϕ_i ) <cit.>, allows us to perform all numerical calculations for the real-valued functions only. Note that real-valued solid harmonics have a number of advantages over often used the complex-valued ones.
The choice of irregular irreducible Cartesian tensors X_n^-(𝐫_i ) (<ref>) as basis of functions for solutions of the exterior steady-state diffusion problems in the configuration manifold Ω^-, implies a very rapid convergence of the numerical solution, as confirmed here with several examples. In addition, contrary to the spherical solid harmonics, there are relatively simple recurcion formulas to calculate analytically the irreducible Cartesian tensors (see Appendix <ref>). Thus, our result is in complete agreement with above Hinsen-Felderhof statement.

Moreover, from pure mathematical viewpoint curvilinear coordinates have one more serious drawback.
The preferred usage of coordinates corresponding to the sink surface is a matter of common knowledge <cit.>. In other words, for a given i-th sink Ω_i it is expedient to choose an appropriate local curvilinear coordinates { O_i; ζ}⊂ℝ^3, where ζ = (ζ_1, ζ_2, ζ_3). Therefore, some coordinate surface of this coordinate system, say at a fixed index m_i ∈1,3 there exists a real constant c_i such that the surface of i-th sink reads:
∂Ω_i = {ζ∈ℝ^3 : ζ_m_i= c_i}.
It is clear that in the Cartesian coordinate system the spherical geometry of the boundary leads to a rather complex formulation of the boundary conditions. So, it seems quite reasonable to use appropriate local spherical coordinates { O_i; r_i, θ_i , ϕ_i } to describe diffusion of B-reactants towards i-th spherical sink in domain Ω_i^-:= ℝ^3\Ω_i. However, the regularity condition for coordinates to be curvilinear is broken at some points. Particularly, in above case of polar spherical coordinates to have one-to-one property in Ω_i^-, Ox_3 axis and half-plane { x_1 ≥ 0, x_2 =0 } should be excluded <cit.>. A number of other advantages of the ICT technique are described also elsewhere <cit.>.


Finally note that from the analytical point of view, obtained in Sec. <ref> solution u_N ( P ) keeps an explicit analytical dependence on P. In turn, from the numerical standpoint, this property allows one to truncate the resolving ISLAE at lower sizes as compared to other numerical methods such as, e.g., a finite elements method <cit.>. 






 §.§ Description of the diffusive interaction



Let us dwell now on several particularly important points on the diffusive interaction which directly follows from our results.

If an array comprises only fully reflecting obstacles we have no diffusive interaction between them. On the other hand in situation when a given array contains at least one absorbing sink the diffusive interaction exists and we should also take into account influence of inert obstacles. 

It is very important to note that here the ICT technique was applied to the diffusion problems under consideration posed in the unbounded domains, whereas the same approach is also applicable to the similar interior boundary value problems describing reactions occurring inside a spherical region (see Refs. <cit.> for details). So, the GMSV in terms of ICT can be very useful not only in the theory of diffusion-influenced reactions, but also in biology, when reactions can take place into bounded domains <cit.>. We note, incidentally, that Smoluchowski's rate constant (<ref>) cannot describe even one-sink reaction inside a larger spherical domain since it does not consider diffusive interaction effect due to the boundary of englobed domain.

In Ref. <cit.> authors concluded that their system of equations (3.97) solves the N-sink diffusion problem under Smoluchowski boundary conditions exactly. However, results obtained in  Sec. <ref> show that the system (3.97) has been derived for the monopole approximation and, therefore, the above conclusion is incorrect.

It is interesting to note that, clustering of sinks was found to reduce the absorbing rate even more significantly compared with nonclustered systems <cit.>. The method gives us possibility to tackle trapping reactions in several clusters of obstacles and sinks in the same way as for separated sinks.

Taking into account our previous results <cit.> and Proposition <ref> we can conclude that diffusive interaction "glues" fully absorbing sinks effectively onto a fully absorbing effective sink. Therefore the whole N-sink array behaves like a single, isolated fully absorbing sink of some characteristic size. 




For a number of applied tasks, one should use the derived microscopic solution (<ref>) (or (<ref>)) along with microscopic trapping rates (<ref>) calculating appropriate configurationally averaged values to obtain effective macroscopic values associated with heterogeneous media comprising different types of sinks <cit.>.

Finally, let us dwell on one more possible application of the approach at issue. In Ref. <cit.> we have established for the first time that even classical passive phoresis of microparticles (i.e., motion due to a constant gradient of some scalar field) should be treated by means of their interaction with so-called surrounding bodies, which are responcible for the external gradient. In this context, it is worth to note that the first and only attempt to apply a similar approach for studing particle interactions for thermophoresis was undertaken by Keh and Chen in Ref. <cit.>. They used so-called polyadic tensors but in fact, these tensors are the same as the ICT up to a constant factors. Again the GMSV in terms of ICT can be used to investigate these quite a challenging problems.





§ CONCLUSIONS AND PERSPECTIVES
























































Motivated through new applications in physics, chemistry and biology this paper deals with the generalized method of separation of variables by means of the irreducible Cartesian tensors formalism extending our previous study on the diffusion-controlled reactions to the general case of diffusion-influenced reactions. We have developed here a  complete modification of the above mentioned method which essentially hinges on a number useful properties of the irreducible Cartesian tensors. Formally we reduced the steady-stae diffusion problem to a set of N Laplace equations with respect to partial solution together with the prescribed N boundary conditions.

In order to determine the place and significance of our research we have presented a rather comprehensive review of literature devoted to the theoretical methods for diffusion-influenced reactions on spherical sinks.
Besides, we paid close attention to the mathematical aspects of the appropriate diffusion problem, namely: (a) rigorous statement of the problem, including terminology clarification, (b) refining solution algorithm and (c) justification of the method convergence. The formulation of the geometric part of the diffusion boundary value problems was performed the most fully than anywhere else. We highlighted the importance of a concept of a smooth manifold for the future theoretical investigations on diffusion-influenced processes occuring in domains with multiply connected boundaries.  We have given a proof of the irregular to regular translation addition theorem in terms of the irreducible Cartesian tensors which seems to be the simplest among previously known theorem for solid harmonics in terms of polar spherical coordinates. For the first time corresponding solution algorithm was divided into eight clear steps. Then we implemented this algorithm step by step in full detail to solve the posed diffusion problem. 
A convergence analysis was conducted to verify the possibility of simple-iteration method usage in order to solve appropriate resolving infinite system of linear algebraic equations. It turns out that the method of reduction may be used to solve numerically the resolving system to any degree of accuracy even when the simple-iteration method is inapplicable. Wherein the desired accuracy of the solution can be achieved just by increasing the truncation degree.
Furthermore, note that the method has advantages such as shorter calculation time and more accurate results compared to purely numerical methods. 

Our results fully confirm known Hinsen-Felderhof statement that it is more preferred to work with harmonic functions using their irreducible Cartesian tensor form throughout diffusion-reaction problems even for ensembles of spherical sinks.

This research provides also a better understanding of the reactants diffusion among arrays comprising arbitrary collection of arbitrary many fully absorbing, reflecting and partially reflecting sinks from both mathematical and physical viewpoints.

The approach developed here is applicable to a wide class of reaction-diffusion systems, providing useful tests for numerical calculations along with approximate analytical estimates of the local concentration field and reaction rates in many specific cases of practical importance. 



Finally it is important to highlight that numerical methods can handle diffusion problems for complex geometries in unbounded domains, yet they require to introduce an articial boundary, reducing the original exterior domain to a bounded domain. The generalized method of separation of variables does not have this drawback and provides simple direct numerical and analytical calculations.
So we can obtain approximations as the truncated point multipole expansions up to the monopole, dipole, quadrupole etc. order. 

It is essential to note that self-phoresis of the Janus particles moving due to a chemical	reaction in the surrounding fluid is investigated by means of monopole approximation <cit.>.
Using the approach considered here this important problem may be attack taking into account multipole corrections <cit.>.

The last but not least: in the present study we clearly show that the theory of 3D differentiable manifolds is an appropriate mathematical tool to describe diffusion and reactions in domains with multiply connected boundaries. In this context, we note that, rigorously posing and solving the appropriate diffusion problems, we, nevertheless, endeavored to avoid considering a number of inherent subtle mathematical points. This issues will be a topic for the future work. 








  
  


§ APPENDIX. BACKGROUND ON THE IRREDUCIBLE CARTESIAN TENSORS




To make this paper maximally self-contained and facilitate applications in diffusion-influenced reactions theory, we provide here necessary mathematical notations, definitions and facts, required in rigorous formulation and applications of the ICT technique.




 §.§ Definitions and examples




It is well-known that there exits two alternative approaches to harmonic functions theory based on either Laplace's representation of solid harmonics with respect to the polar spherical coordinates or Maxwell's theory of poles given in terms of the successive derivatives of fundamental solution written for a Cartesian coordinate system <cit.>. The latter approach may be formulated with the aid of the irreducible Cartesian tensors.

Generally speaking, the irreducible tensors transform according to irreducible representation of the special rotation group  SO( 3). In other words the irreducible tensors like spherical harmonics form a basis in the representation space of group SO( 3) (see <ref>).	
To our knowledge, the ICT simple mathematical description suitable for applied physicists has been provided in Ref. <cit.>. Moreover, for many years there has been a significant number of publications devoting to the use of the ICT formalizm in a wide variety of applications (see, e.g., Refs. <cit.> and references therein).
As regards adopted notations and definitions, in this paper we basically follow known books by Hess <cit.> and Snider <cit.>, where  comprehensive discussions on the various facets of irreducible Cartesian tensors are presented.

A real tensor field of n-th order defined on ℝ^3 we denote by T_n (𝐫) in
coordinate-free notation. The same tensor field can be also described by its components with respect to a fixed orthonormal basis  {𝐞_α}_α = 1^3 given in { O; 𝐫}, i.e. T_n (𝐞_γ _1, … , 𝐞_γ _n)  = T_γ _1…γ _n(𝐫). Here and in the following, T_γ _1…γ _n(𝐫) are 3^n components of n-rank Cartesian tensor field T_n (𝐫). It is clear that the action of the nabla operator with ∂ _r_μ := ∇ _μ
to a field, which is a tensor of rank l, yields a tensor of rank l+1. We recall in passing that one should not distinguish between co- and contravariant index due to duality of Cartesian basis <cit.>.

The product symbol ⊙^q represents a full tensorial contraction of multiplicity q or shortly the q-fold contraction between two tensors with the convention S_0 ⊙ T_0 :=   T_0. Particularly, for tensors S_n and T_m and q ≤ n, m we have a tensor 

    S_n ⊙^q  T_m := ∑_γ _1…∑_γ _q S _γ _1 …γ _n⊙^q T_μ _1 …μ _m

of rank n+m-2q, whereas it is a scalar if n=m=q.

Any tensor, which cannot be reduced to a tensor of lower rank is known as an irreducible tensor <cit.>.
While a general (reducible) n-th order Cartesian tensor T_n has 3^n components, corresponding irreducible tensor besides the rank n is likewise characterized by its weight j <cit.>. It is common knowledge that latter tensor has 2j+1 independent components, which form the basis of the j weight irreducible representation of the full 3D special rotation group SO ( 3) including mirror <cit.>. 

Furthermore, among all ICT obtained from a given n rank tensor by means of reductions there exists one unique irreducible tensor with the same weight and rank n=j. It is known in literature as ICT in natural form <cit.>. We deal with irreducible tensors in natural form because they best suited to apply within the scope of the GMSV.

Thus, for n rank Cartesian tensors our focus lies on their irreducible parts in natural form at that we use notations T_n ^
	[-.04in].01in.05in.15in.01in[-0.04in].01in
	.05in and T_γ _1…γ _n
^[-.04in].01in.05in.34in.01in
	[-0.04in].01in.05in, where symbol (  ...  )
^[-.04in].01in.05in.25in.01in
	[-0.04in].01in.05in stands for the 
irreducible part of a Cartesian tensor T_n of rank n ≥ 2  with respect to the special rotation group SO(3) <cit.>.


The irreducible tensor in natural form is appeared to be fully symmetric and traceless, therefore, this allows us to use rather simple definition of the irreducible Cartesian tensors <cit.>.

	A Cartesian tensor of rank n ≥ 2 is called irreducible Cartesian tensor, if it is: (i) totally symmetric under an arbitrary permutation of the indices; and (ii) traceless under the contraction of any pair of indices.


An important point is that according to definition (<ref>) the ICT can be expressed in terms of successive partial derivatives of the generating function 1/r for all points in ℝ^3 \{0}.	
 
	An irreducible Cartesian tensors of rank n (when n ≥ 2)  may be defined as <cit.> 
	
    r_γ _1… r_γ _n^
    			[-.04in].01in.05in.45in.01in[-0.04in].01in
    			.05in=( -1) ^n/( 2n-1) !!r^2n+1∂
    		_r_γ _1…∂ _r_γ _n( 1/r)   .

	Here r_γ _ν are Cartesian coordinates of a vector 𝐫 (γ _ν =1,3),
	and ( 2n-1) !!:=1· 3·…·( 2n-1) indicates the double factorial for odd natural numbers with the convention (-1)!!:=1.


	Check that the tensors introduced by Eq. (<ref>) are symmetric for any pair of indexes and traceless.
	
	(i)  Indeed,  Eq. (<ref>) involves directly that the introduced tensor is symmetric with respect to permutation of its indexes 
	
    r_γ _1…r_α…r_β…r_γ _n^
    			[-.04in].01in.05in1.15in.01in[-0.04in].01in
    			.05in = r_γ _1…r_β…r_α…r_γ _n^
    			[-.04in].01in.05in1.15in.01in[-0.04in].01in
    			.05in .

	(ii) Moreover, due to harmonicity of the generating function 		
	
    ∇^2 ( 1/r)  = δ_αβ⊙^2 ∂_r_α∂_r_β( 1/r)= 0       r>0

	tensors (<ref>) are also traceless with respect to any pair of their indexes
	
    δ_αβ⊙^2  r_γ _1…r_α…r_β…r_γ _n^
    			[-.04in].01in.05in1.15in.01in[-0.04in].01in
    			.05in = 0   ,

	where δ_αβ is the Kronecker delta. Thus Eq. (<ref>) satisfies both points of definition (<ref>).


	Note furthermore that traces cannot be defined for any scalars a_0 ∈ℝ (tensors of rank 0) or vectors (r_1, r_2, r_3) (tensors of rank 1)). Nethertheless they are harmonic and, therefore, irreducible by definition <cit.>.  


In applications are also used the multipole potentials defined by: X_n^∓(𝐫)=X^∓ _γ _1 …γ _n(  r) <cit.>, where components read

    X^- _γ _1 …γ _n(  r):=( -1) ^n∂
    	_r_γ _1…∂ _r_γ _n( 1/r) ,  
    X^+ _γ _1 …γ _n(  r):= ( 2n-1)!!r_γ _1… r_γ _n^
    		[-.04in].01in.05in.45in.01in[-0.04in].01in
    		.05in .

Owing to definition  they are connected by the relation

    X_n^+(𝐫) =  r^(2n+1)X_n^-(𝐫)  .

For reasons explained below (see <ref>), it is quite natural to call potentials X_n^-  and  X_n^+ irregular and regular ICT, respectively <cit.>. However, throughout this paper we use ICT defined by Eq. (<ref>), since they are most extensively employed in applications.

Denote by d^2 r̂ the magnitude of the surface
element on the unit sphere in global coordinates 

    ∂Ω_1 ≡∂Ω( 0; 1 ) := {𝐫∈ℝ^3: r̂=1 } .

Then for all n ≥ 0 we can write the following orthogonality condition <cit.>

    ∮_∂Ω_1r̂_γ _1…r̂_γ _n^
    		[-.04in].01in.05in.47in.01in[-0.04in].01in
    		.05in·r̂_μ _1…r̂_μ _k^
    		[-.04in].01in.05in.49in.01in[-0.04in].01in
    		.05in d^2 r̂ 
    	= 4π n!/( 2n+1 )!!δ_n kΔ^(n) _γ _1 …γ _n, μ
    		_1  …μ _n ,

where in case of sphere 𝐧̂ = 𝐫̂  (in component notation r̂_α=r_α/r) is a unit vector identifying a point on ∂Ω_1 and . r_γ _1…r_γ _n^
	[-.04in].01in.05in.48in.01in[-0.04in].01in
	.05in| _∂Ω_1 = r̂_γ _1…r̂_γ _n^
	[-.04in].01in.05in.47in.01in[-0.04in].01in
	.05in is the restriction of the ICT to the unit sphere (<ref>). In Eq. (<ref>)
Δ^(n) _γ _1 …γ _n, μ_1  …μ _k stands for the 2n rank projector, which project any tensor of rank n into its irreducible part. Particularly, three first projectors are <cit.>

    Δ^(0) = 1 ,   Δ^(1) = δ_αβ , 
    Δ^(2) _αβ, α' β' = 1/2( δ_αα'δ_ββ'+δ_αβ'δ_βα') - 1/3δ_αβδ_α' β' .







 §.§ Efimov's recurrent formula


In <ref> we have already mentioned that an important advantage of the ICT over solid harmonics is relatively simple recurcion formulas obtained by rather easy successive differentiation of generating function. Interested readers are encouraged to refer to works, where explicit recurrent expressions for the ICT are comprehensively presented <cit.>. Here let us dwell briefly on a little-known, but simple and powerful approach proposed by Efimov <cit.> that allowed us to calculate ICT. 


To derive the relevant Efimov's recurrent formula first we introduce an auxiliary tensor field

    Γ _γ
    		_1 ... γ _n(  r):= ( -1) ^nX^+ _γ _1 …γ _n(  r) .


Inversion of the radius vector 𝐫→𝐫 r^-2 leads to an isomorphic  mapping  of the  gradient operator in the inverted space of harmonic functions (with opposite sign), i.e. the vector operator 𝐃=-r^-1∇ _1/r r. Direct calculations yield the expression

    𝐃:=2𝐫( 𝐫·∇)
    	-r^2∇ +𝐫 .

With the help of this operator it may be shown that the following recursion holds 

    Γ _0=1,   Γ _γ _n γ _1...γ _n-1^( n) =D_γ _n·Γ
    _γ _1...γ _n-1^( n-1)      n∈ℕ .

Taking this relations into account one easily derives the desired Efimov's formula <cit.>

    Γ _γ _1 ...γ _n=D_γ _1D_γ _2...D
    	_γ _n· 1  .



Hence, in the explicit form the first several irreducible Cartesian tensors of ranks n corresponding to the monopole (n=0), dipole (n=1), quadrupole (n=2) and octupole (n=3) potentials read as follows:

    a_0^
    		[-.04in].01in.05in.13in.01in[-0.04in].01in
    		.05in := a_0   , 
    	r_γ _1^
    		[-.04in].01in.05in.13in.01in[-0.04in].01in
    		.05in := r_γ _1 , 
    	r_γ _1r_γ_2^
    		[-.04in].01in.05in.3in.01in[-0.04in].01in
    		.05in = r_γ _1r_γ _2  - 1/3r^2 δ_γ _1 γ _2 , 
    r_γ _1r_γ _2r_γ _3^
    		[-.04in].01in.05in.43in.01in[-0.04in].01in
    		.05in = r_γ _1r_γ _2r_γ _3 - 1/5r^2 (r_γ _1δ_γ _2 γ _3+r_γ _2δ_γ _1 γ _3+r_γ _3δ_γ _1 γ _2)   .













 §.§ Connection with solid harmonics



As indicated in Sec. <ref>, solid harmonics were thoroughly used previously to investigate many-sink effects during diffusion-influenced reactions. So, it is expedient to point to an intimate relationship between them and the ICT <cit.>.

Recall that real solid harmonics ψ_nm(r, θ, ϕ), where the integer m such that - n ≤ m ≤ n and n=0,∞, are harmonic functions in spherical coordinates { O; r, θ, ϕ} <cit.>. For brevity we shall denote them as ψ_n(r, θ, ϕ).
Further, according to their behavior at the origin O, two kinds of solid harmonics are distinguished:

(1)  irregular at O solid harmonics ψ_n^-(r, θ, ϕ) has a pole of order n+1 singularity at O;

(2) regular at O solid harmonics ψ_n^+(r, θ, ϕ) are real-valued harmonic homogeneous polynomials of order n.

Hence, irregular solid harmonics ψ_n^- automatically satisfies regularity condition at infinity (<ref>) and, therefore, they are regular at infinity. 

Clearly, irregular X_n^- (<ref>) and regular  X_n^+ (<ref>) ICT correspond to irregular and regular solid harmonics ψ_n^- and ψ_n^+, respectively.

Represent tensors X_n^±(  r) in the form most closely resemble classical irregular ψ_n^-(r, θ, ϕ), regular ψ_n^+(r, θ, ϕ) solid and real spherical Y_n(θ, ϕ) harmonics:

    X^- _γ _1 …γ _n(  r)=r^-(n+1) Y_γ _1…γ _n(𝐫̂)   , 
    X^+ _γ _1 …γ _n(  r)=r^n Y_γ _1…γ _n(𝐫̂)  , 
    
    	Y_γ _1…γ _n(𝐫̂):=( 2n-1)!!r̂_γ _1…r̂_γ _n^
    		[-.04in].01in.05in.45in.01in[-0.04in].01in
    		.05in .

It is apparent that the ICT  Y_γ _1…γ _n(𝐫̂) like spherical harmonics depend on the direction of vector 𝐫 only <cit.>.

As noted above, the ICT similar to solid harmonics form a basis in the representation space of the group SO( 3), therefore, there exists a non-singular linear transformation and its inverse in space ℝ^2n+1 which map these basises into each other <cit.>. Using known Sylvester's theorem <cit.> one can find the desired connection between tensors X_n^± and classical real solid harmonics ψ_n^± explicitly.  We emphasize here that original Sylvester's proof is rather sophisticated <cit.>, therefore, interested readers are encouraged to refer to a simple proof of Sylvester's theorem proposed by Backus in Ref. <cit.>. 

Thus, it can be shown that Sylvester's theorem implies

    1/r^n+1Y_n(θ, ϕ) = ∑_m=0^n∑_l=0^m  c_lm∂_r_1^l∂_r_2^m-l∂_r_3^n-m( 1/r) , 
    where Y_n(θ, ϕ) = ∑_m=0^n  Y_n^m(θ, ϕ)  , 
    
    	Y_n^m(θ, ϕ) :=( a^m_n cos mϕ +  b^m_n sin mϕ)P^m_n (cosθ)

is the real spherical harmonic of degree n and order m.  In Eq. (<ref>) a^m_n, b^m_n and c_lm are some known real coefficients and P^m_n (cosθ) are associated Legendre functions <cit.>.

Particularly, for the axially symmetric case (m=0), taking all derivatives along O r_3 axis, we have r_γ _1 = r_γ _2= … = r_γ _n = r_3 and formula (<ref>) boils down to well-known relation <cit.>

    1/r^n+1 P_n( cosθ) = ( -1) ^n/n!∂^n _r_3( 1/r)  .









