Competitive effects between stationary chemical reaction centres: a theory based on off-center monopoles.

The subject of this paper is competitive effects between multiple reaction sinks. A theory based on off-center monopoles is developed for the steady-state diffusion equation and for the convection-diffusion equation with a constant flow field. The dipolar approximation for the diffusion equation with two equal reaction centres is compared with the exact solution. The former turns out to be remarkably accurate, even for two touching spheres. Numerical evidence is presented to show that the same holds for larger clusters (with more than two spheres). The theory is extended to the convection-diffusion equation with a constant flow field. As one increases the convective velocity, the competitive effects between the reactive centres gradually become less significant. This is demonstrated for a number of cluster configurations. At high flow velocities, the current methodology breaks down. Fixing this problem will be the subject of future research. The current method is useful as an easy-to-use tool for the calibration of other more complicated models in mass and/or heat transfer.


I. INTRODUCTION
In this paper, chemical reaction systems are studied in which one reacting component (or a catalytic site) is immobilized on a stationary macroscopic object, while a second molecular component moves-by diffusion and/or by some non-stochastic form of motion-towards the static object. This is a good model for many chemical reaction systems in industrial practice (e.g., heterogeneous catalysis in fixed and fluidized beds, combustion of liquid and solid fuels, polymerization reactions, crystal growth) and in biological systems (enzyme-substrate reactions). An interesting special case is the evaporation and subsequent combustion of a spray of finely dispersed liquid fuel droplets, a process in which the transfer of both heat and mass is involved. [1][2][3][4] In systems where multiple stationary reaction centres are in close proximity to each other, the reactive centres compete for the diffusing molecular agent. As a result, the reaction rate per site is lower than in the case of single sites. Much of the work on this subject has been limited to systems in which only diffusive flows are taken into account (e.g., Refs. 3 and 4). Papers on this subject in which non-stochastic flows are included (e.g., Refs. 1 and 2) depend on numerical simulations rather than on analytical or approximate-analytical methods, on account of the complexity of the equations involved. In the current paper, we consider a simple case of constant flow which is amenable to analytical treatment, albeit at the expense of physical realism.
The stationary macroscopic reaction centres are assumed to be spheres for which the spatial coordinates are known a) biello@math.ucdavis.edu b) Author to whom correspondence should be addressed. Electronic mail: rsamsonjvlk@gmail.com and fixed. The theory is developed for an arbitrary but finite number of spheres. The discussion will focus on examples with a relatively small number of spheres.
In the current paper, the solution of the diffusion and the convection-diffusion equation is approximated by introducing fictitious image charges inside the spherical reaction centres, expanding the charge distribution in terms of multipoles and truncating the expansion at the dipolar level. This method has a long history in mathematical physics (see Ref. 17). In recent years, it has again become the subject of intense research activity following the seminal work of Greengard and co-workers on so-called fast multipole moments. This theory has now been developed for conductors, 11,[26][27][28] for dielectric media, [29][30][31][32] for the steady-state 33 and the time-dependent heat equation, 34 for molecular dynamics, 26,35 and in numerous other applications. It is remarkable that the recent flurry of activity in applying image-charge concepts to dielectric media [29][30][31][32] comes more than 130 years after the basic ideas were originally published. 36 In Appendix A, some remarks will be made concerning the relation between our work and the method of images (MOIs) and fast multipoles.

II. DESCRIPTION OF THE MODEL FOR THE LAPLACE EQUATION
We start with the case of purely diffusive flows. In the steady state, the concentration c of the mobile species satisfies the time-independent diffusion equation for ⃗ r ∈ R 3 , with the decay condition at large distances, and an additional internal boundary condition, such as c (S i ) = 0, for i = 1, . . . , N, where {S i |i = 1, . . . , N } are N identical spheres of radius a centred at locations ⃗ r i . Equation (3) represents the case of fully absorbing boundaries. Local angle coordinates are introduced on each sphere so that the position vector on the surface of a sphere is given by wherer i is a unit radial vector of sphere i. A familiar strategy in meeting the surface boundary condition-Eq. (3)-consists of the introduction of fictitious source terms in the Laplace equation, leading to the Poisson equation In the terminology of electrostatics, the q j are charges located at positions ⃗ R j , which are all within the boundaries of one of the spheres. The strengths and the locations of the charges are so chosen that the boundary conditions are satisfied. This approach is related to the well-known method of images (see Refs. 11 and 17, and Appendix A). In general, an infinite number of charges are needed to satisfy boundary condition equation (3) exactly. The crucial finding of the current method is that, in the linear approximation which is outlined below, a consistent (and accurate) solution arises by limiting the charges to one single, off-center charge per sphere (or, equivalently, as we show in Appendix C, to the zeroth and first moments of the charge distribution). The off-center location of the charge gives rise to associated higher-order multipoles. It will be shown that a multipolar expansion that is truncated at the linear (dipolar) level is sufficient to give very accurate results.
The infinite sum in Eq. (5) is replaced by a finite sum with one charge per sphere, where, in general, ⃗ R j ⃗ r j . Using the Green's function, we find the concentration field as where is the [1/distance] Green's function of the Laplace operator. A solution of Poisson equation (5) with source term (6) consists of a set of charges and their coordinates N), such that the concentration field vanishes everywhere on the surface of each sphere, i.e., c (⃗ r i + ar i (θ, ϕ)) = 0, ∀θ, ϕ, i = 1, . . . , N.
In an exact theory, this equation would be satisfied exactly. In an approximate theory, such as the present one, it is a priori clear that this requirement cannot be satisfied exactly but only approximately. An explicit demonstration of this statement is given in Appendix A (see Eq. (A4)), where it is shown that the surface integral of the concentration is not exactly zero for the special case of two touching spheres. Defining d ij = ⃗ d ij being the center-to-center distance between spheres i and j, and s j = ⃗ s j being the off-center shift of charge q j within sphere j, we find

III. SOLUTION STRATEGY
The denominators in Eq. (11) can be expanded using the familiar multipolar expansion where and P l (x) are Legendre polynomials. We shall drop terms of quadrupolar and higher order (terms with l ≥ 2). The term in Eq. (11) under the summation sign (which represents the potential on the ith sphere due to the charges in each of the other spheres) will be expanded in terms of ar i − ⃗ s j /d ij which is clearly <1.
Regarding the second term in Eq. (11) (which represents the potential on the ith sphere due to the single charge within it), it is clear that s i is always <a, therefore we will expand this term in the ratio s i /a. However, we cannot guarantee "a priori" that s i /a is very small. Therefore, we would expect the Taylor series expansion of this term to converge more slowly than that due to the charges in the other spheres.
Substituting the series expansions (12) in (11) to linear order inr i and using the following definitions: we obtain the following equality: Here, Q i is the (non-dimensionalized) charge in sphere i, ⃗ M i is the dipole moment associated with charge Q i and the off-center shift ⃗ s i , δ ij is the Kronecker delta, and α ij and ⃗ A ij are, respectively, a purely symmetric and a purely anti-symmetric matrix. In Eq. (15), Q i and ⃗ M i are the unknowns, while α ij and ⃗ A ij are known variables, which depend on the coordinates of the spheres.
Since Eq. (15) must be approximately true for every point r i (θ, ϕ) on the surface of sphere i, it follows that terms that are independent of and linear inr i must equate to zero separately. In other words, we have a system of equations A solution of Eq. (16) consists of the set Here, U is the vector (1, 1, . . . , 1) ∈ R N , all components of which are equal to unity. Furthermore, the convention of implied summation over repeated indices has been used; superscripts ν = 1, 2, 3 refer to spatial coordinates; subscripts i, j, k = 1, . . . , N refer to sphere labels. L is a symmetric matrix and therefore has real eigenvalues and can be diagonalized by a rotation. Explicitly it can be written as From Eq. (18), we find the charges Q j by inverting the matrix L.
For practical applications (in, for example, mass-or heattransfer problems), this is the most important result so far. In order to find the absorption rate of a system of absorbing objects (e.g., N spheres), we should consider the surface integral of the normal gradient of the concentration of mobile species over all the reactive centres. Applying the Gauss divergence theorem to Poisson equation (5), we find a relation between the integrated normal gradient of the concentration and the average source strength, where Here, S denotes a surface integral over all the spheres andQ is an average source strength (averaged over all N absorption centres). In general,Q will turn out to be <1 (as will be demonstrated further on), reflecting the fact that each reactive site competes for the diffusing species and thereby diminishes the absorbing power of all the other centres in the assembly.

IV. THE LAPLACE EQUATION FOR TWO SPHERES
Consider two spheres with their centers along the x-axis at a center-to-center distance equal to d 12 . Assume that sphere #1 lies to the right of sphere #2, i.e.,d 12 = +ê x (the unit vector along the positive x-axis). Let ∆ ≡ d 12 / (2a), i.e., ∆ ≥ 1, and let β ≡ 1/(2∆), then so that where We also need the antisymmetric matrix A. From Eq. (14) follows The average charge is the same for both spheres (as it must be for reasons of symmetry), the dipole moments are 27) and the locations of the centers of charge are This means that the charges are shifted away from the centers of the spheres in the direction opposite of the other sphere.

A. Comparison of the model with the exact solution for two spheres
The Laplace equation for two spheres is a rare case where an exact solution exists, thanks to the separability of the Laplace operator in so-called bispherical coordinates. 5 This solution has been known for over a century; 6 there is a considerable body of literature on this subject. 7-11 With the same boundary conditions as in Eqs.
(2) and (3), the surface-integral of the normal gradient of the concentration is 12Q where µ is defined by cosh (µ) = ∆. Figures 1 and 2 compare exact result (29) with the result obtained in Sec. III, where β ≡ 1/(2∆). Plotted is the fractional deviation between these two functions (on the y-axis) as a function of ∆ (on the x-axis) for two almost touching spheres ( Figure 1) and for two more widely separated spheres ( Figure 2). For two exactly touching spheres (∆ = 1),Q exact = ln(2) = 0.6931 . . ., whileQ model = 16/23 = 0.6957 . . ., i.e., a fractional difference of 0.0036. It is seen that the dipolar result is a remarkably good approximation of the exact value, up to the smallest distances between the spheres, where the test is at its harshest.
Another way of testing the relation between current model (30) and exact result (29) is the following. One would expect that a Taylor expansion of 1/Q exact in terms of 1/∆ must give the same first three terms as the analogous expansion of 1/Q model , namely, non-vanishing constant, linear, and quartic terms, but no contributions from the quadratic and cubic terms. This is proven in the following. Consider the equalities where Using Eq. (31), it is straightforward (although rather tedious) to show that 1/Q exact has the following expansion in powers of x: However, what we are after is an expansion in powers of 1/∆. From Eq. (32) and cosh (µ) = ∆, it follows that x can be expanded in powers of 1/∆. Substituting that expansion in Eq. (33), we obtain i.e., through order ∆ −4 , the dipolar model and the exact result coincide. The deviation between the two expressions is of order ∆ −6 .

A. Assemblies with full permutational symmetry
In this section, we consider groups of more than two spheres. Symmetric clusters with full permutational symmetry are interesting special cases. In such clusters, every member can be interchanged with any other member without any effect on the absorptive properties of either member. Examples of such symmetric arrangements are the equilateral triangle, the tetrahedron, the planar square, and the cube. A counterexample is the body-centered cube. In the latter case, the sphere in the center of the cube is much more shielded off from the mobile species at close packing conditions than a sphere at one of the vertices of the cube; hence, the central sphere and one of the peripheral spheres cannot be interchanged without affecting their individual absorptive fluxes. In Sec. V C, some body-centered configurations will be considered.
For the assemblies with full permutational symmetry, the interaction matrix L (see Eq. (19)) has interesting special properties. The ith row of L represents the "interactions" of the sphere labelled i with all the other spheres in the assembly. The property of mutual interchangeability implies that rows only differ by a permutation. Consequently, all row-sums of L are equal. Since L is symmetric, also its column-sums are all equal. A matrix with this property is called "semimagical." Such matrices have a very useful property: if L is semi-magical with row-sum Λ, then the inverse of L is also semi-magical with row-sum 1/Λ. 13 Consequently, to find the non-dimensional absorption fluxQ for such permutationally symmetric configurations, we do not even have to invert L explicitly: it is sufficient to calculate its row-sum.
In Appendix B, we list the relation betweenQ and the intersphere distance ∆ for a number of symmetric assemblies that have semi-magical interaction matrices.

B. Comparison of model results with numerical simulations
For the two-sphere system, an exact solution was available with which the model result could be compared. For the larger multi-sphere clusters considered here, this is not the case and we have to rely on other methods to judge the accuracy of the model results.
For this purpose, we have used a commercially available software package: COMSOL multiphysics, a general-purpose finite-element solver. The Laplace equation was solved in a domain consisting of a cluster of N small spheres near the origin of the domain, all of which are enclosed by a much larger sphere. The large sphere is a necessary artefact, because the software cannot deal with an infinitely large domain.
The presence of the large sphere is a complicating factor. Without it the integrated flux for a single small sphere at the origin would simply be 4π. With a large concentric sphere of radius R 2 around the central small sphere (of radius R 1 ), the integrated flux is 4π/(1 − R 1 /R 2 ); i.e., fluxes calculated in this setup are systematically larger than fluxes calculated in an infinite medium by a factor equal to approximately (1 + R 1 /R 2 ). The question which ratio R 2 /R 1 to use in the numerical simulations is a dilemma. On the one hand, one would like to have R 2 as large as possible; however, this choice is in conflict with the requirement to have a sufficiently fine mesh on the surface of the small sphere. By trial and error, we found that a ratio R 2 /R 1 = 100 is a reasonable compromise between these conflicting requirements. At this ratio, the results were better than for smaller ratios; at larger ratios, the software often would not run stably at fine grid sizes.
In a comparative test between the numerical scheme and the dipole method, the most sensible test condition is to pack the spheres in the cluster very densely. Under this condition, the artificial gradient-enlarging tendency in the numerical method due to the presence of the large sphere is as small as possible (R 1 /R 2 being as small as possible), while the dipole approximation is being tested under the most severe conditions.
In Table I, calculated results for the dipolar model (see Appendix B) and numerical results (COMSOL calculated data) are compared for various cluster geometries at ∆ = 1.25, i.e., for a very crowded cluster geometry.
As can be seen from Table I, the fractional deviations between the dipole model and the numerical scheme are of the same order as for the doublet system (<0.01).
Calculations were also carried out at larger ∆-values. The deviations were usually not more than a few percent at most. The numerical fluxes were consistently larger than the dipolar fluxes and the deviations were very often roughly in line with the error estimate based on the "enclosing sphere effect" discussed before. We conclude from this that the discrepancy in those runs is mostly attributable to errors in the numerical method rather than in the dipolar method. Hence, we conclude that the runs at larger ∆-values are not very useful for testing the validity of the dipole method.
Summarizing, this comparison confirms that-just as in the case of two spheres-also for larger assemblies, the dipolar approximation is remarkably accurate: at small intersphere separations, the fractional error is probably <0.01; at relatively large intersphere separations, the fractional error is probably ≪0.01.

C. Assemblies with less symmetry
The validity of the dipolar model is completely independent of the degree of symmetry of the assembly.
Consider the addition of one extra sphere to the center of one of the previously discussed permutationally symmetric clusters. For example, instead of the doublet, consider now three spheres in a row. In this assembly, the terminal spheres are of course not interchangeable with the central sphere. Inserting a central sphere in increasingly larger clusters, e.g., a tetrahedron or a cube, the sphere in the centre will be ever more sterically hindered by its surroundings, leading to lower and lower diffusive fluxes to the central sphere. We cannot use the simplifying feature anymore of semi-magical matrices, but this is not essential for the applicability of the dipolar method.
Body-centered clusters were considered in a recent paper 22 by Eun, Kekenes-Huskey, and McCammon. These authors used numerical methods (finite elements) to solve the Laplace equation and calculate diffusive flux rates to the central sphere. In Table II  a ∆ is the center-to-center separation between nearest-neighbour spheres divided by the sphere diameter (i.e., ∆ ≥ 1). bQ is the non-dimensional flux rate of the mobile species to an arbitrarily chosen sphere in the cluster. aQ is the non-dimensional flux rate of the mobile species to the pertinent sphere; either to one of the peripheral spheres or to the central sphere. b ∆ is the center-to-center separation between nearest-neighbour spheres divided by the sphere diameter (i.e., ∆ ≥ 1).
The tabulated data in Table II were kindly provided to us by Ref. 37.
For all but the closest packings, as the distances within the cluster increase, the fractional deviation between the two methods quickly drops off to less than 0.01, which is quite satisfactory. At or near closest packing however, especially in the larger clusters, the fractional deviation between the two methods is considerable (0.1 to 0.2). It is quite likely that this is related to mesh resolution issues in the numerical method. At close packing conditions, it is very difficult to have a fine enough resolution of the mesh near the surface of the central sphere. Dr. Eun and Dr. Kekenes-Huskey communicated to us that trials to refine their mesh at these close packing conditions led to intractable problems. In view of these problems, the disagreement between the two methods under close packing conditions is not surprising.
Generally speaking, it appears that our approach is both more efficient (in terms of computational speed) and more accurate than the most commonly used numerical methods (e.g., finite-elements).

VI. THE CONVECTION-DIFFUSION EQUATION
In this section, we discuss the situation where there is a non-stochastic flow superimposed on diffusive (i.e., stochastic) motion of the mobile species. Traditionally, this is modelled by the so-called convection-diffusion equation. In this equation, the macroscopic flow is governed by a given velocity field ⃗ u (⃗ r). The relevant Poisson-type equation in the presence of a continuous source field q (⃗ r) is Here, D is the diffusion constant of the mobile species in the flow medium. The same boundary conditions are assumed as before (Eqs. (2) and (3)). If the Green's function for this equation is known, we can-at least in principle-apply the same procedure as before. Although the Green's function is not known for a general flow field, in the special case of a uniform (constant and unidirectional) flow, a Green's function can be constructed, as shown below. In Sec. VII B, we discuss the consequences of the assumption of constancy of the flow field.
With the trial function ϕ (⃗ r) satisfies the modified Helmholtz equation The spatial z-coordinate has been chosen parallel to the flow direction. The Green's function for the modified Helmholtz equation is well known. 14 Collecting factors, we find the following Green's function for Eq. (36) (see also Refs. 15 and 16): This Green's function is not symmetric in its variables, i.e., G (⃗ r,⃗ r ′ ) G (⃗ r ′ ,⃗ r). This lack of symmetry is a consequence of the nonequivalence of up-and downstream. A source at point ⃗ r ′ does not have the same effect at field point ⃗ r as the reverse situation (source at ⃗ r, field point at ⃗ r ′ ). Upstream of a boundary, the flow forces the solute close to the boundary, causing a relatively steep concentration gradient at the boundary; downstream, we have the reverse situation. This lack of symmetry of the Green's function does not necessarily cause any problems. Due to the vanishing of the concentration on the spherical surfaces (Eq. (3)), relation (20) between the integrated surfacegradient of the concentration and the source strength continues to hold in the presence of a convection term (This equality would still hold for a non-zero surface concentration, if the normal component of the velocity at the boundaries were zero.)

A. Dipolar expansion for the convection-diffusion equation in the case of relatively weak flows (i.e., low Péclet numbers)
The treatment of Secs. II and III is now extended to the convection-diffusion equation. This development will turn out to be valid for weakly advective flows only.
As before (see Eq. (11)), we aim to satisfy approximately a complete-absorption boundary condition on the surface of each single one of the N spheres, for ∀i = 1, . . . , N and ∀ (θ, ϕ). The set {⃗ r i } denotes the centers of the N spheres; the set q j ,⃗ s j denotes the charges q j and the locations (off-center shifts ⃗ s j ) within each sphere of those charges. The Green's function G, although not symmetrical in its variables, is translationally invariant, i.e., and r = |⃗ r |.
As before, we separate the term j = i from the sum in Taking ⃗ S = −⃗ s i and ⃗ R = ar i in the first term on the right-hand side of Eq. (43), and ⃗ S = ar i − ⃗ s j and ⃗ R = ⃗ d ij in the second term, and neglecting second-and higher-order derivatives, we obtain Defining the (non-dimensional) charge and the dipole moment in the sphere labelled i, Eq. (45) becomes We now switch to non-dimensional variables noting thatσ ij =d ij , defining furthermore and neglecting all terms of order (r) 2 and higher, we rewrite Eq. (47) as follows: Notice that in the limit that η → 0 (the diffusive limit), ⃗ Γ ij →σ ij /σ ij and ∆ ij → 1/σ ij . As before, we require that terms that do not depend onr i and terms that are linear inr i approximately satisfy Eq. (51) separately. This leads to the following system of equations: (∀i = 1 . . . N), where summation over the repeated index j is implied and where without summation over repeated indices. Equation (52) is the analogon in the convection-diffusion case of Eq. (16) for the purely diffusive case. As before, we have found a closed system of equations relating the charges and associated dipole moments among the set of N spheres. In contrast to the diffusive case however, the coefficient matrices in Eq. (52) do not have the property of being either purely symmetric or purely anti-symmetric; however, that does not stand in the way of the system being solvable, at least in principle.
We are principally interested in calculating the charges Q j since they are needed to calculate the mass flux to the reactive centers, see Eq. (40). The charges are obtained from Eq. (52) as follows: Here, the subscripts i, j, k, l are sphere labels; the superscripts µ, ν are Cartesian indices; the index µ = 3 is along the flow directionû.
In the limit that η → 0, the L-matrix in Eq. (58) is identical to the L-matrix in Eq. (19).
It is useful to rescale the charges such that isolated charges (for the case that all intersphere separations σ ij tend to infinity) are equal to unity. With the current definition of the charges, this is not the case. When σ ij → ∞, ∆ ij → 0, and it follows that where Q 0 denotes the charge strength of an isolated sphere. DefiningQ we now have properly rescaled charges so thatQ i = 1 for isolated spheres. In the development above, expressions were linearised in terms of the parameter η. The question arises: at which value of η is the error still acceptable? We can make a qualitative estimate of the error by plotting the relative error in Eq. (50) in the downflow direction (i.e., forû ·r = −1) as a function of η; see Figure 3. Between η = 0.05 and 0.15, the fractional error increases from 0.001 to 0.01. Clearly, at values of η larger than, say, 0.3, the error becomes unacceptable.

B. Model results (low-η case)
Normalized flux rates were calculated for a number of different cluster configurations (see Table III), 1. a pair of spheres aligned along the flow direction; 2. a pair of spheres at right angles to the flow direction; 3. three spheres arranged in an equilateral triangle in a plane of the flow (e.g., base line of the triangle perpendicular to the flow; altitude line to the base of the triangle along the flow); 4. the same as item 3, with all sides of the triangle perpendicular to the flow; 5. four spheres arranged in a square in a plane of the flow (e.g., one side perpendicular to the flow; a second side along the flow); 6. the same as item 5, all sides of the square perpendicular to the flow; 7. a regular tetrahedron (four spheres) with one of its planes perpendicular to the flow; 8. a cube (eight spheres) with one of its planes perpendicular to the flow.
Two parameters were varied: ∆ and η.
In contrast to our previous calculations in Secs. V A and V B, the centers in a cluster are not necessarily interchangeable. We have to distinguish carefully between centers at up-and downstream locations. In general, downstream sites have lower absorption rates than upstream sites, since their feed stream may be depleted to a serious extent by sites that are further upstream. We shall call this phenomenon the "wake depletion" effect.
From Table III, the following qualitative trends are inferred.

At fixed η,Q is a monotonically increasing function of
∆; in other words: the larger the distances between the reactive sites, the more they behave like isolated centres (no interference effects). There are a few isolated entries in Table III where this is not the case; see the italicized entries with question marks. In these cases, the entry at ∆ = 1.5 is smaller than the one at ∆ = 1. This may be due to the fact that at ∆ = 1.5, the upper hemisphere of the downstream centers is even more starved of solute than at ∆ = 1, without proper compensation from a richer supply to the lower hemisphere. Conversely, these results may also be influenced by numerical errors due to truncation of higher powers of η in the model. 2. At fixed ∆, for upstream centers,Q is a monotonically increasing function of η: as the flow becomes more intense, the upstream centres become less and less sensitive to mutual interference. If the flow is strong enough, it will carry the solute straight up to the surface of an upstream center, regardless of the absorbing pull of laterally displaced neighbouring sites. This rule may or may not apply to downstream centers. More often than not, the reverse tendency holds for downstream centers. This is due to "wake depletion." 3. A cluster that is oriented perpendicularly to the flow suffers less from interference effects than the same cluster in an orientation along the flow direction. This is logical: the flow streamlines have more direct access to the absorbing surfaces in the former orientation. Conversely, in orientations parallel to the flow, the downstream centers are strongly influenced by "wake depletion." In general, within a cluster, the centers that are located further upstream aQ is the non-dimensional flux rate of the mobile species to the pertinent sphere. b ∆ is the center-to-center separation between nearest-neighbour spheres divided by the sphere diameter. c η = au/2D is the Péclet number (apart from a factor 1 /4) for a sphere of radius a in a flow of velocity u; D being the diffusion constant.
have larger flux rates than their downstream neighbors.
(Strangely, this is not the case for the entries in Table III pertaining to the tetrahedron at ∆ = 1 and η = 0.3.) 4. At the same value of ∆ and η, up to moderately large values of η (η 0.3), large clusters generally have a lower Q than small clusters. This stands to reason: in larger clusters, there is generally more steric hindrance than in smaller clusters. At larger values of η however, this is not necessarily the case, as will be shown in Sec. VI D.

C. The high-Péclet case
In Sec. VI A, it was pointed out that the development there is only valid for low Péclet numbers (say, η ≤ 0.3). This makes the applicability of this theory to problems of practical significance somewhat limited. In typical gas-phase problems, η is often in the range 0.1-1; there, the low-Pe theory could be applied, albeit with some reservations as regards the numerical accuracy. However in typical liquid-phase problems, η is usually in the range 10-100 and the low-Pe theory certainly does not apply. Clearly, there is a need for a model that would also cover the high end of the Pe-range, or better yet, a model that is applicable for all values (small and large) of η.
In the high-Pe case, the solution of the convectiondiffusion equation with a constant flow field for one single sphere with Dirichlet boundary conditions is, as far as we know, an unsolved problem. The existing related papers on this subject (e.g., Refs. 23-25) do not precisely cover the problem of our interest, either because the flow obstacle considered is not a sphere 23,24 or because the boundary conditions are not of the Dirichlet type 24 or because the flow field is not constant. 25 In principle, the high-Pe problem lends itself well to a treatment with the well-known methods of perturbation theory, using the inverse Péclet number as a small parameter. However, obtaining a solution that is valid in all regions around the sphere (both in the region of incidence of the flow near the South Pole, in the wake zone extending upwards from the North Pole and in the equatorial region) is not trivial. A full treatment of this subject will be the subject of a future publication.
If we are prepared to take a step back in terms of mathematical rigor and numerical accuracy by considering the monopole approximation (rather than the dipole approximation that has been the focus of our attention thus far), we are in the comfortable position of having a model that applies across the entire Péclet range, as will be demonstrated in Sec. VI D.

D. The monopole approximation for the convection-diffusion equation
Without going into the details of the derivation (which should be rather obvious from the earlier developments), we shall simply state the central result for the monopolar approximation to the convection-diffusion equation. In this section, the monopoles are assumed to be located at the centers of the spheres and are consequently not associated with higher multipoles. In this approximation, the induced charges on the N spheres are obtained as follows: All symbols are as defined in Sec. VI A. Some calculated results are collected in Table IV for values of η ranging from 0 to 100 for three different configurations: 1. configuration #1: a pair of spheres aligned along the flow direction; aQ is the non-dimensional flux rate of the mobile species to the pertinent sphere. b ∆ is the center-to-center separation between nearest-neighbour spheres divided by the sphere diameter. c η = au/2D is the Péclet number.
2. configuration #5: four spheres arranged in a square in a plane of the flow (e.g., one side perpendicular to the flow; a second side along the flow); 3. configuration #8: a cube (eight spheres) with one of its planes perpendicular to the flow.
To a large extent, Table IV (monopole results) exhibits  the same general characteristics as Table III (dipole results). Additional comments regarding Table IV are as follows.
At large values of η (η ≥ 10), the interference effects on the upstream centers are negligible. The flow carrying solute molecules to the absorbers is so strong that interference effects due to neighboring spheres are negligible. This is not the case for downstream centers that are located in the direct wake of upstream centers. Their absorption is strongly affected by the "wake depletion" effect.
Interestingly, at large η, the wake effect is so strong that interference effects from laterally displaced neighboring centers are negligible. This explains why the three configurations considered here have identical absorption rates at high flow rates (η ≥ 10). It is as if the cube (for example) is split up in four independent doublets along the flow direction, with little or no mutual interference between the four laterally displaced doublets.

A. Comparison with the existing state-of-the-art
The solution of the Laplace equation with the help of fictitious charges using multipole expansions is a timehonoured subject in mathematical physics. 11,17,[26][27][28][29][30][31][32][33][34] At first sight, it would seem highly improbable that anything novel can be added to such a well-established subject. Indeed, strictly speaking, we do not believe that we have added anything really new to the subject that was not known already to experts in the field, in the sense that everything that we have presented is somehow implicitly contained in earlier work. As an example of this, in Sec. VI A, we demonstrated that our dipolar result for the doublet coincides to order ∆ −4 with a series expansion of the (known) exact result.
In spite of this, there are a few remarkable observations regarding the current method: 1. in the existing literature, we have not come across any explicit documentation of our results for the doublet and the larger clusters (Eqs. (30) and (B1)-(B4)); 2. For the doublet system, we have shown that the agreement between the current model and the exact solution is surprisingly good (fractional error lower than 0.004). Moreover, we have presented numerical evidence that for larger clusters, the results are of comparable quality as those of the doublet, even at the closest packing conditions.
The fact that such high accuracy can be achieved with a multipolar expansion that stops at the second (dipolar) term has surprised us. This situation is in stark contrast with the method of images which is known to converge exceedingly slowly for closely packed spheres (see, e.g., Ref. 11 and Appendix A below). Although our method has certain similarities with the MOI, there are also some significant differences. The MOI consists of an infinite, iterative procedure, where in each iteration-painstakingly-one charge of a specified strength is put at a specified position inside one of the spheres. In our method, one single charge, representing a superposition of all the charges within the sphere, is located at an off-center position inside one of the spheres. This was shown to lead to surprisingly accurate results for the integrated flux (=normal gradient of the potential). We posit that this procedure is akin to a type of renormalization procedure and that this is the reason for the much faster rate of convergence of our results than those of the conventional MOI (see Appendixes A and C for more details on this point).
As a bonus, the extension of our method from two to N spheres does not pose any special complications, while the conventional MOI leads to very complicated nested iteration schemes in which reflections of reflections of reflections (etc.) have to be taken into account. These complications are simply absent in our method. Here, the result for any number of spheres to any desired multipolar order (monopolar, dipolar, etc.) is found "in one hit" through straightforward matrix methods.

B. Critical comments on the convection-diffusion results
Section VI of this paper deals with the convectiondiffusion model where the flow field is a prescribed constant, unidirectional field everywhere in space. There are two aspects in which this model is lacking in realism.
1. The effect of the boundaries on the flow field is not taken into account. In real life, inevitably some slowing-down of the flow takes place near a solid boundary. Hence, in reality, mass transfer will be more diffusion-and less flow-driven than suggested by our convection-diffusion results. Since the hydrodynamic boundary-layer thickness is almost always larger than or about equal to (in liquids and gases, respectively) the mass-transfer boundary-layer thickness, 18 the current model almost certainly overestimates the effects of forced flow.
On the other hand, if some form of forced convection is at work, a purely diffusive model totally ignores the effects of convection. In this sense, our two models, one with a (constant) flow and the other without any flow at all, could be regarded as an upper and a lower bound for the effect of a realistic flow field (one with flow attenuation near the solid boundaries) in mass or heat transfer problems. Verification of this statement by numerical simulation studies is needed. Clearly, this topic needs further work.

C. The monopolar approximation
Although the method exposed in this paper is based on monopoles only, the off-center positioning of these charges inside the sphere results in associated higher-order multipoles.
In the narrower sense intended in the current section, the monopolar approximation refers to charges that are located at the centre of each sphere. This leads to much simpler results than in the treatment presented in Secs. I-VI of this paper. For ease of reference, we list the most important monopolar results below.
In Sec. VI D, we already discussed the monopolar approximation for the convective-diffusion equation. The results for the diffusive limit are simply obtained from Eq. (57) by taking the limit η → 0. For the doublet (two spheres at a center-to-center distance ∆), the average flux is equal to (compare with the dipolar result Eq. (30)) For the larger clusters studied in Sec. V and Appendix B, the monopolar results are identical to the dipolar results (Eqs. (B1) through (B4)) when terms of order β 4 are ignored (here, β ≡ 1/(2∆)). The fractional error of Eq. (58) versus the exact doublet result (Eq. (29)) is 0.038 for touching spheres (2/3 versus ln(2)) and < or ≪ 0.038 for more widely separated spheres. We see that the dipolar result (fractional error ≤0.0036) is about one order of magnitude better than the monopolar result. For larger clusters, we expect the same relative error as for the doublet.

D. Usefulness of the model
We believe that the primary usefulness of the current model is as a tool for researchers doing numerical simulations in heat and mass transfer. If any competitive effects between nearby boundaries play a role in their computations, researchers might find it handy to have a quick-and-dirty method for estimating this competitive effect on the flux of heat or a chemical species to two or more competing absorbing bodies. The current method provides such a tool. Most practitioners will probably find the current model far easier to use than some of the more sophisticated theoretical papers in the literature. In the current model, the math involved is no more complicated than inverting a simple N-by-N matrix, where N is the number of competing absorbers.

E. Extension of the method to other partial differential equations (PDEs)
In principle, the methods of this paper can be extended to any other linear partial differential equation for which the Green's function is known. Multipole expansions have been used extensively in laminar-flow hydrodynamics (the linearised Navier-Stokes equation). For this case, the Green's function (a tensor) is well known. 19 A considerable amount of work has been done in this field (see, e.g., Refs. 20 and 21).

VIII. CONCLUSIONS
For the steady-state diffusion equation and for the convection-diffusion equation with a constant flow field, a solution method was presented that is based on a multipolar expansion in terms of fictitious source terms. The aim of this work is to quantify the competitive effects among a set of discrete absorbers, which are modelled as spherical bodies. It was found that the inclusion of one single, off-centre monopole per sphere yielded surprisingly accurate results. A closed system of equations for the charges and their locations within each sphere is attained by expanding the boundary conditions on each sphere up to the dipolar term.
For the Laplace equation with two spherical absorbers, the dipolar approximation was compared with the known exact solution in terms of bispherical coordinates. The accuracy of the dipolar model turns out to be remarkably good: the fractional error is about 0.004 for two touching spheres and < or ≪ 0.004 for spheres that are more widely separated.
For the Laplace equation with three or more spherical absorbers, the dipolar solution was compared to numerical simulations. The accuracy of the dipolar model appears to be comparable to the doublet case.
For both the Laplace equation and the convectiondiffusion equation, the use of this method only involves straightforward matrix methods with matrices of dimension N, where N is the number of absorbing macroscopic centres.
The convection-diffusion equation at high values of the Péclet number has not yet been solved. This will be the subject of future work.
The method is a simple tool that can be used by researchers who want to calibrate complex (numerical) models in which competitive effects among multiple absorbing surfaces play a role. One of us (R.S.) gratefully acknowledges stimulating discussions with Professor Dick Bedeaux of NTNU (the Norwegian University of Science and Technology) about an earlier version of this paper.
Both of us furthermore wish to thank Caroline van Dullemen. She introduced us to each other and she graciously rented her Amsterdam apartment (two floors up from René's) to Joseph during his recent sabbatical year in Amsterdam, thereby greatly enhancing our speed of communication.

APPENDIX A: COMPARISON OF THE DIPOLAR APPROXIMATION WITH THE METHOD OF IMAGES
In this appendix, the dipolar approximation is compared with the MOI for two spheres that are just touching (i.e., the center-to-center distance between the spheres is one diameter). Specifically, we shall compare the surface-integral of the concentration (the electric potential in the terminology of electrostatics) of the two methods. If the methods were exact, the surface-integral would be exactly zero, on account of boundary condition (3), but of course this will turn out not to be the case.
The formulation of the MOI to be employed here is the conventional one, as described, e.g., in Ref. 17. At the end of this appendix, we shall add a few words referring to a more up-to-date formulation of the MOI based on the work of Cheng,11 in which the concept of fast multipoles 26-28 is incorporated as well.
In Sec. IV, we found a solution in terms of two source terms Q 1 and Q 2 (given by Eq. (26)), located at ⃗ R 1 and ⃗ R 2 (given by Eq. (28)). These locations do not coincide with the centers of the spheres but they do fall within the boundaries of the spheres. In line with boundary conditions (2) and (3), the concentration field is given by We wish to calculate the integral norm where ⃗ r 1 is on the surface of sphere S 1 and the integral covers the surface of the sphere. Using the well-known equality (see, e.g., Ref . 17) 1 4π it follows that where ρ 2 is the distance from the charge in sphere #2 to the center of sphere #1 and the values of Q and ρ 2 found in Sec. IV for ∆ = 1 were used (Q = 16/23 and ρ 2 /a = 9/4). The MOI solution for two spheres is given as an infinite sum over image charges in the two spheres, is an image charge situated at ⃗ R 1 (k) inside sphere i (i = 1, 2). The Q's and R's are defined through an iterative algorithm, as specified, e.g., in Ref. 11. For two touching spheres, the charges Q (k) are given by the series:  Note that in the MOI, successive charges are of alternating sign. They are situated along the radius pointing towards the touching point, while the two charges in the dipolar case are along the radius pointing away from the touching point. Notice also that the charge in the dipolar case is smaller (by a factor equal to 16/23, in the case of touching spheres) than the primary charge at the origin in the MOI.
It is striking that one single image charge (per sphere) in the dipolar method produces the same result as 207 charges (per sphere) in the MOI in the case of touching spheres. If the spheres were further apart, the MOI would converge faster, and consequently, the discrepancy between the two methods would be less striking than in the case of touching spheres. For example, for ∆ = 2, nine terms suffice to get an equivalent result in the MOI expansion as in the dipolar case.
So far, we have compared our model to the conventional version of the MOI, as developed by Lord Kelvin almost 150 years ago. Recent advances following the work of Greengard et al. [26][27][28] on fast multipoles have improved the rate of convergence of MOI-like schemes by adding image multipoles to the image charges. Cheng 11 has developed the theory for groups of spherical conductors, some of which may be in close proximity to each other. He has found that the inclusion of higher multipoles works best for spheres that are not too close to other spheres. For spheres that have one or several close neighbours, the inclusion of higher multipoles is not very effective. In that case, convergence can be achieved most effectively by including more image charges. For configurations where some but not all of the spheres are closely bunched together, this hybrid approach can lead to significant computational gains.
Just like in the current paper, Cheng has also considered the case of two spheres in close proximity. The comparison of his results with our model shows pros and cons for both methods. For two spheres at a separation of ∆ = 1.2, Cheng needs to include 15 image charges per sphere to keep his error level below 1.E-3. As ∆ → 1, the required number of images in Cheng's calculations increases dramatically. Our theory, with only one image charge per sphere, gives convergent (and remarkably accurate!) results down to ∆ = 1.
To be fair, there are aspects in which Cheng's method unquestionably is superior to ours. The great advantage of his method is that he has a systematic strategy for reducing the numerical error in his computations, by either including more image charges or by increasing the number of included multipoles. In our model in its present form, this is an option we do not have.

APPENDIX B: SPECIFIC FORMULAS FOR SOME ASSEMBLIES WITH FULL PERMUTATIONAL SYMMETRY
In this appendix, we list the relation betweenQ and ∆ for a number of symmetric assemblies that have semi-magical interaction matrices L (see Eq. (19) and Sec. V A). Here,Q is the average source strength (the non-dimensional flux rate), ∆ is the (non-dimensional) distance between nearest neighbours, and β ≡ 1/(2∆).

APPENDIX C: DIGRESSION ON THE COMPRESSION OF ALL THE CHARGES WITHIN A SPHERE TO ONE SINGLE CHARGE PER SPHERE
An exact solution of the Poisson equation, Eq. (5), requires a continuous charge distribution within each sphere. As is standard in the literature, we have approximated the continuous charge density function as an infinite sum of Dirac delta functions.
Next, in going from Eq. (5) to Eq. (6), we replaced the infinite sum by a finite sum (consisting of one single term). Although the use of one single charge per sphere looks like an additional approximation, it will be argued that it is not an additional approximation but an exact result within the context of a linear theory (i.e., a model that only contains terms up to linear order in the Taylor expansion of the Green's function).
Assuming an infinite number of charges per sphere, boundary condition equation (11) can be written as a sum over spheres j and a sum over charges within each sphere k, where q j k and ⃗ s j k are the magnitude and the location of the kth charge in the jth sphere. Now, for each sphere, i = 1 . . . N split the sum over spheres into the term with j = i and group the other terms together, (C2) Then, following the main body of the paper, we now expand the sums up to dipolar (i.e., linear) order. The appropriate expansion parameters are, respectively, s ik /a < 1 for the first term and ar i − ⃗ s j k  d ij < 1 for the second term, yielding From this point onwards, it is clear that we can perform the infinite sums, assuming they converge, which they do. Defining we find which is simply a rearrangement of Eq. (15). The matrices α and ⃗ A were defined in Eq. (14). This shows that the approximation outlined in the main body of the paper can either be seen as the consistent linear approximation describing one off-center charge per sphere or-alternatively-as the consistent linear approximation describing the total charge and the first moment of the charge distribution within each sphere.