# The generalized method of separation of variables for diffusion-influenced
reactions: Irreducible Cartesian tensors technique

Sergey D. Traytak<EMAIL_ADDRESS>Semenov Federal Research Center for
Chemical Physics, Russian Academy of Sciences, 4 Kosygina St., 119991 Moscow,
Russian Federation

###### Abstract

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.

††preprint: AIP/123-QED

##  I Motivation

It is presently accepted that for the most part, chemical reactions in micro-
heterogeneous liquid media are contact diffusion-influenced Eyring, Lin, and
Lin (1980); Rice (1985); Ovchinnikov, Timashev, and Belyy (1989); Barzykin,
Seki, and Tachiya (2001); Dagdug, Peña, and Pompa-García (2024). 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. III). 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 Traytak (1995, 2008) 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. Traytak (2003). 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. Traytak and Grebenkov (2018); Grebenkov
and Traytak (2019).

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 Rayleigh (1892) (see Sec. V).
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 Traytak (1991,
1992, 2003); Traytak and Grebenkov (2018); Grebenkov and Traytak (2019) 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. III) 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.

##  II 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: $\mathbb{N}$, $\mathbb{Z}$, $\mathbb{R}$ and
$\mathbb{R}_{+}$ 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 ${\mathbb{R}}^{3}$. As is customary, let $\partial\Omega$
denotes the closed boundary of a 3D domain (connected open subset)
$\Omega\subset{\mathbb{R}}^{3}$ such that its closure is
$\overline{\Omega}=\Omega\cup\partial\Omega$.

We shall study real-valued functions $f:\Omega\to\mathbb{R}_{+}$. Let
$m\in\mathbb{N}\cup\\{0\\}$ (particularly $C(\Omega):=C^{0}(\Omega)$) then
$C^{m}(\Omega)$ stands for the class of $m$ times continuously differentiable
functions on $\Omega$ and $C^{m}(\overline{\Omega})$ denotes a space of
functions having $m$ continuous derivatives in $\Omega\subset\mathbb{R}^{3}$
and continuously extendable to the closure $\overline{\Omega}$ of $\Omega$. A
function of class $C^{\infty}(\Omega)$ is an infinitely differentiable
functions. In the following, the term smooth function on $\Omega$ refers to a
$C^{\infty}(\Omega)$-function.

Let $\mathbf{a}=({a}_{1},{a}_{2},{a}_{3})\in\mathbb{R}^{3}$ be a fixed
translation vector (bold symbols are used for vectors or tensors in
${\mathbb{R}}^{3}$). For any point $\mathbf{x}\in\mathbb{R}^{3}$ define the
translation operator (and its inverse one) $T_{\mathbf{a}}$
($T^{-1}_{\mathbf{a}}$)$:\mathbb{R}^{3}\to\mathbb{R}^{3}$ by the relations

$T_{\mathbf{a}}\mathbf{x}=\mathbf{x}+\mathbf{a}\,,\quad
T^{-1}_{\mathbf{a}}\mathbf{x}=\mathbf{x}-\mathbf{a}\,.$ (1)

It is well-known that translations form a subgroup of the 3D Euclidean motion
group $E(3)$ such that: corresponding composition operator
$T_{\mathbf{c}}=T_{\mathbf{a}}T_{\mathbf{b}}$ with translation vector
$\mathbf{c}=\mathbf{a}+\mathbf{b}$; and identity operator
$I:=T_{\mathbf{0}}=T^{-1}_{\mathbf{a}}T_{\mathbf{a}}$ Miller (1977).

###  II.1 A model description

To study diffusion-influenced reactions theoretically so-called trapping model
is often used den Hollander and Weiss (1994). 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 Beenakker and Ross
(1985)), quiescent, isotropic and homogeneous, therefore mathematically it is
well modelled by the 3D Euclidean space ${\mathbb{R}}^{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  Weiss (1986). 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 Rice (1985).

In this paper we deal with irreversible bulk contact diffusion-influenced
reactions described by the simplest reaction scheme Rice (1985); Ovchinnikov,
Timashev, and Belyy (1989)

$\displaystyle A+B\stackrel{{\scriptstyle k}}{{\longrightarrow}}A\cdot
B\stackrel{{\scriptstyle k_{in}}}{{\longrightarrow}}A+Products\,.$ (2)

Here $A\cdot 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\cdot B$ to form inert products, respectively
Rice (1985). Clearly, the rate $k_{in}$ describes surface reactivity of sinks
$A$ Torquato (2002), which we assume to be spherically symmetric.

Within the scope of scheme (2) the macroscopic action mass law equations for
the bulk concentrations of reactants $c_{A}\left(t\right)$ and
$c_{B}\left(t\right)$ are

$\displaystyle
t_{r}{\dot{c}}_{B}\left(t\right)=-c_{B}\left(t\right)\,,\quad{\dot{c}}_{A}\left(t\right)=0\,;$
(3)

where $t_{r}:=\left(kc_{A}\right)^{-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}\left(0\right)\ll c_{B}\left(0\right)$ holds Ovchinnikov, Timashev, and
Belyy (1989). To put it otherwise, reactions (2) are the catalytic type
Bénichou  _et al._ (2005), when activity of the catalyst $A$ remains
unchanged. In this connection, reactions (2) are called sometimes catalytic
degradation reactions.

###  II.2 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 (2). 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
Brailsford (1976); Rice (1985).

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 $\phi_{A}$ in dilute enough sink ensembles the
following condition holds $\phi_{A}:=4\pi c_{A}R^{3}/3\ll 1$ Kipriyanov,
Gopich, and Doktorov (1998). In other words diffusive interaction accepted to
be weak and, therefore, it is commonly assumed that sinks mostly interact with
the unbounded host medium Wagner, Kampmann, and Voorhees (2001).

Smoluchowski’s theory corresponds to the limiting case of infinite dilution
($\phi_{A}\to 0$) Venema and Bedeaux (1989) and, therefore, only one-sink
system should be treated.

Let $\Omega_{1}\subset{\mathbb{R}}^{3}$ be the test sink domain and consider
its compliment
$\Omega_{1}^{-}:=\mathbb{R}^{3}\backslash\overline{\Omega}_{1}$. It is
expedient to use the standard polar spherical coordinates
$\\{O;r,\theta,\phi\\}$ ($\theta$, $\phi$ are the polar and azimuthal angles,
respectively) coordinate systems attached to the origin $O({\mathbf{0}})$
coinciding with the test sink center. Consider a current point
$P\in\Omega_{1}^{-}$ with position vector $\mathbf{r}$ such that
$r=\|{\mathbf{r}}\|$ its radial distance, where $\|\cdot\|$ stands for the
common Euclidean norm. Plainly, for the spherical sink domain and relevant
reaction surface we have: $\Omega_{1}:=\left\\{r<R\right\\}$ and
$\partial\Omega_{1}=\left\\{r=R\right\\}$.

######  Definition II.1.

For structureless point-like $B$ particles the configuration space is a 3D
domain composed of all $B$ particle positions $\mathbf{r}$.

Since diffusion of $B$ particles occur in the exterior of the sink, domain
$\Omega_{1}^{-}$ is the configuration space and clearly $B$ particles local
concentration is spherically-symmetric $n_{1}\left(r\right)$. Often function
$n_{1}\left(r\right)$ is scaled with its bulk concentration $c_{B}$ so that
$n_{1}\left(r\right)/c_{B}$ treated Rice (1985); Eun, Kekenes-Huskey, and
McCammon (2013). 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

$\displaystyle
u_{1}:\Omega_{1}^{-}\to\left(0,1\right)\,,\quad\mbox{where}\quad
u_{1}\left(r\right)=1-\frac{n_{1}\left(r\right)}{c_{B}}\,.$ (4)

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}\left(r\right)$ 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 Rice (1985)

$\displaystyle\nabla^{2}_{r}u_{1}=0\,\quad\mbox{in}\quad\Omega_{1}^{-}\,,$ (5)
$\displaystyle\left.\left[4\pi
R^{2}\,D\partial_{r}u_{1}-{k}_{in}\left(1-u_{1}\right)\right]\right|_{r=R}=0\,,$
(6) $\displaystyle\left.u_{1}\right|_{r\rightarrow\infty}\rightarrow 0\,,$ (7)

where $D$ is the translation diffusion coefficient and
$\nabla_{r}^{2}:=2{r^{-1}}\partial_{r}+{\partial_{r}}^{2}$ is the radial part
of Laplacian. For brevity, henceforth $\partial_{\varsigma}$ stands for the
partial derivative $\partial/\partial{\varsigma}$ with respect to the
independent variable $\varsigma$. Note in passing that in physical literature
condition (7) is commonly referred to as the outer boundary condition Lee and
Karplus (1987); Eun, Kekenes-Huskey, and McCammon (2013).

It has been well-known that solution to this problem is Rice (1985)

$\displaystyle u_{1}\left(r\right)=1-\frac{R}{r}J_{0}^{1}\,.$ (8)

In this formula we introdused dimensionless rate
$J_{0}^{1}:=k_{CK}/k_{S}=\kappa_{1}/\left(1+\kappa_{1}\right)<1$, where
$\kappa_{1}:=k_{in}/k_{S}$. Hereinafter $J_{0}^{1}$ is the Collins-Kimball
rate constant $k_{CK}\left(k_{in}\right)$ normalized by the Smoluchowski rate
constant $k_{S}$ given by

$\displaystyle k_{S}=\lim_{k_{in}\to\infty}k_{CK}=4\pi
RD=\Phi_{1}\left(R\right)/c_{B}\,,$ (9)

where $\Phi_{1}\left(R\right)$ is the total flux of $B$-particles on the test
sink of radius $R$ Rice (1985).

We face up to another important aspect. Lee and Karplus highlighted the fact
that in terms of non-normalized concentration $n_{1}\left(r\right)$ the outer
boundary condition is time dependent Lee and Karplus (1987), i.e.

$\displaystyle\left.n_{1}\right|_{r\rightarrow\infty}\rightarrow
c_{B}\left(t\right)\,.$ (10)

To clarify the matter we note that the steady state reaction regime occurs
under the following condition $t_{D}\ll t\lesssim t_{r}$ Traytak (2013), where
$t_{D}=R^{2}/3D$ being the characteristic encounter time Doktorov (2020) and
the frequency of reaction (3) takes the form $t_{r}^{-1}=4\pi DRJ_{1}c_{A}$.

Strictly speaking, one should consider the above boundary value problem
(5)-(7) as a quasi-steady state one Beenakker and Ross (1985). That is
consistent with the fact that the characteristic time scale of the change in
$c_{B}\left(t\right)$ (see Eq. (3)) is $t_{r}\gg t_{D}=t_{r}J_{1}\phi_{A}$.
Really, one can see that "fast" (microscopic) time evolution relaxation
towards a quasi-steady-state occurs for times $t={\cal{O}}(t_{D})$ while
changing of bulk concentration $c_{B}\left(t\right)$ takes place at "slow"
(macroscopic) times $t={\cal{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.

##  III 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. (9), leads to the total flux of $B$’s
into the whole ensemble of $N$ identical uncharged sinks as the additive sum
$\Phi_{N}=N\Phi_{1}\left(R\right)$ Richards (1986). 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 (2) 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 Frisch and Collins
(1952); Ham (1958); Deutch, Felderhof, and Saxton (1976); McDonald and
Strieder (2004); Bluett and Green (2006); Eun, Kekenes-Huskey, and McCammon
(2013); Biello and Samson (2015); (2) diffusive (diffusional) interaction
Carstens, Williams, and Zung (1970); Beenakker and Ross (1985); Dubinko  _et
al._ (1989); Traytak (1991, 1992); Marsh (1995); Wagner, Kampmann, and
Voorhees (2001); Galanti  _et al._ (2020); Cedeno  _et al._ (2022); (3) effect
of multiparticle interaction Labowsky (1978); Marberry, Ray, and Leung (1984);
Zinchenko (1994); (4) Laplace’s interaction Sangania and Mo (1996); (5)
concentration effects Kim, Lee, and Kim (2014); (6) collective effects Carrier
_et al._ (2016); Stark (2018); Masoud, Howell, and Stone (2021); (7) reactive
interference Lee and Lee (2020); (8) shielding effects Schofield  _et al._
(2020); (9) cooperative effects Masoud, Howell, and Stone (2022) and even (10)
chemical interaction Stark (2018). 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" Baydoun and Zagrebnov (2011)) best describes the essence of the
interaction under consideration for the steady-state diffusion, heat transfer
and Stokes hydrodynamics Kim and Karrila (2013). 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 Galanti
_et al._ (2020).

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 Frisch and Collins
(1952). Among early approaches it is also worth mentioning Ham’s theory based
on the Wigner-Seitz cell model Ham (1958). Applying a variational procedure
Reck and Prager had first established rigorous upper and lower bounds on the
rate of a diffusion-controlled reactions Reck and Prager (1965). Detailed
exposition of this approach with applications can be found in the book Ref.
Torquato (2002).

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 Borzilov and Stepanov
(1971); Weins and Chan (1973); Brailsford (1976); Deutch, Felderhof, and
Saxton (1976); Felderhof and Deutch (1976); Marqusee and Ross (1984); Rice
(1985); Mattern and Felderhof (1986); Traytak (1989); Marsh (1995); Strobe,
Reiss, and Heinig (1996); Mandyam  _et al._ (1998); Tsao, Lu, and Tseng
(2001a); Biello and Samson (2015); Stark (2018). 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" Zheng and Chiew
(1989). Using the moment scattering, Bonnecaze and Brady derived analytically
the reaction rate for cubic arrays of sinks up to the quadrupole level
Bonnecaze and Brady (1991). 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 Rice (1985); Bluett and Green (2006) (see also Sec. X).

For $N>2$ the problem was solved by the images and reflections methods
analytically Deutch, Felderhof, and Saxton (1976); Labowsky (1978). It is
noteworthy that the connection between the method of images and the method of
reflections was discovered in Ref. Traytak (2003). 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 Traytak (2003).

Concerning theoretical methods to attack the multi-particle diffusion problem
Ratke and Voorhees (Ref. Ratke and Voorhees (2002)) 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 Marqusee and Ross (1984), irreducible Cartesian
tensor Beenakker and Ross (1986) or using boundary integrals and multipole
expansions Akaiwa and Voorhees (1994).". On this head we note that multiple
scattering theory Marqusee and Ross (1984) based on the microscopic monopole
approximation and multipole expansions method was performed in spherical
coordinates Akaiwa and Voorhees (1994). In its turn the irreducible Cartesian
tensor approach Beenakker and Ross (1986) is nothing more but a scalar version
of well-known induced force method (see Sec. XI).

Later on influence of many neighboring sinks on diffusion-controlled reactions
was theoretically investigated in Ref. Eun, Kekenes-Huskey, and McCammon
(2013). 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 Eun, Kekenes-Huskey, and McCammon (2013).
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 Vazquez (2010); Lee
and Lee (2020) 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…" Cedeno  _et al._ (2022). 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…" Cedeno  _et al._ (2022).
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 Bobbert, Wind, and Vlieger (1987).

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.
Wang  _et al._ (2022) 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.
Calef and Deutch (1983), whereas the current status of research on this
subject one can find in the recent survey Galanti  _et al._ (2020).

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.

##  IV 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. Eun, Kekenes-Huskey, and McCammon (2013) (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 Galanti  _et al._ (2020), so first of all we
specify the geometrical part of the problem.

###  IV.1 Domains definition

Consider a collection of arbitrary but finitely many $N\geq 1$ microscopic
spherical sinks of different radii $R_{i}$, immersed within an unbounded host
fluid medium ${\mathbb{R}}^{3}$. Hereafter the sinks are labelled sequentially
by the index $i=\overline{1,N}$, where symbol $\overline{l,m}$ means that all
integers (including $\\{\infty\\}$) from $l$ to $m$ are taken their values
successively, i.e. $\overline{l,m}:=l,l+1,\ldots,m$, where
$l(<m\leq\infty),m\in\mathbb{Z}$.

Geometrically we treat the $i$-th sink as a 3D body occupying a spherical
domain (open ball) ${\Omega}_{i}\subset{\mathbb{R}}^{3}$. Denote by
$\overline{\Omega}_{i}$ the corresponding closed bounded spherical domain,
supposing that all $\overline{\Omega}_{i}$ are non-intersecting or non-
touching, i.e. $\overline{\Omega}_{i}\cap\overline{\Omega}_{j}=\varnothing$
for $i(\neq j)=\overline{1,N}$. So the total domain occupying by array of $N$
sinks is $\Omega^{+}:=\bigcup_{i=1}^{N}{\Omega}_{i}$. Clearly this domain has
multiply connected boundary
$\partial\Omega^{+}=\bigcup\limits_{i=1}^{N}{\partial\Omega}_{i}$, where
corresponding boundary $\partial\Omega_{i}$ ($i$-th sink reaction surface Rice
(1985); Ovchinnikov, Timashev, and Belyy (1989)) is the $i$-th connected
component of the total boundary $\partial\Omega^{+}$.

On the other hand, it is apparent that the diffusion problem should be posed
on the whole 3D physical space $\mathbb{R}^{3}$ excepting the total domain
occupying by sinks, i.e., in sink-free region Torquato (2002). Hence the
required domain outside all sinks is the exterior domain with respect to
sinks:
$\Omega^{-}:=\mathbb{R}^{3}\backslash\overline{\Omega}^{+}\subset\mathbb{R}^{3}$
and the whole physical space $\mathbb{R}^{3}$ can be naturally partitioned
into two complementary subdomains $\Omega^{+}$ and $\Omega^{-}$.

Plainly, $\Omega^{-}$ is the set of accessible configurations for $B$’s so,
according to above definition II.1, we term this domain as a configuration
space. It is evident that diffusive interaction is highly affected by geometry
of the configuration space $\Omega^{-}$.

######  Example IV.1.

For instance consider two sinks labeled $i$ and $j$ with radii $R_{i}$ and
$R_{j}$, as depicted in Fig. 1. Here we placed a particle $B$ is into some
current point $P$ of the configuration space
$\Omega^{-}:=\mathbb{R}^{3}\backslash\left(\overline{\Omega}_{i}\cup\overline{\Omega}_{j}\right)$.

Figure 1: Definition sketch for two-sink array. The configuration space
$\Omega^{-}\subset\mathbb{R}^{3}$ corresponds to the given microstructure
$X^{(2)}$ with a $B$ particle located at point $P\in\Omega^{-}$. Cartesian
coordinate systems: global $\\{O;{\mathbf{r}}\\}$ (green); local in $i$-th
sink $\\{O_{i};{{\mathbf{r}}_{i}}\\}$ and in $j$-th sink
$\\{O_{j};{\mathbf{r}}_{j}\\}$ (blue); sink radii $R_{i}$ and $R_{j}$ (red)

Let $\mathbf{r}:=\overrightarrow{OP}$ be a vector that represents the position
of a current point $P$ in $\Omega^{-}$ with respect to some fixed point $O$.
Introduce a global Cartesian coordinate system, taking point $O$ as the
origin: $\\{O;\mathbf{r}\\}$ and attach the standard orthonormal Cartesian
basis $\\{{\mathbf{e}}_{\alpha}\\}_{\alpha=1}^{3}$ to this origin. The
relaited components of the position vector ${\mathbf{r}}$ we denote by
$r_{\alpha}$ ($\alpha=1,2,3$) or in ordered triple form by
$\left(r_{1},r_{2},r_{3}\right)$.

Using the global coordinates, we designate by $\mathbf{r}^{i}_{0}$ the
position of the center of $i$-th sink and by
$\mathbf{L}_{ij}:=\overrightarrow{O_{j}O_{i}}=\mathbf{r}^{i}_{0}-\mathbf{r}^{j}_{0}$
the separation vector between the two sinks centers, pointing from $j$-th sink
toward $i$-th sink (see Fig. 1).

Clearly, any $i$-th sink spherical domain is:
$\Omega_{i}=\\{\mathbf{r}\in{\mathbb{R}}^{3}:\|\mathbf{r}-\mathbf{r}^{i}_{0}\|<R_{i}\\}$
with the boundary
$\partial\Omega_{i}=\\{\mathbf{r}\in{\mathbb{R}}^{3}:\|\mathbf{r}-\mathbf{r}^{i}_{0}\|=R_{i}\\}$.

We shall consider diffusion of $B$’s in the configuration manifold
$\Omega^{-}$ formed by an arbitrary but finite ensemble of $N$ spherical sinks
with fixed radii $\left\\{R_{i}\right\\}_{i=1}^{N}$ centered at positions
$\left\\{{\mathbf{r}}^{i}_{0}\right\\}_{i=1}^{N}$ immersed in the 3D host
medium. Hence the geometry of this manifold is completely determined by
$N$-sink configuration Torquato (2002) or briefly microstructure Brown and
Peterseim (2016):

$\displaystyle
X^{(N)}:=\left\\{\left({\mathbf{r}}^{i}_{0},R_{i}\right)\right\\}_{i=1}^{N}\,.$
(11)

In Ref. Traytak (2003) 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
V.1) we consider here a general formalism based on the manifold conception
Traytak (2003). We shall show here that the manifold conception arises simply
for the above microstructure (11). 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 Sawyer (1971), 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."

###  IV.2 Manifold structure on $\Omega^{-}$

Let us introduce the manifold structure on the configuration space
$\Omega^{-}$. First let us introduce open complements of the sink domains
$\Omega_{i}$ in the physical space ${\mathbb{R}}^{3}$, which play an important
role in our study:
$\Omega^{-}_{i}={\mathbb{R}}^{3}\backslash\overline{\Omega}_{i}$. It is clear
that for all arrays ($N\geq 1$) we have:
$\Omega^{-}=\bigcap\limits_{i=1}^{N}\Omega_{i}^{-}\subseteq\Omega_{i}^{-}$.

Configuration space $\Omega^{-}$ is a 3D smooth manifold (strictly speaking
sub-manifold in $\mathbb{R}^{3}$) Tu (2011) since: (a) there is its coverage
$\bigcup\limits_{i=1}^{N}\Omega_{i}^{-}\supset\Omega^{-}$; (b) $\Omega^{-}$ is
a metric space as a subset of $({\mathbb{R}}^{3},\|\cdot\|)$; (c) there exists
smooth homomorphisms $\varphi_{i}:\Omega_{i}^{-}\to{\mathbb{R}}^{3}$; (d) For
intersection
$\Omega^{-}_{ij}:=\Omega^{-}_{i}\cap\Omega^{-}_{j}\neq\varnothing$ the overlap
map:
$\varphi_{ij}:=\varphi_{j}\circ{\varphi_{i}}^{-1}:\varphi_{i}\left(\Omega^{-}_{ij}\right)\to\varphi_{j}\left(\Omega^{-}_{ij}\right)$
is smooth for each pair of indices $i(\neq j)=\overline{1,N}$.

A coordinate chart (or just a chart) on $\Omega^{-}$ is a pair
$\left(\Omega^{-}_{i},\varphi_{i}\right)$. A collection of all charts
$\left(\Omega^{-}_{i},\varphi_{i}\right)$ ($i=\overline{1,N}$) is called a
smooth atlas on $\Omega^{-}$, which defines the manifold structure on
$\Omega^{-}$. For short in the following we shall simply term $\Omega^{-}_{i}$
the $i$-th chart. It is clear that boundary $\partial\Omega^{-}$ is a 2D
manifold with empty boundary.

Consider the position vector ${\mathbf{r}}$ of point $P\in\Omega^{-}_{i}$,
then
$\varphi_{i}\left(\mathbf{r}\right)=(r_{1}^{i},r_{2}^{i},r_{3}^{i})\in\mathbb{R}^{3}$.
These $3$ real numbers are called local coordinates of point $P$ in the chart
$\Omega^{-}_{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}\\}\equiv\\{O_{i};{\mathbf{r}}_{i}\\}$
with the same orientation as the global one and the origin
$O_{i}({\mathbf{r}}^{i}_{0})$ at the center of the $i$-th sink (see Fig. 1).
So, for any fixed $i$ we can choose standard local orthonormal Cartesian basis
$\\{{\mathbf{e}}_{\alpha}^{i}\\}_{\alpha=1}^{3}$ attached to the origin
$O_{i}({\mathbf{r}}^{i}_{0})$, such that $r^{i}_{\alpha}$ are coordinates in
this basis. Introdused local coordinates are useful for computation of the
diffusion field in arrays wth microstructure (11).

Let us take above position vector ${\mathbf{r}}$ for the point $P$ and
consider it with respect to two local Cartesian coordinate systems
$\\{O_{i};{{\mathbf{r}}_{i}}\\}$ and $\\{O_{j};{\mathbf{r}}_{j}\\}$.
Corresponding point $P$ position vectors are

$\displaystyle\varphi_{i}\left(\mathbf{r}\right)={\mathbf{r}}_{i}=\mathbf{r}-\mathbf{r}^{i}_{0}\quad\mbox{and}\quad\varphi_{j}\left(\mathbf{r}\right)=\mathbf{r}_{j}=\mathbf{r}-\mathbf{r}^{j}_{0}\,.$
(12)

It is also evident that
$\varphi_{i}^{-1}\left(\mathbf{r}_{i}\right)=\mathbf{r}_{i}+\mathbf{r}_{0}^{i}=\mathbf{r}$.
Perform the transition between these local coordinates given by the overlap
map
$\varphi_{ij}=\varphi_{j}\circ{\varphi_{i}}^{-1}:{\mathbf{r}}_{i}\rightarrow{\mathbf{r}}_{j}$
with the rule (see Fig. 1):

$\displaystyle{\mathbf{r}}_{j}=\varphi_{ij}\left({\mathbf{r}}_{i}\right)=\varphi_{j}\left(\left({{\varphi_{i}}^{-1}\left({\mathbf{r}}_{i}\right)}\right)\right)=\mathbf{r}_{i}+\mathbf{L}_{ij}\,,\quad\mbox{where}\quad\mathbf{L}_{ij}=\mathbf{r}^{i}_{0}-\mathbf{r}^{j}_{0}\,.$
(13)

Hereafter $\mathbf{L}_{ij}$ being the vector, connecting centers of sinks $i$
and $j$. In this work we also use the corresponding unit vector
${\widehat{\mathbf{L}}}_{ij}:={\mathbf{L}}_{ij}/L_{ij}$, where $L_{ij}$ is the
distance between above centers
$L_{ij}=\|{\mathbf{L}}_{ij}\|=\|\mathbf{r}^{i}_{0}-\mathbf{r}^{j}_{0}\|=\left\|{\mathbf{r}}_{i}-{\mathbf{r}}_{j}\right\|$.
Transformations (12) and (13) are clearly illustrated by the following
triangle commutative diagram:

${\\{O_{i};{\mathbf{r}}_{i}\\}}$${\\{O_{j};{\mathbf{r}}_{j}\\}}$${\\{O;{\mathbf{r}}\\}}$${\varphi_{i}}^{-1}$$\varphi_{ij}$$\quad\varphi_{j}$

######  Remark IV.1.

Note that simple linear connection between local coordinates (13) 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 $\Omega^{-}$ can be called
the configuration manifold.

An additional point to emphasize is that theoretical calculations are often
used the distance between surfaces of domains $\Omega_{i}$ and $\Omega_{j}$,
i.e. $h_{ij}:=L_{ij}-\left(R_{i}+R_{j}\right)$ (see Fig. 1). Hence the
condition describing configurations of $N$ sinks, which are not intersecting
or touching may be written explicitly as

$\displaystyle h_{ij}>0\quad\mbox{for all}\quad i(\neq j)=\overline{1,N}\,.$
(14)

###  IV.3 Governing equation

To describe the steady state microscopic diffusion-reaction model introduce a
local concentration (number density) of $B$-reactants,
$n:\Omega^{-}\rightarrow\mathbb{R}_{+}$ depends on a given microstructure
$X^{(N)}$ (11) and also different reactivities of sinks. To describe the
latter it is convenient to introduce the set of sinks reactivities
$k_{in}^{(N)}:=\left\\{{k}^{i}_{in}\right\\}_{i=1}^{N}$, where rate
${k}^{i}_{in}$ corresponds to the $i$-th sink reactivity. Thus, in a fixed
global Cartesian coordinate system $\\{O;\mathbf{r}\\}$, extending notations
of Refs. Brailsford (1976); Torquato (2002), one can write this function as
$n\left(P;X^{(N)},\,k_{in}^{(N)}\right)$. However, for short we denote it just
by $n_{N}\left(P\right)$. Moreover, to investigate the effects of diffusive
interaction as in the case of one sink Eq. (4) we introduce normalized
concentration of $B$ particles related to their rescaled local concentration

$\displaystyle u_{N}:\Omega^{-}\to\left(0,1\right)\,,\quad\mbox{where}\quad
u_{N}\left(P\right)=1-\frac{n_{N}\left(P\right)}{c_{B}}\,.$ (15)

Assuming that rescaled Fick’s local diffusion flux constitutive relation holds
Rice (1985)

$\displaystyle\mathbf{j}=D{\bm{\nabla}}u_{N}\left(P\right)\quad\mbox{in}\quad\Omega^{-}$
(16)

the rescaled local concentration $u_{N}\left(P\right)$ is governed by the
$B$-particles conservation law which leads to the Laplace equation

${\bm{\nabla}}^{2}u_{N}\left(P\right)=0\quad\mbox{in}\quad\Omega^{-}\,.$ (17)

Hereinafter as usual
${\bm{\nabla}}:=\sum_{\alpha=1}^{3}{\mathbf{e}}_{\alpha}{\partial_{\alpha}}$
and ${\bm{\nabla}}^{2}$ stand for the gradient operator and the Laplacian on
$\mathbb{R}^{3}$, respectively.

To describe diffusion-reaction phenomena, Eq. (17) should be subjected to the
appropriate boundary conditions on the sink surfaces $\partial\Omega_{i}$ and
condition at infinity.

###  IV.4 General boundary conditions

The present work is concerned with the whole possible range of the intrinsic
rate constants $0\leq{k}^{i}_{in}<+\infty$ and in general case we can impose
the local form of inhomogeneous Robin boundary conditions

$\left.\left[4\pi
R_{i}^{2}D\,\hat{\mathbf{n}}_{i}\cdot{\bm{\nabla}}u_{N}-{k}^{i}_{in}\left(f_{i}\left(P_{i}\right)-u_{N}\right)\right]\right|_{\partial\Omega_{i}}=g_{i}\left(P_{i}\right)\,,\quad
i=\overline{1,N}\,;$ (18)

where $\hat{\mathbf{n}}_{i}$ being the outer-pointing unit normal vector with
respect to $\partial\Omega_{i}$ (see Fig. 2) and point $P_{i}\in\Omega_{i}$.
We assume for simplicity that $f_{i}\,,g_{i}\in C(\partial\Omega_{i})$.

Figure 2: A geometrical sketch for three kinds of sinks in the configuration
manifold $\Omega^{-}$: partialy reflecting $\Omega_{1}$ (gray), fully
absorbing $\Omega_{2}$ (black) and fully reflecting $\Omega_{3}$ (white);
$\hat{\mathbf{n}}_{i}$ are the corresponding outward to $\Omega^{-}$ unit
normal vectors; and small green circles depict $B$ particles

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

(i) In the case when $k_{in}^{i}\to 0$ conditions (18) give inhomogeneous
Neumann boundary conditions

$\displaystyle 4\pi
R_{i}^{2}D\left.\hat{\mathbf{n}}_{i}\cdot{{\bm{\nabla}}u_{N}}\right|_{\partial\Omega_{i}}=g_{i}\left(P_{i}\right)\,.$
(19)

(ii) Provided $k_{in}^{i}\rightarrow\infty$, conditions (18) are reduced to
inhomogeneous Dirichlet boundary conditions

$\displaystyle\left.u_{N}\right|_{\partial\Omega_{i}}=f_{i}\left(P_{i}\right).$
(20)

For uniqueness of a solution to Laplace’s equation (17) under conditions (18)
one should require satisfying the regularity conditions at infinity:
$u_{N}\left(P\right)={\cal{O}}(\|{\mathbf{r}}\|^{-1})$ and
${{\bm{\nabla}}u_{N}}\left(P\right)={\cal{O}}(\|{\mathbf{r}}\|^{-2})$ as
$\|{\mathbf{r}}\|\rightarrow\infty$ Tikhonov and Samarskii (1963); Krutitskii
(1999). These conditions are equivalent to the following regularity condition
at infinity

$\displaystyle\left.u_{N}\left(P\right)\right|_{\|{\mathbf{r}}\|\rightarrow\infty}\rightrightarrows
0\,$ (21)

and it uniquely determines the function $u_{N}\left(P\right)$, which is
harmonic in the infinite domain $\Omega^{-}$ Tikhonov and Samarskii (1963).
Hereafter the sign $\rightrightarrows$ denotes the uniform limit.

Thus, the exterior inhomogeneous Robin boundary value problem (17), (18) and
(21) is formulated. The existence and uniqueness of the classical solution to
the posed problem $u_{N}\left(P\right)$ are a matter of common knowledge
Tikhonov and Samarskii (1963). Furthermore, it is well-known that the
smoothness class of solutions $u_{N}\left(P\right)$ depends on a boundary
operator: (a) for the Robin and Neumann problems $u_{N}\left({P}\right)\in
C^{2}\left(\Omega^{-}\right)\cap C^{1}\left(\overline{\Omega}^{-}\right)$; (b)
for the Dirichlet problem $u_{N}\left({P}\right)\in
C^{2}\left(\Omega^{-}\right)\cap C\left(\overline{\Omega}^{-}\right)$
Ladyzhenskaya and Uralt’seva (1968); Krutitskii (1999).

###  IV.5 Boundary conditions for diffusion-influenced reactions

It must be emphasized that the general boundary value problem (17), (18) and
(21) describes a vast amount of diffusion-reaction phenomena, occurring in
nature and industry (see Sec. I). 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}\left(P_{i}\right)\equiv 1$ and
$g_{i}\left(P_{i}\right)\equiv 0$, describing the microscopic Collins-Kimball
diffusion-influenced reactions Rice (1985).

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 Rice (1985). Provided magnitudes of ${k}^{i}_{in}$ are
positive and finite, we deal with diffusion-influenced reactions on so-called
partially reflecting (absorbing) sinks Bedeaux, Kapral, and Mazur (1977);
Lebenhaft and Kapral (1979). In the case of very large values of
${k}^{i}_{in}$ (i.e., as ${k}^{i}_{in}\to+\infty$), the sinks are called fully
(perfectly) absorbing, and the corresponding reactions, occurring on them, are
known as fully diffusion-controlled Rice (1985). Physicall that means that the
activation barrier appeared to be very small Kumpulainen  _et al._ (2017).
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 $\left\\{1<i\leq N\right\\}$ Doktorov (2020). At first
glance it would seem that this case is described by boundary conditions (18)
automatically. However, solution of the boundary value problem (17), (18) and
(21) may be faced subtle mathematical difficulties Krutitskii (1999); Traytak
and Price (2007).

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.

######  Definition IV.1.

If the boundary conditions differ on different parts of a given connected
component of a boundary $\partial\Omega_{i}$ these conditions are referred to
as proper mixed boundary conditions.

Physically the proper mixed boundary condition means that a given connected
component $\partial\Omega_{i}$ is characterized by a prescribed heterogeneous
reactivity. So from analytical viewpoint one should define a piecewise
constant function ${k}^{i}_{in}:\partial\Omega_{i}\rightarrow[0+\infty)$.

######  Definition IV.2.

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:
$\left.{k}^{i}_{in}\right|_{\partial\Omega_{i}}=const$. Figure 2 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. Traytak
and Price (2007); Holcman and Schuss (2015).

In concluding this section, we note, that crowding and geometrical constraints
effects on diffusion-influenced reactions Lee and Lee (2020) may be treated if
one considers microstructures $X^{(N)}$ containing $1\leq N_{0}\leq N$ inert
obstacles. The presence of fully reflecting obstacles might strongly influence
the reaction rate when the number of these obstacles is large enough Eun,
Kekenes-Huskey, and McCammon (2013); Lee and Lee (2020) (see also numerical
calculations of the reaction rates for the mixed Dirichlet-Robin problem at
$N=3$ in Sec. X).

###  IV.6 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 $\left(O_{i};\mathbf{r}_{i}\right)$ at the $i$-th sink
of radius $R_{i}$ as a characteristic unit:

$\displaystyle{\bm{\xi}}_{i}=\mathbf{r}_{i}/R_{i}\,,\quad{\bm{\xi}}^{i}_{0}=\mathbf{r}^{i}_{0}/R_{i}\quad\xi^{i}_{\nu}=r^{i}_{\nu}/R_{i}\,,\quad\xi_{i}=\|{\bm{\xi}}_{i}\|\,.$
(22)

Besides these dimensionless local coordinates for any pair of sinks $i$ and
$j(\neq i)$ the following important parameters Traytak (2006)

$\displaystyle\varepsilon_{ji}=R_{i}/L_{ij}<1\,,\quad\varepsilon_{ij}=R_{j}/L_{ij}<1$
(23)

are naturally arisen. Clearly, parameters (23) are totally determined by a
given microstructure $X^{(N)}$ (11) and in general case they obey the
condition $\varepsilon_{ij}\neq\varepsilon_{ji}$. Moreover, the unit sphere
centered at a point ${\bm{\xi}}^{i}_{0}$ we henceforth designate by

$\displaystyle\partial\Omega_{i}^{1}\equiv\partial\Omega_{i}\left({\bm{\xi}}^{i}_{0};1\right):=\\{{\bm{\xi}}_{i}\in\Omega_{i}^{-}:\xi_{i}=1\\}\,.$
(24)

Here it should be especially noted that performed normalization (22) holds
only locally in charts $\Omega_{i}^{-}$ and does not hold in the whole
manifold $\Omega^{-}$. Note also that in the interest of readability for all
charts $\Omega_{i}^{-}$, transformed according to Eq. (22), we retain the same
designations.

It is expedient now to write the diffusion problem immediately in the local
coordinates $\\{O_{i};{\bm{\xi}}_{i}\\}$. However, sometimes to emphasize the
fact that we deal with the field on a manifold coordinateless form
$u_{N}\left(P\right)$ will be used. Clearly, the basic equation governing
exterior steady-state Eq. (17) under partially reflecting boundary conditions
(18) and regularity condition at infinity (21) with respect to the local
concentration $u_{N}\left({\bm{\xi}}_{i}\right)$ ($i=\overline{1,N}$) and
dimensionless variables (22) reads

$\displaystyle{\mathbf{\nabla}}_{{\bm{\xi}}_{i}}^{2}\,u_{N}\left({\bm{\xi}}_{i}\right)=0\quad\mbox{in}\quad\Omega_{i}^{-}\,,$
(25)
$\displaystyle-\left.\partial_{\xi_{i}}u_{N}\left({\bm{\xi}}_{i}\right)\right|_{\partial\Omega_{i}^{1}}=\kappa_{i}\left[1-\left.u_{N}\left({\bm{\xi}}_{i}\right)\right|_{\partial\Omega_{i}^{1}}\right]>0\,,$
(26)
$\displaystyle\left.u_{N}\left({\bm{\xi}}_{i}\right)\right|_{\xi_{i}\rightarrow\infty}\rightrightarrows
0\,.$ (27)

Hereafter ${\bm{\nabla}}_{{\bm{\xi}}_{i}}^{2}$ stands for the corresponding
dimensionless Laplacian written in local coordinates
$\\{O_{i};{\bm{\xi}}_{i}\\}$.

Besides we introduce the corresponding dimensionless set of sinks
reactivities: $\kappa^{(N)}:=\left\\{\kappa_{i}\right\\}_{i=1}^{N}$ with
$\kappa_{i}=k_{in}^{i}/k_{S}^{i}\in\left[0,\infty\right)$, where
$k_{S}^{i}=4\pi R_{i}D$ is the Smoluchowski rate constant (9) for the fully
absorbing $i$-th sink.

Thus, from a mathematical viewpoint, we deal with the exterior Robin boundary
value problem (25)-(27) in an unbounded 3D manifold $\Omega^{-}$ with
$N$-connected boundary $\partial\Omega^{-}$. An additional point to emphasize
is that the boundary conditions (26), being posed on all connected components
$\partial\Omega_{i}$ of the boundary $\partial\Omega^{-}$, are integral, which
reflects the important fact that diffusive interaction between sinks are not
pairwise additive.

###  IV.7 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}:={\Phi_{i}}/{c_{B}}$, where $\Phi_{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}\left(P\right)$, the total microscopic trapping rate on the reaction
surface $\partial\Omega_{i}$ can be calculated straightforwardly

$k_{i}=\oint_{\partial\Omega_{i}}\left.\hat{\mathbf{n}}_{i}\cdot\mathbf{j}\right|_{\partial\Omega_{i}}dS_{i}\,,$
(28)

where $dS_{i}$ is the $i$-th surface element. Then, recasting (28) with
respect to $u_{N}\left({\bm{\xi}}_{i}\right)$ in an exterior neighborhood of
the $i$-th sink surface, one can easily obtain the desired microscopic
trapping rate (28) by the surface integral over the unit sphere

$\displaystyle
k_{i}(X^{(N)};\kappa^{(N)})=-\oint_{\partial\Omega_{i}^{1}}\left.\partial_{\xi_{i}}u_{N}\right|_{\partial\Omega_{i}^{1}}d{\hat{\mathbf{\xi}}}_{i}\,.$
(29)

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 (29) is sought
in the form

$\displaystyle
k_{i}(X^{(N)};\kappa^{(N)})=k_{S}^{i}J_{0}^{i}\left(\kappa_{i}\right)J_{i}(X^{(N)};\kappa^{(N)})\,,$
(30)

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

$\displaystyle
J_{0}^{i}\left(\kappa_{i}\right)=\kappa_{i}/\left(1+\kappa_{i}\right)\quad\mbox{for}\quad
i=\overline{1,N}\,.$ (31)

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 (30) the magnitudes $J_{i}$ are commonly called the
screening coefficients Cossali and Tonini (2020) (rate correction factors
Labowsky (1978)). It is absolutely clear that due to dependence on
configuration of sinks $X^{(N)}$ and their reactivities $\kappa^{(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

$\displaystyle 0<J_{i}<1\,,\quad\lim_{\varepsilon_{ij}\to
0}J_{i}=1\quad\mbox{for all}\quad i(\neq j)=\overline{1,N}\,.$ (32)

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

It is significant that, generally speaking, the microscopic trapping rate (29)
for $N\geq 2$ does not coincide with the reaction rate constant $k$ defined by
Eq. (2), which is a macroscopic value Bressloff (2020). So, investigating
diffusive interaction effects in systems with $N\geq 2$ sinks we cannot use
the term "rate constant" for $k_{i}$ given by Eq. (29) any more Traytak
(1992).

##  V 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 (25)-(27).

However, before starting to discuss this topic, we briefly highlight the
application of the boundary integral equations method. Mikhlin (1964);
Grebenkov and Traytak (2019) One can reduce the above diffusion problem
(25)-(27) 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.
Traytak (2005) Despite the fact that this method is also applicable to the
mixed boundary conditions Krutitskii (1999) its implementation in the
diffusion problem at issue involves a number of disadvantages compared to the
approach that would be considered here.

###  V.1 On the method of separation of variables

First recall some well-known mathematical facts on the standard method of
separation of variables (MSV) Tikhonov and Samarskii (1963). 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;\mathbf{r}\\}$ in a
bounded or unbounded domain $\Omega\subset\mathbb{R}^{3}$.

######  Definition V.1.

A domain $\Omega$ is said to be canonical if any solution
$u\left(\mathbf{r}\right)$ to a boundary value problem in $\Omega$ may be
derived by the standard MSV.

Thus, in other words, there exists an orthogonal functional basis of solutions
$\\{\varphi_{n}\left(\mathbf{r}\right)\\}_{n=0}^{\infty}$ (called basis
solutions) to a boundary value problem in a canonical domain
$\Omega\subset\mathbb{R}^{3}$ and a set of real coefficients
$\\{C_{n}\\}_{n=0}^{\infty}$ such that

$\displaystyle
u\left(\mathbf{r}\right)=\sum_{n=0}^{\infty}C_{n}\varphi_{n}\left(\mathbf{r}\right)\quad\mbox{in}\quad\Omega\,,$
(33)

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 $\mathbb{R}^{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
$\Omega_{i}$ and their complements $\Omega_{i}^{-}$ ($i=\overline{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.

###  V.2 Generalized method of separation of variables

The generalized method of separation of variables (GMSV) Traytak (2003);
Traytak and Grebenkov (2018); Grebenkov and Traytak (2019); Galanti  _et al._
(2020) naturally stems from the standard method of separation of variables and
its incipience go back to the classic works by Maxwell Maxwell (1873) and
Basset. Basset (1887) 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
Rayleigh (1892). 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" Nikolaev and Tanchik (2016);
Miroshnikov  _et al._ (2022). Guz’ and Golovchan treated the method as
essentially a particular case of the "method of series" Guz’ and Golovchan
(1972). Finally, in micromechanics of heterogeneous materials mainly term
"multipole expansion approach" is used Kushch (2013).

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" Ivanov (1970). 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 Ivanov (1970). 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
Ivanov (1970). 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.

Figure 3: A block diagram of the generalized method of separation of variables
for diffusion problems in the manifold $\Omega^{-}$

The basic idea of the GMSV consists in reducing the original boundary value
problem posed in the configuration manifold $\Omega^{-}$ to the corresponding
$N$ boundary value problems in charts $\Omega^{-}_{i}$. These $N$ problems are
coupled through the boundary conditions, prescribed on the whole disconnected
boundary $\partial\Omega^{-}$. Thus, the GMSV is actually the method of
separation of variables on the 3D manifold Traytak (2003).

###  V.3 Outline of approach

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

  1. 1.

Proof of the well-posedness of the solution to the posed boundary value
problem;

  2. 2.

Deconposition of the desired solution in the form of the sum of partial
solutions by means of the general linear superposition principle;

  3. 3.

Determination of the appropriate basis solutions to the diffusion boundary
value problem posed in the given canonical domains with disconnected boundary;

  4. 4.

Determination of relevant addition theorems to obtain local regular expansions
of partial solutions in a neighborhood of corresponding connected components
of bounday;

  5. 5.

Application of the addition theorems in order to satisfy the boundary
conditions;

  6. 6.

Reduction of the problem to a resolving self-consistent infinite system of
linear algebraic equations;

  7. 7.

Solution of the resolving infinite system of linear algebraic equations;

  8. 8.

Calculation of the microscopic reaction rate.

In Fig. 3 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. 3) 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.

##  VI 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 Grebenkov and Traytak (2019). 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. Galanti  _et al._
(2016); Grebenkov and Traytak (2019); Galanti  _et al._ (2020). 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 Appendix. Background on the irreducible Cartesian tensors). The wide
use of solid harmonics in the theory of diffusion-influenced reactions is a
matter of common knowledge (see, e.g., Refs. Venema and Bedeaux (1989); Tsao,
Lu, and Tseng (2001b); Tsao (2002); McDonald and Strieder (2003, 2004) and
references cited there). In the meantime, studies devoted to similar
applications of ICT are only very few Lebenhaft and Kapral (1979); Beenakker
and Ross (1986); Traytak (1991, 1992, 2001, 2003). 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. XI). 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 (11) Erdmann, Kröger, and Hess (2003); Baranov and Martynov
(2010); Agre (2011); Gus’kov (2013); Singh, Ghose, and Adhikari (2015); Hess
(2015); Mane (2016); Bechinger  _et al._ (2016); Manna, Kumar, and Adhikari
(2017); Andrews (2020); Turk, Singh, and Adhikari (2022).

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 (25)-(27) Lebenhaft and Kapral
(1979). Subsequently, an explicit form of the GMSV with the aid of ICT was
carried out in our paper Ref. Traytak (1991), 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 Schmitz and
Felderhof (1978). Although Ref. Schmitz and Felderhof (1978) 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 $\xi_{i}$ and corresponding Cartesian local coordinates
$\xi_{\gamma_{m}}^{i}$ in every $i$-th chart ($i=\overline{1,N}$) Schmitz and
Felderhof (1978)

$\displaystyle
u_{N}\left({\bm{\xi}}_{i}\right)=\sum\limits_{n=0}^{\infty}w_{n}\left({\xi}_{i}\right)C_{\gamma_{1}\ldots\gamma_{n}}^{i}\odot^{n}\mathop{{\xi}_{\gamma_{1}}^{i}{\ldots\xi}_{\gamma_{n}}^{i}}\limits^{\rule[-2.8903pt]{0.72229pt}{3.61371pt}\rule{33.96698pt}{0.72229pt}\rule[-2.8903pt]{0.72229pt}{3.61371pt}}\quad\mbox{in}\quad\Omega_{i}^{-}\,.$
(34)

Hereinafter
$\mathop{{\xi}_{\gamma_{1}}^{i}{\ldots\xi}_{\gamma_{n}}^{i}}\limits^{\rule[-2.8903pt]{0.72229pt}{3.61371pt}\rule{33.96698pt}{0.72229pt}\rule[-2.8903pt]{0.72229pt}{3.61371pt}}$
is the ICT defined by Eq. (A.1). The explicit form of functions
$w_{n}\left({\xi}_{i}\right)$ may be found by substituting in Laplace’s
equation (25) while the unknown tensor coefficients
$C_{\gamma_{1}\ldots\gamma_{n}}^{i}$ to be determined from the boundary
conditions (26).

Note in passing that so-called self-consistent field method that is another
version of the above method of Cartesian ansatz Gus’kov (2013). A method of
reflections with the help of ICT used in Ref. Keh and Chen (1995) 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 Kim (1987); Tsao, Lu, and Tseng
(2001a).

In Ref. Beenakker and Ross (1986) 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 Traytak (2003, 2005).

##  VII The addition theorem in terms of the irreducible tensors

Due to the fact that the Laplace equation (25) remains invariant under the
action of the group of Euclidean motions Miller (1977) irregular solid
harmonics (see definitions in Appendix Appendix. Background on the irreducible
Cartesian tensors) written in local coordinates of any chart $\Omega_{i}^{-}$
($\Omega_{j}^{-}$) can be recast in terms of regular solid harmonics written
in local coordinates of another chart $\Omega_{j}^{-}$ ($\Omega_{i}^{-}$).

One of the key point of the GMSV algorithm is the transformation of the solid
harmonics under action of subgroup of translations (1) Steinborn and
Ruedenberg (1973); Galanti  _et al._ (2020). Provided charts $\Omega_{i}^{-}$
and $\Omega_{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$\to$R TAT).

######  Remark VII.1.

On the contrary, in the case when forms of domains $\Omega_{i}^{-}$ and
$\Omega_{j}^{-}$ are different (like the exteriors of a sphere
$\Omega_{1}^{-}$ and a prolate spheroid $\Omega_{2}^{-}$ see Ref. Traytak and
Grebenkov (2018)) 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$\to$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.

###  VII.1 The degenerate translation addition theorem

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

Assume that $r>\|\mathbf{r}^{i}_{0}\|$, where $r=\|{\mathbf{r}}\|$ is the
distance between the global origin $O$ and a point $P$ (see Fig. 1).

Thereby, we can write the known Taylor’s expansion

$\displaystyle\frac{1}{\|{\mathbf{r}}_{i}\|}=\frac{1}{\|{\mathbf{r}}-\mathbf{r}^{i}_{0}\|}=\sum_{k=0}^{\infty}\frac{(-1)^{k}}{k!}\left(\mathbf{r}^{i}_{0}\cdot{\mathbf{\nabla}}\right)^{k}\left(\frac{1}{r}\right)\,.$
(35)

Clearly, series (35) converges absolutely and uniformly for all $r>0$ Tikhonov
and Samarskii (1963).

Function $1/r$ is often termed the generating function. Traditionally, Eq.
(35) is referred to as the multipoles expansion of the (global) fundamental
solution of the Laplacian at the pole point $\mathbf{r}^{i}_{0}$ in powers of
$r^{-1}$ (see, e.g., Ref. Tikhonov and Samarskii (1963)).

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

######  Theorem VII.1.

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

$\displaystyle\frac{1}{\|\mathbf{r}-\mathbf{r}^{i}_{0}\|}=\sum\limits_{k=0}^{\infty}\omega_{k}r^{-(2k+1)}\mathop{\hskip
5.69046pt{r}_{\gamma_{1}}{\ldots
r}_{\gamma_{k}}}\limits^{\rule[-2.8903pt]{0.72229pt}{3.61371pt}\rule{32.52127pt}{0.72229pt}\rule[-2.8903pt]{0.72229pt}{3.61371pt}}\odot^{k}\mathop{r_{\gamma_{1}}^{0i}{\ldots
r}_{\gamma_{k}}^{0i}}\limits^{\rule[-2.8903pt]{0.72229pt}{3.61371pt}\rule{32.52127pt}{0.72229pt}\rule[-2.8903pt]{0.72229pt}{3.61371pt}}\,,$
(36)

where $\omega_{k}:=\left(2k-1\right)!!/k!$, $r_{\gamma_{\nu}}$ and
$r_{\gamma_{m}}^{0i}$ are Cartesian coordinates of vectors $\mathbf{r}$ and
$\mathbf{r}^{i}_{0}$, respectively.

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

###  VII.2 The translation addition theorem

It has been known that I$\to$R TAT holds true for solid harmonics
$\psi^{\pm}_{nm}\left({\mathbf{r}}\right)$ defined in polar spherical
coordinate systems (see, e.g., Ref. Grebenkov and Traytak (2019) and
references therein). Here we shall present a derivation of the similar I$\to$R
TAT in terms of ICT $\bm{X}^{\pm}_{n}$ (see Appendix Appendix. Background on
the irreducible Cartesian tensors), which was given in Refs. Traytak (1991,
1992) without a proof.

To calculate the local concentration $u_{N}\left({P}\right)$ at any given
point $P\in\Omega^{-}\subset\Omega^{-}_{i}\cap\Omega^{-}_{j}$ with position
vector ${\mathbf{r}}$ we consider functions connected with every coordinate
chart $\left(\Omega^{-}_{i},\varphi_{i}\right)$:
$u^{\prime}\left(\mathbf{r}_{i}\right)=u_{N}\left({P}\right)$ for
$\mathbf{r}_{i}=\varphi_{i}\left({\mathbf{r}}\right)$ ($i=\overline{1,N}$)
(see IV.2). In any other chart $\left(\Omega_{j}^{-},\varphi_{j}\right)$ at
the same point $P\in\Omega^{-}$ we have alternatively

$\displaystyle
u^{\prime\prime}\left(\mathbf{r}_{j}\right)=u^{\prime}\left(\mathbf{r}_{i}\right)=u_{N}\left({\mathbf{r}}\right)\,,$
(37) $\displaystyle
u^{\prime\prime}\left(\varphi_{ij}\left({\mathbf{r}}_{i}\right)\right)=u^{\prime}\left(\varphi_{ji}\left({\mathbf{r}}_{j}\right)\right)=u_{N}\left({\mathbf{r}}\right)\,.$
(38)

For our problem the overlap map $\varphi_{ij}$ is a simple translation (13) so
let us start with an important general

######  Definition VII.1.

Any two given local curvilinear coordinate systems
$\\{O_{i};{\bm{\zeta}}_{i}\\}$ and $\\{O_{j};{\bm{\zeta}}_{j}\\}$ are said to
be consistent if there exists a translation (1) with vector ${\mathbf{a}}$
such that
$\mathbf{\bm{\zeta}}_{j}=T_{\mathbf{a}}{\bm{\zeta}}_{i}={\bm{\zeta}}_{i}+{\mathbf{a}}$.

This definition immediately brings to

######  Proposition VII.1.

Any two local Cartesian coordinate systems $\\{O_{i};{\mathbf{r}}_{i}\\}$ and
$\\{O_{j};{\mathbf{r}}_{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};{{\mathbf{r}}_{i}}\\}$
and $\\{O_{j};{\mathbf{r}}_{j}\\}$ depicted in Fig. 1 are consistent. Note in
passing that some authors prefer to use term "parallel coordinate systems"
rather than "consistent coordinate systems" Steinborn and Ruedenberg (1973).

######  Definition VII.2.

A translation addition theorem for a given smooth function
$u:\Omega^{-}\to\mathbb{R}_{+}$ is a formula expressing its value
$u(\mathbf{r}+\mathbf{a})$ in terms of values $u(\mathbf{r})$ and
$u(\mathbf{a})$, where $\mathbf{r},\mathbf{r}+\mathbf{a}\in\Omega^{-}$, and
also their derivatives.

######  Definition VII.3.

The translation addition theorem I$\to$R yields an explicit formula for
converting an irregular ICT $\bm{X}_{n}^{-}$ given in one origin $O_{j}$ into
a local expansion with respect to the corresponding regular ICT
$\bm{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.

######  Theorem VII.2.

The irregular to regular translation addition theorem for the irreducible
Cartesian tensors holds true

$\displaystyle r_{j}^{-\left(2n+1\right)}\mathop{{r}_{\gamma_{1}}^{j}{\ldots
r}_{\gamma_{n}}^{j}}\limits^{\rule[-2.8903pt]{0.72229pt}{3.61371pt}\rule{36.135pt}{0.72229pt}\rule[-2.8903pt]{0.72229pt}{3.61371pt}}=\sum\limits_{k=0}^{\infty}U^{ij}_{\gamma_{1}\ldots\gamma_{n}\mu_{1}\ldots\mu_{k}}\left({\widehat{\mathbf{L}}}_{ij}\right)\odot^{k}\mathop{{r}_{\mu_{1}}^{i}{\ldots
r}_{{\mu_{k}}}^{i}}\limits^{\rule[-2.8903pt]{0.72229pt}{3.61371pt}\rule{36.135pt}{0.72229pt}\rule[-2.8903pt]{0.72229pt}{3.61371pt}}\,,$
(39)

where the transformation matrix elements read

$\displaystyle
U^{ij}_{\gamma_{1}\ldots\gamma_{n}\mu_{1}\ldots\mu_{k}}\left({\widehat{\mathbf{L}}}_{ij}\right)=\sigma_{kn}L_{ij}^{-\left(n+k+1\right)}\Lambda_{\gamma_{1}\ldots\gamma_{n}\mu_{1}\ldots\mu_{k}}\left({\widehat{\mathbf{L}}}_{ij}\right)\,,$
$\displaystyle\sigma_{kn}={\left(-1\right)^{n}}\frac{\left[2\left(k+n\right)-1\right]!!}{{k!\left(2n-1\right)!!}}\,,$
$\displaystyle\Lambda_{\gamma_{1}\ldots\gamma_{n}\mu_{1}\ldots\mu_{k}}\left({\widehat{\mathbf{L}}}_{ij}\right)=\mathop{\hat{L}_{\gamma_{1}}^{ij}{\ldots}\hat{L}_{\gamma_{n}}^{ij}\hat{L}_{\mu_{1}}^{ij}{\ldots}\hat{L}_{\mu_{k}}^{ij}}\limits^{\rule[-2.8903pt]{0.72229pt}{3.61371pt}\rule{75.16031pt}{0.72229pt}\rule[-2.8903pt]{0.72229pt}{3.61371pt}}\,,\quad\hat{L}_{\mu_{m}}^{ij}=L_{\mu_{m}}^{ij}/L_{ij}\,,\quad
m\in\mathbb{N}\,$

wherein series (39) converges uniformly and absolutely if inequalities (23)
hold true.

###### Proof.

Rewriting Eq. (13) in the relevant local Cartesian coordinates of $j$-th and
$i$-th spheres we have

$\displaystyle
r_{\gamma_{\nu}}^{j}=r_{\gamma_{\nu}}^{i}+L_{\gamma_{\nu}}^{ij}\,,$ (40)

where $L_{\gamma_{\nu}}^{ij}$ are coordinates of vector ${\mathbf{L}}_{ij}$
($\gamma_{\nu}=\overline{1,3}$). Linear dependence (40) leads to the evident
relations $\partial_{r_{\gamma_{\nu}}^{j}}=\partial_{r_{\gamma_{\nu}}^{i}}$.
Using them in Eq. (A.1), one readily obtains

$\displaystyle r_{j}^{-\left(2n+1\right)}\mathop{r_{\gamma_{1}}^{j}{\ldots
r}_{\gamma_{n}}^{j}}\limits^{\rule[-2.8903pt]{0.72229pt}{3.61371pt}\rule{32.52127pt}{0.72229pt}\rule[-2.8903pt]{0.72229pt}{3.61371pt}}=\frac{\left(-1\right)^{n}}{\left(2n-1\right)!!}\partial_{r^{i}_{\gamma_{1}}}\ldots\partial_{r^{i}_{\gamma_{n}}}\left(\frac{1}{\left\|\mathbf{r}_{i}-\mathbf{L}_{ij}\right\|}\right)$

or, taking into account that
$\partial_{r_{\gamma_{\nu}}^{i}}=-\partial_{L_{\gamma_{\nu}}^{ij}}$, one has
the relation

$\displaystyle r_{j}^{-\left(2n+1\right)}\mathop{r_{\gamma_{1}}^{j}{\ldots
r}_{\gamma_{n}}^{j}}\limits^{\rule[-2.8903pt]{0.72229pt}{3.61371pt}\rule{32.52127pt}{0.72229pt}\rule[-2.8903pt]{0.72229pt}{3.61371pt}}=\frac{1}{\left(2n-1\right)!!}\partial_{L_{\gamma_{1}}^{ij}}\ldots\partial_{L_{\gamma_{n}}^{ij}}\left(\frac{1}{\left\|\mathbf{r}_{i}-\mathbf{L}_{ij}\right\|}\right)\,.$
(41)

Theorem VII.1 yields

$\displaystyle\frac{1}{\|\mathbf{r}_{i}-\mathbf{L}_{ij}\|}=\sum\limits_{k=0}^{\infty}\frac{\left(-1\right)^{k}}{k!}\partial_{L_{\gamma_{1}}^{ij}}\ldots\partial_{L_{\gamma_{k}}^{ij}}\left(\frac{1}{L_{ij}}\right)\odot^{k}\mathop{r_{\gamma_{1}}^{i}{\ldots
r}_{\gamma_{k}}^{i}}\limits^{\rule[-2.8903pt]{0.72229pt}{3.61371pt}\rule{36.135pt}{0.72229pt}\rule[-2.8903pt]{0.72229pt}{3.61371pt}}\,.$
(42)

Substituting this expansion into Eq. (41) and taking into consideration the
uniform convergence of the series (42) we arrive at the required translation
addition theorem (39). ∎

Note that, following Miller, the transformation matrix
$U^{ij}_{\gamma_{1}\ldots\gamma_{n}\mu_{1}\ldots\mu_{k}}$ is commonly termed
the mixed-basis matrix Miller (1977).

For further investigation of the problem (25)-(27) it is expedient to recast
expression (39) in a dimensionless form

$\displaystyle\xi_{j}^{-\left(2n+1\right)}\mathop{{\xi}_{\gamma_{1}}^{j}{\ldots\xi}_{\gamma_{n}}^{j}}\limits^{\rule[-2.8903pt]{0.72229pt}{3.61371pt}\rule{36.135pt}{0.72229pt}\rule[-2.8903pt]{0.72229pt}{3.61371pt}}=\sum\limits_{k=0}^{\infty}\widetilde{U}^{ij}_{\gamma_{1}\ldots\gamma_{n}\mu_{1}\ldots\mu_{k}}\odot^{k}\mathop{{\xi}_{\mu_{1}}^{i}{\ldots\xi}_{{\mu_{k}}}^{i}}\limits^{\rule[-2.8903pt]{0.72229pt}{3.61371pt}\rule{36.135pt}{0.72229pt}\rule[-2.8903pt]{0.72229pt}{3.61371pt}}\,,$
(43)

where dimensionless mixed-basis matrix elements are

$\widetilde{U}^{ij}_{\gamma_{1}\ldots\gamma_{n}\mu_{1}\ldots\mu_{k}}:=\sigma_{kn}\varepsilon_{ij}^{n+1}\varepsilon_{ji}^{k}\Lambda_{\gamma_{1}\ldots\gamma_{n}\mu_{1}\ldots\mu_{k}}\left({\widehat{\mathbf{L}}}_{ij}\right)\,.$

Particularly dimensionless form of the degenerate translation addition theorem
(36) reads

$\displaystyle\xi_{j}^{-1}=\sum\limits_{k=0}^{\infty}\widetilde{U}^{ij}_{\mu_{1}\ldots\mu_{k}}\odot^{k}\mathop{{\xi}_{\mu_{1}}^{i}{\ldots\xi}_{{\mu_{k}}}^{i}}\limits^{\rule[-2.8903pt]{0.72229pt}{3.61371pt}\rule{36.135pt}{0.72229pt}\rule[-2.8903pt]{0.72229pt}{3.61371pt}}\,.$
(44)

The above presented proof of the I$\to$R TAT in terms of ICT seems to be the
simplest among previously known theorems for solid harmonics in terms of polar
spherical coordinates Caola (1978); Weniger and Steinborn (1985). It is
important to stress that using connection between spherical solid harmonics
and ICT (see Appendix LABEL:sec:appendixCon) the proved addition theorem (39)
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 Traytak and Price (2007), when intrinsic rates are functions on the
angular local coordinates $\kappa_{i}\left(\theta_{i},\phi_{i}\right)$, we
need the rotational addition theorem for the ICT. The latter theorem may be
proved in entirely similar way.

##  VIII Solution to the problem

In this section, we implement the GMSV algorithm described in Sec. V using ICT
technique. Moreover we shall develop here the version of the GMSV elaborated
previously in Refs. Traytak (1991, 1992, 2001).

###  VIII.1 The GMSV algorithm implementation

For clarity sake we organize our treatment in according to the method outline
V.3 distinguishing 8 steps inherent the GMSV.

Step 1. The mixed boundary value problem (25)-(27) is well posed and has a
unique classical solution Courant and Hilbert (1953); Tikhonov and Samarskii
(1963); Krutitskii (1999).

Step 2. First, we can apply superposition principle. Decompose the general
solution $u_{N}\left(P\right)$ of above diffusion problem (25)-(27) with
respect to the $i$-th sink, considering the set of partial solutions to Eq.
(25), which are defined on the charts:
$\left\\{u_{i}\left({\bm{\xi}}_{i}\right)\right\\}_{i=1}^{N}$
($u_{i}:\Omega_{i}^{-}\to\left(0,1\right]$). Thus, in the local Cartesian
coordinates $\\{O_{i};{\bm{\xi}}_{i}\\}$ we can write

$\displaystyle
u_{N}\left({\bm{\xi}}_{i}\right)=u_{i}\left({\bm{\xi}}_{i}\right)+\sum\limits_{j\left(\not=i\right)=1}^{N}u_{j}\left({\bm{\xi}}_{j}\right)\quad\text{in}\quad\Omega^{-}_{i}\,.$
(45)

Here, by virtue of Eqs. (13) and (23), one has evident relations:
$\varepsilon_{ij}{\bm{\xi}}_{j}=\varepsilon_{ji}{\bm{\xi}}_{i}+{\widehat{\mathbf{L}}}_{ij}$
(see Fig. 1). From the physical standpoint representation (45) 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}\left(P\right)$ in a vicinity of this sink.

Clearly, by definition functions $u_{i}\left({\bm{\xi}}_{i}\right)$ are
harmonic in charts $\Omega^{-}_{i}$, i.e.

$\displaystyle{{\bm{\nabla}}}_{{\bm{\xi}}_{i}}^{2}u_{i}\left({\bm{\xi}}_{i}\right)=0\,,$
(46)
$\displaystyle\left.u_{i}\left({\bm{\xi}}_{i}\right)\right|_{\xi_{i}\rightarrow\infty}\rightrightarrows
0\,.$ (47)

Therefore, the original diffusion boundary value problem (25)-(27) in the
manifold $\Omega^{-}$ (intersection of all charts $\Omega^{-}_{i}$) is reduced
to $N$ coupled problems for $u_{i}\left({\bm{\xi}}_{i}\right)$ in simpler
charts $\Omega^{-}_{i}$, with the advantages of the local Cartesian
coordinates $\\{O_{i};{\bm{\xi}}_{i}\\}$.

Alternatively, for inhomogeneous partially reflecting boundary conditions (26)
in an exterior $\epsilon$-neighborhood of the boundary $\partial\Omega_{i}$:
$\Omega^{-}_{i}\left(\epsilon\right):=\\{{\bm{\xi}}_{i}:1<\xi_{i}<1+\epsilon\\}\subset\Omega^{-}_{i}$
it is convenent to recast solution $u_{N}\left({\bm{\xi}}_{i}\right)$ as

$\displaystyle
u_{N}\left({\bm{\xi}}_{i}\right)=u_{\left(0\right)}^{i}{\left(\xi_{i}\right)}+\delta
u_{i}\left({\bm{\xi}}_{i}\right)\,\quad\text{in}\quad\Omega^{-}_{i}\left(\epsilon\right)\,,$
(48)

where $u_{\left(0\right)}^{i}{\left(\xi_{i}\right)}$ is the unperturbed one-
sink solution and $\delta u_{i}\left({\bm{\xi}}_{i}\right)$ is a perturbation
due to influence of each from other $N-1$ sinks for
${j\left(\not=i\right)=\overline{1,N}}$. Owing to Eq. (26) we assume that
these solutions obey inhomogeneous and homogeneous partially reflecting
boundary conditions, respectively

$\displaystyle-\left.\left(\partial_{\xi_{i}}u_{\left(0\right)}^{i}-\kappa_{i}u_{\left(0\right)}^{i}\right)\right|_{\partial\Omega_{i}^{1}}=\kappa_{i}\,,$
(49) $\displaystyle\left.\left(\partial_{\xi_{i}}\delta u_{i}-\kappa_{i}\delta
u_{i}\right)\right|_{\partial\Omega_{i}^{1}}=0\,.$ (50)

Therewith the unperturbed by diffusive interaction solution

$\displaystyle
u_{\left(0\right)}^{i}{\left(\xi_{i}\right)}=J_{0}^{i}{\xi_{i}}^{-1}\,$ (51)

obeys inhomogeneous boundary condition (49). Here $J_{0}^{i}$ is the $i$-th
Collins-Kimball rate (31).

One important point to emphasis is that Eqs. (45) and (48) are global and
local representations of the same function defined on the manifold:
$u_{N}:\Omega^{-}\to\left(0,1\right]$.

Step 3. It can be shown that tensors fields introduced by Eqs. (A.5) and (A.6)
are families of irregular
$\\{\mathbf{X}_{n}^{-}\left({\bm{\xi}}_{i}\right)\\}_{n=0}^{\infty}$ and
regular $\\{\mathbf{X}_{n}^{+}\left({\bm{\xi}}_{i}\right)\\}_{n=0}^{\infty}$
solid Cartesian harmonics, which form desired basis functions (a complete and
orthogonal sets) in domains $\Omega^{-}_{i}$ and $\Omega^{+}_{i}$,
respectively. Thus, we can represent a partial solution
$u_{i}\left({\bm{\xi}}_{i}\right)$ with the help of the irregular ICT,
satisfying regularity condition at infinity

$\displaystyle
u_{i}\left({\bm{\xi}}_{i}\right)=\sum\limits_{n=0}^{\infty}\xi_{i}^{-\left(2n+1\right)}A_{\gamma_{1}...\gamma_{n}}^{i}\odot^{n}\mathop{{\xi}_{\gamma_{1}}^{i}{\ldots\xi}_{\gamma_{n}}^{i}}\limits^{\rule[-2.8903pt]{0.72229pt}{3.61371pt}\rule{36.135pt}{0.72229pt}\rule[-2.8903pt]{0.72229pt}{3.61371pt}}\quad\mbox{in}\quad\Omega^{-}_{i}\,.$
(52)

On the other hand, in an exterior neighborhood of $i$-th sink boundary
$\partial\Omega_{i}$ we can represent the perturbation function in Eq. (45) as
a series with respect to the regular ICT

$\displaystyle\sum\limits_{j\left(\not=i\right)=1}^{N}u_{j}\left({\bm{\xi}}_{j}\right)=\sum\limits_{n=0}^{\infty}B_{\gamma_{1}\ldots\gamma_{n}}^{i}\odot^{n}\mathop{{\xi}_{\gamma_{1}}^{i}{\ldots\xi}_{\gamma_{n}}^{i}}\limits^{\rule[-2.8903pt]{0.72229pt}{3.61371pt}\rule{33.96698pt}{0.72229pt}\rule[-2.8903pt]{0.72229pt}{3.61371pt}}\,.$
(53)

In Eqs. (52) and (53) $A_{\gamma_{1}\ldots\gamma_{n}}^{i}$ and
$B_{\gamma_{1}\ldots\gamma_{n}}^{i}$ are tensor coefficients to be determined
from the boundary conditions (26).

Step 4. Theorem VII.2 in the dimensionless form (43) gives the required
translation addition theorem to tackle the diffusion problem (25)-(27).

Step 5. Now, by means of dimensionless translation addition theorem (43), we
satisfy the Robin boundary conditions (26). Hence, taking advantage of the ICT
linear independence, we can express the coefficients
$A_{\gamma_{1}\ldots\gamma_{n}}^{i}$ in terms of
$B_{\gamma_{1}\ldots\gamma_{n}}^{i}$

$\displaystyle A_{0}^{i}=J_{0}^{i}\left(1-B_{0}^{i}\right)\quad\mbox{for}\quad
n=0\,,$ (54) $\displaystyle
A_{\gamma_{1}\ldots\gamma_{n}}^{i}=\frac{n-\kappa_{i}}{1+n+\kappa_{i}}B_{\gamma_{1}\ldots\gamma_{n}}^{i}\quad\mbox{for}\quad
n\in{\mathbb{N}}\,.$ (55)

With the aid of these connections formula (48) for the perturbed diffusion
field in an exterior neighborhood of the surface $\partial\Omega_{i}$
($i=\overline{1,N}$) yields

$\displaystyle\delta
u_{i}\left({\bm{\xi}}_{i}\right)=\sum\limits_{k=0}^{\infty}\left[1+\left(\frac{k-\kappa_{i}}{1+k+\kappa_{i}}\right)\xi_{i}^{-\left(2k+1\right)}\right]B_{\gamma_{1}\ldots\gamma_{k}}^{i}\odot^{k}\mathop{{\xi}_{\gamma_{1}}^{i}{\ldots\xi}_{\gamma_{k}}^{i}}\limits^{\rule[-2.8903pt]{0.72229pt}{3.61371pt}\rule{33.96698pt}{0.72229pt}\rule[-2.8903pt]{0.72229pt}{3.61371pt}}\,.$
(56)

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

And, in addition, using Eq. (52) under relationships (54) and (55) for domaind
$\Omega^{-}_{j}$ we immediately write down an auxiliary relation

$\displaystyle\sum\limits_{j\left(\not=i\right)=1}^{N}u_{j}\left({\bm{\xi}}_{j}\right)=\sum\limits_{j\left(\not=i\right)=1}^{N}u_{\left(0\right)}^{j}{\left(\xi_{j}\right)}$
$\displaystyle+\sum\limits_{j\left(\not=i\right)=1}^{N}\sum\limits_{k=0}^{\infty}\left(\frac{k-\kappa_{j}}{1+k+\kappa_{j}}\right)\xi_{j}^{-\left(2k+1\right)}B_{\gamma_{1}\ldots\gamma_{k}}^{j}\odot^{k}\mathop{{\xi}_{\gamma_{1}}^{j}{\ldots\xi}_{\gamma_{k}}^{j}}\limits^{\rule[-2.8903pt]{0.72229pt}{3.61371pt}\rule{32.52127pt}{0.72229pt}\rule[-2.8903pt]{0.72229pt}{3.61371pt}}\,.\qquad$
(57)

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

Step 6. Now in the local Cartesian coordinates $\\{O_{i};{\bm{\xi}}_{i}\\}$
one can carry out a self-consistent procedure to find unknown tensor
coefficients $B_{\gamma_{1}\ldots\gamma_{k}}^{i}$ appearing in formulas (56)
and (57). First we apply the dimensionless translation addition theorem (43)
including its degenerate form (44) to the right-hand side of Eq. (57) in order
to recast it in the local coordinates $\\{O_{i};{\bm{\xi}}_{i}\\}$. Then with
allowance made for expression (53) and orthogonality property of the ICT on
the unit spheres $\partial\Omega_{i}^{1}$ (A.9) we derive the required self-
consistent infinite system of linear algebraic equations (ISLAE) of the II
kind with respect to unknown tensor coefficients
$B_{\gamma_{1}\ldots\gamma_{k}}^{i}$. We should note that these infinite
systems are commonly referred to as the resolving ISLAE Kushch (2013); Traytak
(2020).

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

(a) For partialy reflecting sinks (when $0<\kappa_{i}<\infty$) the resolving
ISLAE is

$\displaystyle
B_{\gamma_{1}\ldots\gamma_{k}}^{i}=B_{\gamma_{1}\ldots\gamma_{k}}^{i\left(0\right)}+\sum\limits_{j\left(\neq
i\right)=1}^{N}\sum\limits_{l=0}^{\infty}\left(\frac{k-\kappa_{j}}{1+k+\kappa_{j}}\right)$
$\displaystyle\times\widetilde{U}^{ij}_{\gamma_{1}\ldots\gamma_{k}\mu_{1}\ldots\mu_{l}}\odot^{l}B_{\mu_{1}\ldots\mu_{l}}^{j}\,,\quad
k=\overline{0,\infty}\,.$ (58)

Hereinafter for short we denoted

$\displaystyle
B_{\gamma_{1}\ldots\gamma_{k}}^{i\left(0\right)}=\sum\limits_{j\left(\neq
i\right)=1}^{N}J_{0}^{j}\widetilde{U}^{ij}_{\gamma_{1}\ldots\gamma_{k}}\,.$
(59)

Another two important ISLAE are special cases of (58).

(b) In case of fully absorbing sinks (when $\kappa_{i}\to\infty$) the
resolving ISLAE is simplified to

$\displaystyle
B_{\gamma_{1}\ldots\gamma_{k}}^{i}=B_{\gamma_{1}\ldots\gamma_{k}}^{i\left(0\right)}-\sum\limits_{j\left(\neq
i\right)=1}^{N}\sum\limits_{l=0}^{\infty}\widetilde{U}^{ij}_{\gamma_{1}\ldots\gamma_{k}\mu_{1}\ldots\mu_{l}}\odot^{l}B_{\mu_{1}\ldots\mu_{l}}^{j}\,,$
(60)

where $k=\overline{0,\infty}$.

(c) For fully reflecting (when $\kappa_{i}\to 0$) the resolving ISLAE (58)
obviously yields

$\displaystyle
B_{\gamma_{1}\ldots\gamma_{k}}^{i}=\left(\frac{k}{1+k}\right)\sum\limits_{j\left(\neq
i\right)=1}^{N}\sum\limits_{l=0}^{\infty}\widetilde{U}^{ij}_{\gamma_{1}\ldots\gamma_{k}\mu_{1}\ldots\mu_{l}}\odot^{l}B_{\mu_{1}\ldots\mu_{l}}^{j}\,,$
(61)

where $k=\overline{1,\infty}$.

Step 7. Consider here the most general case of arrays with partialy reflecting
sinks. In other words to find $B_{\gamma_{1}\ldots\gamma_{k}}^{i}$ the system
(58) requires the inversion of an corresponding infinite-dimensional matrix
Galanti  _et al._ (2016, 2020). So we deal with a typical problem which may be
solved by methods of functional analysis Kantorovich and Krylov (1964);
Kantorovich and Akilov (1982). For instance, provided the ISLAE (58) is of
normal Poincare-Koch type it implies that the Fredholm-Hilbert alternative
holds true Kantorovich and Krylov (1964). Hence, a unique solution of (58)
exists and it may be found by the method of reduction to any degree of
accuracy Traytak (2020). The latter means that the resolving ISLAE (58) should
be truncated and then inverted obtained finite system numerically to get
approximations of the coefficients $B_{\gamma_{1}\ldots\gamma_{k}}^{i}$. Under
some conditions the resolving ISLAE (58) possesses solution which may be
obtained by means of simple iterations Galanti  _et al._ (2016) (see also
discussion below in Sec. IX). Obviously the zero and the first iterations for
$B_{0}^{i}$ are

$\displaystyle B_{0}^{i\left(0\right)}=\sum\limits_{j\left(\neq
i\right)=1}^{N}J_{0}^{j}\widetilde{U}^{ij}_{0}\,,$ (62) $\displaystyle
B_{0}^{i(1)}=B_{0}^{i\left(0\right)}-\sum\limits_{j\left(\neq
i\right)=1}^{N}J_{0}^{j}\widetilde{U}^{ij}_{0}B_{0}^{j\left(0\right)}-\sum\limits_{j\left(\neq
i\right)=1}^{N}\sum\limits_{l=1}^{\infty}J_{0}^{j}\widetilde{U}^{ij}_{\mu_{1}\ldots\mu_{l}}\odot^{l}B_{\mu_{1}\ldots\mu_{l}}^{j\left(0\right)}\,.$
(63)

Here and below $\widetilde{U}^{ij}_{0}$ is the monopole approximation to the
mixed-basis matrix. Interestingly that in the course of solving the ISLAE (58)
by simple iterations angular dependence arises starting from terms
${\mathcal{O}}\left(\varepsilon^{4}\right)$.

Step 8. Once the local concentration $u_{N}\left(P\right)$ 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 (48) of function
$u_{N}\left(P\right)$. Due to orthogonality condition for the ICT (A.9)
explicit surface integration on the unit sphere in the general formula (29)
will yield the desired rates (30).

######  Remark VIII.1.

Quite apparently, functions
$\mathop{{\hat{\xi}}_{\gamma_{1}}^{i}\ldots{\hat{\xi}}_{\gamma_{n}}^{i}}\limits^{\rule[-2.8903pt]{0.72229pt}{3.61371pt}\rule{33.96698pt}{0.72229pt}\rule[-2.8903pt]{0.72229pt}{3.61371pt}}$
are harmonics polynomianls defined on the $i$-th unit sphere
$\partial\Omega_{i}^{1}$ and, in accordance with the Weierstrass approximation
theorem, any function $w_{i}{\left({\hat{\bm{\xi}}}_{i}\right)}\in
C(\partial\Omega_{i}^{1})$ may be expanded in the series

$\displaystyle
w_{i}\left({{\hat{\bm{\xi}}}}_{i}\right)=\sum\limits_{n=0}^{\infty}\alpha_{\gamma_{1}\ldots\gamma_{n}}^{i}\odot^{n}\mathop{{\hat{\xi}}_{\gamma_{1}}^{i}\ldots{\hat{\xi}}_{\gamma_{n}}^{i}}\limits^{\rule[-2.8903pt]{0.72229pt}{3.61371pt}\rule{32.52127pt}{0.72229pt}\rule[-2.8903pt]{0.72229pt}{3.61371pt}}\,,$
(64)

where $\alpha_{\gamma_{1}\ldots\gamma_{n}}^{i}$ are some constant tensor
coefficients to be determined by means of the orthogonality condition (A.9).
One can see that expression (64) makes it possible to apply the ICT method
being considered to the diffusion problems under general Robin boundary
conditions (18).

###  VIII.2 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
$\kappa^{(N)}$ we have proved the following assertions.

(a) On global representation of the solution $u_{N}\left(P\right)$.

######  Theorem VIII.1.

Using the local Cartesian coordinates $\\{O_{j};{\bm{\xi}}_{j}\\}$
($j=\overline{1,N}$) the solution $u_{N}\left(P\right)$ ($P\in\Omega^{-}$, see
Fig. 1) to the exterior Robin boundary value problem (25)-(27) may be
repesented as follows:

$\displaystyle u_{N}\left(P\right)=\sum\limits_{j=1}^{N}$
$\displaystyle\left[J_{0}^{j}\left(1-B_{0}^{j}\right){\xi_{j}}^{-1}+\sum\limits_{k=0}^{\infty}\left(\frac{k-\kappa_{j}}{1+k+\kappa_{j}}\right)\right.$
$\displaystyle\left.\quad\times\xi_{j}^{-\left(2k+1\right)}B_{\gamma_{1}\ldots\gamma_{k}}^{j}\odot^{k}\mathop{{\xi}_{\gamma_{1}}^{j}{\ldots\xi}_{\gamma_{k}}^{j}}\limits^{\rule[-2.8903pt]{0.72229pt}{3.61371pt}\rule{32.52127pt}{0.72229pt}\rule[-2.8903pt]{0.72229pt}{3.61371pt}}\right]\,,$
(65)

where tensor coefficients $B_{\gamma_{1}\ldots\gamma_{k}}^{i}$
($k=\overline{0,\infty}$) are solution of the resolving ISLAE (58).

(b) On local representation of the solution $u_{N}\left(P\right)$.

######  Theorem VIII.2.

Solution $u_{N}\left({\bm{\xi}}_{i}\right)$ to the exterior Robin boundary
value problem (25)-(27) in a small exterior $\epsilon$-neighborhood of the
unit sphere $\partial\Omega_{i}^{1}$ with respect to the local Cartesian
coordinates $\\{O_{i};{\bm{\xi}}_{i}\\}$ reads

$\displaystyle
u_{N}\left({\bm{\xi}}_{i}\right)=J_{0}^{i}{\xi_{i}}^{-1}+\sum\limits_{k=0}^{\infty}\left[1+\left(\frac{k-\kappa_{i}}{1+k+\kappa_{i}}\right)\xi_{i}^{-\left(2k+1\right)}\right]$
$\displaystyle\times
B_{\gamma_{1}\ldots\gamma_{k}}^{i}\odot^{k}\mathop{{\xi}_{\gamma_{1}}^{i}{\ldots\xi}_{\gamma_{k}}^{i}}\limits^{\rule[-2.8903pt]{0.72229pt}{3.61371pt}\rule{33.96698pt}{0.72229pt}\rule[-2.8903pt]{0.72229pt}{3.61371pt}}\quad\mbox{in}\quad\Omega^{-}_{i}\left(\epsilon\right)\quad
i=\overline{1,N}\,;\quad$ (66)

where tensor coefficients $B_{\gamma_{1}\ldots\gamma_{k}}^{i}$
($k=\overline{0,\infty}$) are solution of the resolving ISLAE (58).

So, if we know the solution $u_{N}\left({\bm{\xi}}_{i}\right)$ in an exterior
neighborhood adjacent to the $i$-th reaction surface (66), we can calculate
the $B$-particles flux on the $i$-th sink and, therefore, the corresponding
absorbing rate and, in this way, screening coefficient (30). Thus, taking
advantage of Theorem VIII.2, one can readily prove an important

######  Corollary 1.

For the exterior Robin boundary value problem (25)-(27) the screening
coefficient in Eq. (30) has the form

$\displaystyle J_{i}=1-B_{0}^{i}\quad\mbox{and}\quad 0\leq B_{0}^{i}<1\,,$
(67)

where coefficients $B_{0}^{i}$ are determined by means of solution to the
resolving ISLAE (58).

It is worth noting here that in Ref. Krasovitov (1991) Krasovitov thoroughly
reproducing the notation, derivations and even misprints of our paper Traytak
(1991). However, he has been the first who kept all terms comprising
$B_{\gamma_{1}...\gamma_{n}}^{i}$ ($n\in\mathbb{N}$) in the series with
respect to orthogonal system of ICT after integration over the unit sphere
(29) (see detailed discussion of these points in Ref. Traytak (2001)).

###  VIII.3 Monopole approximation

As we already noted above in Sec. III 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

######  Definition VIII.1.

Monopole approximation corresponds to truncations of series (52) and (53) for
any partial solutions $u_{i}\left({\bm{\xi}}_{i}\right)$ ($i=\overline{1,N}$)
such that: $A_{0}^{i},B_{0}^{i}\neq 0$ and
$A_{\gamma_{1}\ldots\gamma_{k}}^{i},B_{\gamma_{1}\ldots\gamma_{k}}^{i}\equiv
0$ for all $k\in\mathbb{N}$.

It is notable that one can see directly from the resolving ISLAE (61) 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.
(67) 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 (60) one infers that for the monopole approximation, values
of the total flux $J_{i}$ may be determined by the much simpler finite
$N\times N$ system of equations:

$\displaystyle
J_{i}=1-\sum\limits_{j=1}^{N}\widetilde{V}^{ij}_{0}J_{j}\,,\quad
i=\overline{1,N}\,;$ (68)
$\displaystyle\widetilde{V}^{ij}_{0}:=\left\\{\begin{array}[]{ll}\widetilde{U}^{ij}_{0}&\quad\mbox{if}\quad
j\neq i\,,\\\ 0&\quad\mbox{if}\quad i=j\,.\end{array}\right.$ (71)

Therefore, expressions (45) and (67) are simplified to

$\displaystyle
u_{N}=u_{\left(0\right)}^{i}+B_{0}^{i}\left(1-\frac{1}{\xi_{i}}\right)\,,$
(72) $\displaystyle J_{i}=1-B_{0}^{i}\,.$ (73)

In the important particular case at $N=2$ system (68) gives

$\displaystyle J_{i}=1-\varepsilon_{ij}J_{j}\,,\quad(i\neq j)=1,2$ (74)

with the first iteration

$\displaystyle J_{i}=1-\varepsilon_{ij}\,,\quad(i\neq j)=1,2\,.$ (75)

It is evident that system (74) has the exact solution

$\displaystyle
J_{i}=\frac{1-\varepsilon_{ij}}{1-\varepsilon_{ji}\varepsilon_{ij}}\,,\quad(i\neq
j)=1,2\,.$ (76)

Note in passing that sometimes in literature Eq. (76) is called the monopole
approximation, however, the monopole approximation is given by Eq. (75).
Besides Eq. (76) may be also obtained by the reflection method for the
monopole approximation Traytak (1992).

###  VIII.4 General remarks on the obtained solution

In Sec. III we already stated that exact solution to the diffusion boundary
value problem (25)-(27) is possible only for two fully absorbing sinks by
means of bispherical coordinates Carstens, Williams, and Zung (1970); Samson
and Deutch (1977); Dubinko  _et al._ (1989); Bluett and Green (2006). Using
the bispherical coordinates for the problem under partially reflecting
boundary conditions (26) 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}\left(P\right)$ from
the boundary value problem (25)-(27) in general case of $N$ sinks
microstructure (11), when $N>2$ is a far too complicated task to be solved
analytically in closed form.

A comparison of the general resolving ISLAE (58) with that (60) 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_{\left(0\right)}^{i}{\left(\xi_{i}\right)}$ (51) is
influenced by a contribution mediated by all other $N-1$ sinks (66), leading,
in it turn, to the screening effects similar to that in an electrostatics
Deutch, Felderhof, and Saxton (1976).

We emphasize again that discussed in Sec. IV 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}$ (67) for the Dirichlet-Robin proper mixed
problem at $N=3$ below in Sec. X.

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 Bénichou
_et al._ (2005).

Finally, it is worth noting that stationary diffusion equation under so-called
conjugate boundary conditions Mityushev and NawalaniecX (2019); Grebenkov and
Traytak (2019); Traytak (2020) may be also tackled with above method quite
similarly.

##  IX 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 (58). Alternative numerical solution of the
aforementioned truncated system of equations corresponding to (58) 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 (25)-(27)
for the case of fully diffusion-controlled reactions.

###  IX.1 Relations with method of reflections

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

$\displaystyle\max_{P\in\overline{\Omega}^{-}}|u_{N}\left(P\right)|\leq\max_{P\in\partial\Omega^{-}}|u_{N}\left(P\right)|\,.$
(77)

Let a harmonic function $v_{N}\left(P\right)$ be an approximation to the exact
solution $u_{N}\left(P\right)$ of the diffusion problem (25)-(27). It
immediately follows from inequality (77) that the maximum absolute error of
the residual for any approximation $|v_{N}\left(P\right)-u_{N}\left(P\right)|$
is achieved on the boundary of the configuration manifold
$\partial\Omega^{-}$. 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 ${\Omega}^{-}$.

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) Smirnov (1964); Traytak (2006); Laurent, Legendre, and
Salomon (2021). 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 Traytak (2006). 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 Kim (1987). We do not intend to discuss
this subject here and so we refer the interested reader to Ref. Laurent,
Legendre, and Salomon (2021).

According to the method of reflections one shold construct a functional
sequence of relevant approximations
$\left\\{v_{m}\left(P\right)\right\\}_{m=0}^{\infty}$ so-called reflections
(see p. 638 of Ref. Smirnov (1964) and Ref. Traytak (2006) for details) such
that this sequence converges in $\Omega^{-}$ uniformly to the exact solution

$\displaystyle v_{m}\left(P\right)\rightrightarrows
u_{N}\left(P\right)\quad\text{as}\quad m\to\infty\,.$ (78)

In Ref. Traytak (2006) we proved the following theorem on the conditions for
uniformly convergence of the method of reflections (78).

######  Theorem IX.1.

For the relationship (78) to hold the following sufficient (a) and necessary
(b) conditions must be satisfied

$\displaystyle\mbox{(a)}\quad
M_{N}:=\max_{i=\overline{1,N}}\sum\limits_{k\left(\neq
i\right)=1}^{N}\left(\frac{\varepsilon_{ik}}{1-\varepsilon_{ki}}\right)<1\,,$
(79) $\displaystyle\mbox{(b)}\quad
m_{N}:=\min_{i=\overline{1,N}}\sum\limits_{k\left(\neq
i\right)=1}^{N}\left(\frac{\varepsilon_{ik}}{1+\varepsilon_{ki}}\right)<1\,.$
(80)

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
Laurent, Legendre, and Salomon (2021).

Plainly, Theorem IX.1 implies the following two assertions.

######  Corollary 2.

For sinks of identical radius $R_{i}=R$ ($i=\overline{1,N}$) conditions (79)
and (80) are simplified to

$\displaystyle\mbox{(a)}\quad
M_{N}=\max_{i=\overline{1,N}}\sum\limits_{k\left(\neq
i\right)=1}^{N}\left(\frac{L_{ik}}{R}-1\right)^{-1}<1\,,$ (81)
$\displaystyle\mbox{(b)}\quad
m_{N}=\min_{i=\overline{1,N}}\sum\limits_{k\left(\neq
i\right)=1}^{N}\left(\frac{L_{ik}}{R}+1\right)^{-1}<1\,.$ (82)

######  Corollary 3.

For the monopole approximation conditions (79) and (80) are reduced to

$\displaystyle\mbox{(a)}\quad
M_{N}=\max_{i=\overline{1,N}}\sum\limits_{k\left(\neq
i\right)=1}^{N}\varepsilon_{ik}<1\,,\quad\mbox{(b)}\quad
m_{N}=\min_{i=\overline{1,N}}\sum\limits_{k\left(\neq
i\right)=1}^{N}\varepsilon_{ik}<1\,.$ (83)

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

$\displaystyle\|\widetilde{V}^{ij}_{0}\|_{\infty}:=\max_{i=\overline{1,N}}\sum_{k=1}^{N}\widetilde{V}^{ik}_{0}<1\,,$
(84)

where $\widetilde{V}^{ij}_{0}$ is the hollow matrix (71).

Moreover, previously in Ref. Traytak (2003) we have shown that this approach
is tightly connected with the corresponding iterative solution within the
scope of the GMSV and proved

######  Theorem IX.2.

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

In this way Theorem IX.2 immediately leads to the following remarkable result

######  Corollary 4.

Conditions (79) and (80) may be used to estimate the convergence of the
simple-iteration method when applied to the resolving ISLAE (60).

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 Galanti  _et
al._ (2016).

###  IX.2 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
(79) is satisfied and it yields inequality
$M_{2}=R_{1}/\left(1-R_{2}\right)<1$. Hence we have
$h=L-\left(R_{1}+R_{2}\right)>0$, where $h$ is the distance between sink
surfaces. In its turn, the latter inequality coincides with non-overlapping
sinks condition (14).

Then, by Theorem IX.1 and Corollary 4 we arrive at the analog of the known
Goluzin theorem Smirnov (1964).

######  Theorem IX.3.

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

To show how one can directly apply the Corollary 4, as examples consider
arrays of $N$ fully absorbing sinks of identical radius $R_{i}=R$
($i=\overline{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. 4 (a)). We have
performed calculations at $N=100$. It turned out that the simple-iteration
method using to solve the resolving ISLAE (60) is not valid when $R/a<4.9$,
and when $R/a<9.4$ its applicability requires additional justification.

Figure 4: The array geometry of fully absorbing sinks located: (a) along a
straight line; (b) at the apices of an equilateral triangle

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 Deutch, Felderhof, and Saxton (1976); Labowsky (1978); Traytak (1992);
Eun, Kekenes-Huskey, and McCammon (2013).

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. 4 (b)). One can easy verify that sufficient condition for the
convergence of the simple-iteration method (81) fails when $1/3<R/a<1/2$. This
is entirely consistent with the conclusion by Kim in Ref. Kim (1987), where
were found numeically that for three sphere system relevant reflections series
will diverge when $0.46\lesssim 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. Mityushev and NawalaniecX (2019); Laurent, Legendre, and
Salomon (2021).

##  X 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.

###  X.1 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.

Figure 5: Cross section through the level surfaces of $u_{2}\left(P\right)$
for decrement step $\vartriangle h=0.05$ ($l_{0}=20$): (a) radius
$R^{\ast}_{1}=3$ and the sink spacing $L^{\ast}=5$; (b) radius
$R^{\ast}_{1}=10$ and the sink spacing $L^{\ast}=12$ Figure 6: Cross section
through the level surfaces of $u_{2}\left(P\right)$ for radius
$R^{\ast}_{1}=100$ and the sink spacing $L^{\ast}=111$: (a) decrement step
$\vartriangle h=0.05$ ($l_{0}=20$); (b) decrement step $\vartriangle h=0.01$
($l_{0}=100$).

By means of Theorem VIII.1 the local concentration field may be provided in
the entire manifold $\Omega^{-}$. For this purpose let us introduce the
concentration level surfaces: $S\left(h\right)\subset\Omega^{-}$ such that
$\left.u_{2}\left(P\right)\right|_{S}=h\in\left(0,1\right)$. As usually for
discretization we perform uniform partition of the concentration range
$\left(0,1\right)$ with a fixed step value $\vartriangle h$, such that $1$
correspond to the reaction surfaces and $0$ to infinity: $h_{l}=l\vartriangle
h$ for $l=\overline{0,l_{0}}$, where $l_{0}=1/\vartriangle 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}\left(P\right)$ on two dimensionless geometrical
parameters in units of the smaller sink: reduced sink 1 radius
$R^{\ast}_{1}:=R_{1}/R_{2}$ and the sink spacing $L^{\ast}=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. 5 and (a) of Fig.
6 we used decrement step $\vartriangle h=0.05$, whereas calculations in plot
(b) of Fig. 6 were performed at $\vartriangle h=0.01$.

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

######  Proposition X.1.

There exists a number $h_{e}\in\left(0,1\right)$ such that one can treat the
domain $\Omega_{e}:=\left\\{P:u_{N}\left(P\right)<h_{e}\right\\}$ as an
effective fully absorbing sink with respect to the scaled concentration
$u^{\ast}_{N}\left(P\right):=u_{N}\left(P\right)/h_{e}$ given on the
configuration manifold
$\Omega^{-}_{e}:=\left\\{P:u^{\ast}_{N}\left(P\right)>1\right\\}$.

For instance, at $N=2$ Fig. 6 shows that one can take $h_{e}=0.95$ to get
snowman-shaped effective fully absorbing sink with the configuration manifold
$\Omega_{e}^{-}=\left\\{P:u_{2}\left(P\right)>0.95\right\\}$.

###  X.2 Calculations of the screening coefficients

Within this subsection to calculate the screening coefficients due to
siffusive interaction we shall use the Corollary 1 of the Theorem VIII.2 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}\geq
J_{2}$ at all sink-sink center separations $L>R_{1}+R_{2}$, wherein $J_{1}\to
1$ as $R_{2}\to 0$. Furthermore, there is no diffusive interaction between
sinks and any number $N_{p}$ ($1\leq N_{p}\leq N$) of isolated points. Really
for the isolated point sinks $R_{j}=\varepsilon_{ij}=0$ and, therefore,
corresponding mixed-basis matrix elements in the resolving ISLAE (58) 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 (60)
converges up to the contact of sinks at $N=2$ Traytak (2006).

For simplicity sake consider here two equal sinks $R_{1}=R_{2}=R$, so we have:
$\varepsilon_{21}=\varepsilon_{12}=\epsilon\in(0,0.5)$ and for screening
coefficients
$J_{1}\left(\epsilon\right)=J_{2}\left(\epsilon\right)=J\left(\epsilon\right)$.
Solving the ISLAE (60) by simple iterations as $\epsilon\rightarrow 0+$ we
find:

$\displaystyle J\left(\epsilon\right)=J_{a}\left(\epsilon\right)+{\cal
O}\left({\epsilon}^{6}\right)\,,$ (85) $\displaystyle
J_{a}\left(\epsilon\right)=1-\epsilon+{\epsilon}^{2}-{\epsilon}^{3}+2{\epsilon}^{4}-3{\epsilon}^{5}\,.$
(86)

It is noteworthy that angular dependence in (86) arises starting from the
terms ${\cal O}\left({\epsilon}^{4}\right)$.

To show that $J_{a}\left(\epsilon\right)$ 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

$\delta_{m}:=\sup_{\epsilon\in(0,0.5)}\delta\left(\epsilon\right)=\lim_{\epsilon\rightarrow
0.5-}\delta\left(\epsilon\right)\,,$

where
$\delta\left(\epsilon\right):=|J\left(\epsilon\right)-J_{a}\left(\epsilon\right)|/J\left(\epsilon\right)$
is the relative error.

In the case under consideration with the aid of bispherical coordinate system
Samson and Deutch obtained the exact reaction rate Samson and Deutch (1977)

$\displaystyle
J\left(\epsilon\right)=\sum\limits_{n=1}^{\infty}\left(-1\right)^{1+n}\frac{\sinh\left(\mu_{0}\right)}{\sinh\left(n\mu_{0}\right)},$
(87)

where $\mu_{0}:=\cosh^{-1}\left(1/{2\epsilon}\right)$. For the contact, i.e.
when $\epsilon\rightarrow 0.5-$ (or $\mu_{0}\rightarrow 0+$) we have fully
absorbing dumbbell system and formula (87) yields the exact limit

$\lim_{\epsilon\rightarrow 0.5-}J\left(\epsilon\right)=\ln 2$

and, therefore, we get desired upper bound for the relative error

$\displaystyle\delta\left(\epsilon\right)<\delta_{m}=1-\frac{21}{32\ln
2}\approx 0.053\,.$ (88)

One can see that relative error decreases rather rapidly when ${\epsilon}$
decreases, e.g., $\delta\left(0.4\right)\approx 0.015$,
$\delta\left(0.3\right)\approx 0.003$ and besides
$J\left(\epsilon\right)>J_{a}\left(\epsilon\right)$ for all
$\epsilon\in(0,0.5)$.

Figure 7: Screening coefficients $J_{i}$ as functions of $\kappa$: long-dashed
line indicates screening coefficients $J_{1}=J_{3}$ for partially reflecting
sinks 1 and 3; the solid line indicates the fully absorbing sink 2 with
$J_{2}$.

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

Finally, we also give a numerical example for screening coefficients (see
Subsection IV.7) in the case of three sinks of the same radii ($R_{i}=R$,
$i=\overline{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. 7). Moreover we suppose that sinks possess different
reactivities, namely, sink 2 is fully absorbing ($\kappa_{2}\to\infty$) while
the other two sinks are of the same, but varied reactivity
($\kappa_{1}=\kappa_{3}=\kappa\in\left[0,\infty\right)$). Thus, the
reactivities set is simplified to
$\kappa^{(3)}:=\left\\{\kappa,+\infty,\kappa\right\\}$. Introduce the
notations
$J\left(\kappa\right):=J_{1}\left(\kappa^{(3)}\right)=J_{3}\left(\kappa^{(3)}\right)$
and $J_{2}\left(\kappa\right):=J_{2}\left(\kappa^{(3)}\right)$ and assume that
$a/R=3$.

Fig. 7 shows that in both sinks 1 (3) and sink 2 cases influence of diffusive
interaction slighty decreases with increasing intrinsic rate constant
$\kappa$. Furthermore, screening coefficient $J\left(\kappa\right)$
($J_{2}\left(\kappa\right)$) is monotonically increasing (decreasing)
approaches corresponding horizontal asyymptotes

$\displaystyle\lim_{\kappa\rightarrow+\infty}J\left(\kappa\right)\approx
0.70\,,\quad\lim_{\kappa\rightarrow+\infty}J_{2}\left(\kappa\right)\approx
0.55\,.$ (89)

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

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

$\lim_{\kappa\rightarrow
0+}J\left(\kappa\right)=0\,,\quad\lim_{\kappa\rightarrow
0+}J_{2}\left(\kappa\right)\approx 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
Kantorovich and Krylov (1964); Kantorovich and Akilov (1982). Note in passing
that, although for the microstructure under consideration when
$\kappa\to\infty$ relation (81) gives $M_{3}=1$, the simple-iteration method
also converges.

##  XI 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 (25) and regularity at infinity (27), while
the boundary conditions (26) 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.

###  XI.1 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" Hinsen and Felderhof (1992).

It appears that applying of irreducible Cartesian tensors technique leads to
rather fast convergence of the local concentration expansion (66). Besides,
this approach, contrary to standard one based on the irregular solid spherical
harmonics ${\psi}^{-}_{n}\left(r_{i},\theta_{i},\phi_{i}\right)$ Galanti  _et
al._ (2016); Grebenkov and Traytak (2019); Galanti  _et al._ (2020), 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
$\bm{X}_{n}^{-}\left(\mathbf{r}_{i}\right)$ (A.5) as basis of functions for
solutions of the exterior steady-state diffusion problems in the configuration
manifold $\Omega^{-}$, 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
XII.2). 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 Traytak and Grebenkov (2018). In
other words, for a given $i$-th sink $\Omega_{i}$ it is expedient to choose an
appropriate local curvilinear coordinates
$\\{O_{i};{\bm{\zeta}}\\}\subset{\mathbb{R}}^{3}$, where
${\bm{\zeta}}=(\zeta_{1},\zeta_{2},\zeta_{3})$. Therefore, some coordinate
surface of this coordinate system, say at a fixed index
$m_{i}\in\overline{1,3}$ there exists a real constant $c_{i}$ such that the
surface of $i$-th sink reads:
$\partial\Omega_{i}=\\{{\bm{\zeta}}\in{\mathbb{R}}^{3}:\zeta_{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},\theta_{i},\phi_{i}\\}$ to describe diffusion of $B$-reactants
towards $i$-th spherical sink in domain
$\Omega_{i}^{-}:={\mathbb{R}}^{3}\backslash\overline{\Omega}_{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 $\Omega_{i}^{-}$, $Ox_{3}$ axis and half-plane
$\\{x_{1}\geq 0,x_{2}=0\\}$ should be excluded Fomenko and Mishenko (2009). A
number of other advantages of the ICT technique are described also elsewhere
Coope, Snider, and McCourt (1965); Hess (2015); Snider (2018).

Finally note that from the analytical point of view, obtained in Sec. VIII
solution $u_{N}\left(P\right)$ 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 Eun, Kekenes-Huskey, and McCammon (2013).

###  XI.2 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.
Galanti  _et al._ (2016); Grebenkov and Traytak (2019); Galanti  _et al._
(2020) 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