Stéphane Balac

Laurent Chupin

Sébastien Martin

Computation of the magnetic potential induced by a collection of spherical particles using series expansions

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

This work is devoted to the study of a way of computing the magnetic field induced by a cluster of metallic particles subjected to the static magnetic field of a Magnetic Resonance Imaging (MRI) device. These magnetic particles locally induce magnetic field inhomogeneities. In some context, magnetic field inhomogeneities are annoying since they perturb the imaging process giving rise to susceptibility artifacts in the image [START_REF] Schenck | The role of magnetic susceptibility in Magnetic Resonance Imaging: MRI magnetic compatibility of the first and second kinds[END_REF]2]. In other contexts, the magnetic field inhomogeneities induced by the metallic particles can be exploited in the imaging process for specific medical diagnostic. For instance, magnetic particles are used for the magnetic labeling of cells and MRI detection. In stem cell therapy, MRI offers the potential of tracking labeled cells in vivo [START_REF] Rogers | Technology insight: in vivo cell tracking by use of MRI[END_REF]. Tumor cell detection can be achieved by labeling the cells with iron oxide-based contrast agents [START_REF] Weis | Labeling of cancer cells with magnetic nanoparticles for MRI[END_REF]. We can also quote the measurement of susceptibility effects due to blood in brain micro-vessel networks and application to the MRI studies with or without contrast agents [START_REF] Chung | Evaluation of the contribution of signals originating from large blood vessels to signals of functionally specific brain areas[END_REF][START_REF] Kim | Cerebral blood volume MRI with intravascular superparamagentic iron oxide nanoparticles[END_REF]. As well, pathological iron deposit in the brain plays a role in neuro-degeneration [START_REF] Zecca | Iron, brain ageing and neurodegenerative disorders[END_REF] and iron has been identified as a potential MRI biomarker for early detection and diagnosis of Alzheimer's disease. Thus, mapping the brain iron content, identifying and quantifying iron deposits using MRI could provide new ways of diagnosing Alzheimer's disease at an early stage [START_REF] Fujiwara | Quantification of iron in the non-human primate brain with diffusion-weighted magnetic resonance imaging[END_REF]. This list of examples is not exhaustive. Over the years, numerical simulation has become an essential tool in the research field of MRI, for a better insight into the physical mechanisms that govern the MR signal, for improvement of simulation sequences and protocols, for safety assessment of medical devices etc. The above mentioned problem is however tricky from a numerical simulation point of view. Indeed, the large amount of particles prevent the use of standard numerical approaches such as the Finite Difference Method [START_REF] Kunz | The finite difference time domain method for electromagnetics[END_REF], the Finite Element Method [START_REF] Jin | The finite element method in electromagnetics[END_REF] or methods based on integral formulations [START_REF] Volakis | Integral equation methods for electromagnetics[END_REF] since they require the meshing of the computational domain, taking into account metallic inclusions, which is prohibitive from a computational point of view. For such a problem, in the existing MRI literature, the calculation of the variations of the static magnetic field is usually done by adding the analytical solutions of the magnetic field due to sources described by simple geometrical shapes: point sources, spheres, (infinite) cylinders. The simplifying assumption is that the field due to each source interacts only with the main magnetic field but does not interact with the fields due to the other sources, see e.g. [START_REF] Zhong | Effects of susceptibility variations on NMR measurements of diffusion[END_REF][START_REF] Boxerman | MR contrast due to intravascular magnetic susceptibility perturbations[END_REF][START_REF] Kiselev | Analytical model of susceptibility-induced MR signal dephasing: Effect of diffusion in a microvascular network[END_REF][START_REF] Kiselev | Effect of magnetic field gradients induced by microvasculature on NMR measurements of molecular self-diffusion in biological tissues[END_REF][START_REF] Pathak | A novel technique for modeling susceptibility-based contrast mechanisms for arbitrary microvascular geometries: The finite perturber method[END_REF]] and references therein. This assumption is referred as the "non-interacting assumption". Sometimes these analytical models are coupled with statistical methods [START_REF] Yablonskiy | Theory of NMR signal behavior in magnetically inhomogeneous tissues: The static dephasing regime[END_REF] or Monte-Carlo methods [START_REF] Boxerman | MR contrast due to intravascular magnetic susceptibility perturbations[END_REF].

The goal of the present work is to study from a mathematical point of view the validity of the "non-interacting assumption". At best, the assumption is justified experimentally in the above mentioned references. Intuition suggests that this assumption is fulfilled when the particles are "far enough" but that it is defective when the particles are "close" to each other. We will quantify in the document how "far" and "close" have to be interpreted.

The features of the particles involved in the above mentioned applications are generally not known precisely. Their shape and size can vary around a mean configuration. However, in order to simplify the mathematical investigations, we assume that all the particles have the same shape and size corresponding to the mean configuration. Namely, we denote by Ω 1 , . . . , Ω N the open sets corresponding to the particles, see Fig. 1. For simplicity and clarity, we assume that the particles are spherical, i.e. the open set Ω j , j ∈ {1, . . . , N }, is a ball centered at position c j ∈ R 3 . We also assume that all the balls have the same radius denoted by ε > 0. The boundary of Ω j is denoted by Σ j . It corresponds to the sphere centered in c j with radius ε. We also denote by Ω c the surrounding area, i.e. Ω c = R 3 \ N j=1 Ω j , and by Σ the boundary of Ω c ,i.e. Σ = N j=1 Σ j . In this study, besides ε, a key parameter is the minimal distance between two particles defined as δ = min i,j∈{1,...,N } i =j

|c i -c j |. (1)
Figure 1. Notations used to describe the cluster of metallic particles.

Moreover, we denote by R > 0 the radius of a ball containing all the particles. Assuming the particles are not in contact to each other and using volume considerations, we have the following constraints:

δ > 2ε and N δ 3 < 8R 3 . (2)
We also assume that the particles are made of the same paramagnetic or weak ferromagnetic material and we denote by µ p > 1 their relative magnetic permeability. As mentioned earlier, the effect of the main static magnetic field H 0 of the MRI device (assumed to be constant in direction and strength in the area where the particles are located) is to generate a secondary magnetic field induced by the metallic particles. As usual in magneto-statics, the mathematical problem can be formulated in term of the scalar magnetic potential ϕ. Our main result states that the error made when approaching the magnetic potential ϕ by the sum of the magnetic potential ϕ c j induced by the particles Ω j , j = 1, . . . , N taken individually, i.e. under the "non-interacting assumption", satisfies, in a convenient norm to be specified later, the following estimate

ϕ - N j=1 ϕ c j W 1 (R 3 ) 4 √ π (µ p -1) 2 µ p + 2 N ε R 3 2 2ε δ 3 B 3 + 3 ln 2R δ R |H 0 |,
where the constant B 3 was found to be B 3 = 142 + 3 ln 2 ≈ 144.1. This result is consistent with intuition: when considering a very large number of particles able to fill the area of interest completely and evenly, the ratio 2ε δ 3 is proportional to the volume fraction of particles. Thus, when the volume fraction is low, that is to say when the particles are spaced at a fairly large distance from each other, computing the magnetic potential by neglecting magnetic field reciprocal interactions between particles gives a quite accurate approximation.

The paper is organized as follows. In Section 2, we introduce and study the magneto-statics problem. In Section 3, we study the accuracy of the "non-interacting assumption" and prove the above mentioned error estimate together with a point-wise error estimate. In Section 4, we study the approach consisting in taking into account pairwise interactions in the computation of the magnetic potential and we prove an error estimate in that case. Finally, Section 5 is devoted to numerical illustrations. [START_REF] Boxerman | MR contrast due to intravascular magnetic susceptibility perturbations[END_REF] 2. The magnetostatic problem 2.1. Formulation for the scalar magnetic potential. The magnetostatic problem described in introduction, expressed for the magnetic field H, reads

   rot H = 0 in R 3 div H = 0 in Ω 1 , . . . , Ω N and Ω c [µH • n] = 0 across Σ 1 , . . . , Σ N (3)
where µ is the piecewise constant function taking the value 1 in Ω c and µ p in Ω j , j = 1, . . . , N . Moreover, we have the following condition at infinity

lim |x|→+∞ x∈R 3 H(x) = H 0 (4)
where H 0 is the applied magnetic field assumed to be constant in strength and direction. By convenience, we denote by | | either the absolute value of a real number or the euclidean norm of a vector in R 3 . In (3) and throughout the paper, we will denote by [u] the jump of the quantity u across one of the boundaries Σ j : [u] = u| Ω c -u| Ω j , j = 1, . . . , N . As well, the quantity n will refer to the outward unit normal to the boundary Σ j .

We reformulate problem (3) in terms of the induced magnetic field H , such that H = H 0 + H . It is natural to assume that H belongs to L 2 (R 3 ) 3 the Hilbert space of square integrable vector functions. Since H is curl free, thanks to Poincaré lemma we can introduce the so-called scalar magnetic potential

ϕ ∈ W 1 0 (R 3 ) such that H = -∇ϕ (5)
where W 1 0 (R 3 ) is the weighted Sobolev space defined by

W 1 0 (R 3 ) = ψ ∈ D (R 3 ) ; ψ 1 + |x| 2 ∈ L 2 (R 3 ) , ∇ψ ∈ L 2 (R 3 ) 3 .
As outlined in [START_REF] Deny | Les espaces du type de Beppo Levi[END_REF] the semi-norm defined by

ψ ∈ W 1 0 (R 3 ) -→ ˆR3 |∇ψ| 2 dx 1 2 is a norm on W 1 0 (R 3 ), denoted by • W 1 (R 3 ) in the sequel. Moreover, all ψ ∈ W 1 0 (R 3 ) satisfies [ψ] = 0 across Σ where Σ denotes any regular surface Σ in R 3 . The magnetic potential ϕ ∈ W 1 0 (R 3 ) satisfies the following problem deduced from (3)    ∆ϕ = 0 in Ω 1 , . . . , Ω N and Ω c µ ∂ϕ ∂n = [µ] H 0 • n across Σ 1 , . . . Σ N (6)
We introduce the following normalization: x i = R x i for i = 1, 2, 3 where x i denotes the dimensionless i-th coordinate. In the dimensionless coordinates system, problem [START_REF] Deny | Les espaces du type de Beppo Levi[END_REF] for the new unknown ϕ defined by

|H 0 |R ϕ( x) = ϕ(R x) for all x ∈ R 3 , reads      ∆ ϕ = 0 in Ω 1 , . . . , Ω N and Ω c µ ∂ ϕ ∂n = [µ] H 0 |H 0 | • n across Σ 1 , . . . , Σ N (7)
where

Ω j = B( c j , ε), c j = c j /R, ε = ε/R and Ω c = R 3 \ N j=1 Ω j .
We also set δ = δ/R. Note that the conditions (2) give rise to the conditions

2 ε < δ < 8 N < 1. (8)
From now on, we will only consider the dimensionless problem [START_REF] Druet | Inégalités de Sobolev optimales et inégalités isopérimétriques sur les variétés[END_REF] and, for convenience, we will drop the tilde symbol introduced to distinguish dimensionless quantities. As well, we do not write down vectors in bold letters anymore. Thus, we have a magnetostatic problem in the form:

     ∆ϕ = 0 in Ω 1 , . . . , Ω N and Ω c µ ∂ϕ ∂n = [µ] h across Σ 1 , . . . , Σ N (9)
for a given h ∈ L ∞ (Σ). The weak formulation of problem [START_REF] Hecht | New development in FreeFem++[END_REF] reads:

find ϕ ∈ W 1 0 (R 3 ) such that for all ψ ∈ W 1 0 (R 3 ) ˆΩc ∇ϕ • ∇ψ dx + µ p N j=1 ˆΩj ∇ϕ • ∇ψ dx = (µ p -1) N j=1 ˛Σj h ψ dσ. (10)
A direct application of Lax-Milgram theorem proves that problem [START_REF] Jin | The finite element method in electromagnetics[END_REF] has a unique solution ϕ ∈ W 1 0 (R 3 ). Moreover, using ϕ as a test function in [START_REF] Jin | The finite element method in electromagnetics[END_REF], we get

ϕ 2 W 1 (R 3 ) (µ p -1) N j=1 h L ∞ (Σ j ) ϕ L 1 (Σ j ) (11)
and from the trace inequality (85) given on Appendix A, we deduce that

ϕ W 1 (R 3 ) A 1 (µ p -1) N ε 3 2 h L ∞ (Σ) , (12)
where h L ∞ (Σ) = sup j∈{1,...,N } h L ∞ (Σ j ) and where A 1 = √ 4π. Applying these results in the special case when h = H 0 |H 0 | • n whose L ∞ -norm equals 1, we have proven the following proposition.

Proposition 1. There exists a unique ϕ ∈ W 1 0 (R 3 ) solution to problem [START_REF] Druet | Inégalités de Sobolev optimales et inégalités isopérimétriques sur les variétés[END_REF]. It satisfies

ϕ W 1 (R 3 ) A 1 (µ p -1) N ε 3 2 , ( 13
)
where

A 1 = √ 4π.
2.2. The one particle case as a reference problem. We consider the following reference problem: find

ϕ ref ∈ W 1 0 (R 3 ) such that    ∆ϕ ref = 0 in B(0, ε) and B(0, ε) µ ∂ϕ ref ∂n = [µ] h across S(0, ε) = ∂B(0, ε) (14)
where

µ = µ p in B(0, ε), µ = 1 in B(0, ε) and h = H 0 |H 0 | • n.
For convenience, we consider a reference frame (0, e 1 , e 2 , e 3 ) such that H 0 defines the direction of e 3 . With this convention, we have h = e 3 • n. Problem ( 14) admits a unique solution ϕ ref ∈ W 1 0 (R 3 ). (Indeed, it corresponds to problem [START_REF] Hecht | New development in FreeFem++[END_REF] in the special case when N = 1 and Ω 1 = B(0, ε).) In this simple geometrical framework, the solution ϕ ref can be computed analytically using the separation of variables method in spherical coordinates. Details are given in appendix C. We have

ϕ ref (x) =        µ p -1 µ p + 2 e 3 • x in B(0, ε) µ p -1 µ p + 2 ε 3 |x| 3 e 3 • x in B(0, ε) (15)
Moreover, the gradient of ϕ ref is given by

∇ϕ ref (x) =        - µ p -1 µ p + 2 e 3 in B(0, ε) - µ p -1 µ p + 2 ε 3 |x| 3 3 e 3 • x x |x| 2 -e 3 in B(0, ε) (16)
From ( 15) and ( 16) we deduce the following estimates. Proposition 2. Let ϕ ref be the solution to problem [START_REF] Kiselev | Analytical model of susceptibility-induced MR signal dephasing: Effect of diffusion in a microvascular network[END_REF]. For all x ∈ B(0, ε), we have

|ϕ ref (x)| µ p -1 µ p + 2 ε 3 |x| 2 and |∇ϕ ref (x)| 2 µ p -1 µ p + 2 ε 3 |x| 3 . ( 17
)
Remark 1. Note that for more general shapes of particles, similar estimates can be obtained using classical results of potential theory, see e.g. [START_REF] Kellogg | Foundations of potential theory[END_REF].

For a spherical particle with radius ε centered at position c j ∈ R 3 , introducing the translation mapping

T c j : x ∈ R 3 -→ x -c j ∈ R 3 , (18)
we readily obtain that the function

ϕ c j = ϕ ref • T c j satisfies    ∆ϕ c j = 0 in Ω j and Ω j µ ∂ϕ c j ∂n = [µ] h across Σ j ( 19
)
The situation is illustrated in Fig. 2. The solution ϕ c j to problem [START_REF] Pathak | A novel technique for modeling susceptibility-based contrast mechanisms for arbitrary microvascular geometries: The finite perturber method[END_REF] is the magnetic potential induced by the metallic particle Ω j submitted to the magnetic field H 0 .

3. Accuracy of the non-interaction assumption 3.1. Approximation under the non-interacting assumption. The solution ϕ to problem [START_REF] Druet | Inégalités de Sobolev optimales et inégalités isopérimétriques sur les variétés[END_REF] can be expressed as

ϕ = N j=1 ϕ c j + ψ. (20)
The first term can be interpreted as the scalar magnetic potential induced by the metallic particles under the assumption that they do not interact together. Since outside a particle, the magnetic field induced by the particle decreases as 1/r 3 where r denotes the distance to the particle center (see Proposition 2), if the particles are at a distance sufficiently large from each other, or if the solid fraction of particles is small, their mutual interaction contribution to the total induced magnetic field can be neglected and it would be reasonable to expect that the approximation ϕ ≈ N j=1 ϕ c j is of good quality. However, if the particles are "not too far" from each other, or if the volume solid fraction of particles is "not too small", this approximation is not anymore accurate. We want to quantify what "not too far" or "not to small" means by studying the behavior of the correction term ψ.

where the data g 0 is defined for all x ∈ Σ i , i ∈ {1, . . . , N }, by

g 0 (x) = - N j=1 j =i ∇ϕ ref (T c j (x)) • n(x). (22)
Proof. By linearity, ψ belongs to W 1 0 (R 3 ) and satisfies the Laplace equation in each domain Ω i , i ∈ {1, . . . , N } and in Ω c , together with the following interface conditions across Σ i , i ∈ {1, . . . , N }:

µ ∂ψ ∂n = µ ∂ϕ ∂n - N j=1 µ ∂ϕ c j ∂n . (23)
Let i ∈ {1, . . . , N }. From ( 7) and ( 19), we have the following condition at the interface

Σ i µ ∂ϕ ∂n = [µ] h and µ ∂ϕ c i ∂n = [µ] h.
We deduce that

µ ∂ψ ∂n = - N j=1 j =i µ ∂ϕ c j ∂n across Σ i . (24)
For any j ∈ {1, . . . , N }, j = i, the potential ϕ c j solution to problem [START_REF] Pathak | A novel technique for modeling susceptibility-based contrast mechanisms for arbitrary microvascular geometries: The finite perturber method[END_REF] is continuous in R 3 and it is regular (regularity C ∞ ) everywhere except on the boundary Σ j . Since for all j = i we have |c i -c j | > 2ε (this condition reflects the fact that the particles Ω i and Ω j with radius ε do not intersect), the normal derivative ∂ϕc j ∂n is continuous across Σ i so that (24) also reads

µ ∂ψ ∂n = -[µ] N j=1 j =i ∂ϕ c j ∂n = -[µ] N j=1 j =i ∇ϕ ref • T c j • n | Σ i . ( 25
)
Proposition 3 is proved.

3.2. W 1 -error estimate under the non-interacting assumption.

Proposition 4. The source term g 0 defined in [START_REF] Schenck | The role of magnetic susceptibility in Magnetic Resonance Imaging: MRI magnetic compatibility of the first and second kinds[END_REF] is such that

g 0 L ∞ (Σ) 2 µ p -1 µ p + 2 2ε δ 3 B 3 + 3 ln 2 δ , (26)
where B 3 = 142 + 3 ln 2 ≈ 144.1.

Proof. The proof relies on Proposition 2. For all x ∈ Σ i , we have

|g 0 (x)| N j=1 j =i ∇ϕ ref (T c j (x)) 2 µ p -1 µ p + 2 N j=1 j =i ε 3 |T c j (x)| 3 . ( 27
)
The quantity T c j (x) is the distance from the center c j of the particle Ω j to the point x ∈ Σ i , see Fig. 3. It is shown in Appendix B, see relation (90), that

N j=1 j =i 1 |T c j (x)| 3 B 3 + 3 ln 2 δ 2 δ 3 . (28)
Combining ( 27) and (28) gives the estimate.

Note that problem [START_REF] Sansone | Orthogonal Functions[END_REF] for ψ has the form of problem [START_REF] Hecht | New development in FreeFem++[END_REF] for ϕ for the boundary data h = g 0 . As a direct consequence of ( 12) and ( 26), we have the following result. Let Ω 1 , . . . , Ω N be N spherical particles in B(0, 1) ⊂ R 3 , centered respectively at c j , j ∈ {1, . . . , N }, with the same radius ε > 0. Let δ be the minimal distance between two particles as defined in (1). Let ϕ c j = ϕ ref • T c j where ϕ ref is given by (15) be the magnetic potential induced by particle Ω j taken in isolation from the others. We have

ϕ - N j=1 ϕ c j W 1 (R 3 ) 2A 1 (µ p -1) 2 µ p + 2 N ε 3 2 2ε δ 3 B 3 + 3 ln 2 δ , (29)
where A 1 = √ 4π and B 3 = 142 + 3 ln 2 ≈ 144.1.

Remark 2. Note that the bound given in relation (29) that corresponds to an estimate of the first order error term in our asymptotic expansion of the magnetic potential ϕ is the product of the term N ε 3 2 of the order of ϕ, see estimate [START_REF] Kiselev | Effect of magnetic field gradients induced by microvasculature on NMR measurements of molecular self-diffusion in biological tissues[END_REF], and a complementary term comparable to 2ε δ 3 1 + ln 2 δ . 3.3. L ∞ -estimate for the error under the non-interacting assumption. In this subsection, we are interested in providing a point-wise estimate for ψ = ϕ -N j=1 ϕ c j . While the estimate of ψ W 1 (R 3 ) was based on Lax-Milgram theorem, our point-wise estimate for ψ relies on integral representation formula for the Laplacian. Namely, the solution ψ to (21) can be expressed, for all x ∈ Ω c , as [5, Thm 2, Chap. XI, p. 120]

ψ(x) = -(µ p -1) N j=1 ˆΣj g 0 (y) G(x, y) dy - ˆΣj ψ(y) ∂ n G(x, y) dy , (30)
where G, the Green kernel of the Laplacian in R 3 , and its normal derivative on Σ j , ∂ n G are given by:

G(x, y) = 1 4π|x -y| , ∂ n G(x, y) = n • (x -y) 4π|x -y| 3 . ( 31
)
We deduce from (30), the following estimate for ψ:

|ψ(x)| (µ p -1) N j=1 g 0 L ∞ (Σ j ) G(x, •) L 1 (Σ j ) + ψ L 4 (Σ j ) ∂ n G(x, •) L 4 3 (Σ j ) , (32)
where the choice of the norms involved in (32) are guided by the regularity of Green's kernel. Indeed, it is known, see [START_REF] Dautray | Analyse mathématique et calcul numérique pour les sciences et les techniques[END_REF]Rem. 4, Chap. XI, p. 127], that G(x, •) is integrable on any boundary Σ j and that ∂ n G(x, •) belongs to L p (Σ j ) for all p < 2. Having a L ∞ (Σ j )-bound on g 0 (see Proposition 3), the way the first integral is bounded is a natural choice. Concerning the second integral, we will use the fact that the Sobolev and trace inequalities implies the following type estimate:

ψ L 4 (Σ j ) ψ W 1 (R 3 )
, and that ψ is already controlled in W 1 (R 3 )norm (see Theorem 1).

Our point-wise estimate for ψ then relies on the estimates provided in the following Lemma.

Lemma 1. Let Ω 1 , . . . , Ω N be N isolated spherical particles in B(0, 1) ⊂ R 3 , with boundaries Σ 1 , . . . , Σ N , respectively centered at c j , j ∈ {1, . . . , N }, with the same radius ε. Let δ be the minimal distance between two particles as defined in (1). For all x ∈ R 3 , we have

N j=1 G(x, •) L 1 (Σ j ) B 1 2 δ 3 ε 2 + M x ε, (33a)
and

N j=1 ∂ n G(x, •) L 4 3 (Σ j ) B 2 2 δ 3 ε 3 2 + M x A 3 ε 3 2 , ( 33b
)
where M x corresponds to the number of particles whose distance at x is less than 2δ (we note that M x 64 and that M x N ) and where B 1 ≈ 8.385, B 2 ≈ 6.611 and A 3 ≈ 0.3439 (all the constants introduced are summarized in table 2 page 36).

Proof. Let us start by proving the estimate for

S 1 = N j=1 G(x, •) L 1 (Σ j ) . For x ∈ Ω c and for j ∈ {1, . . . , N }, the quantity G(x, •) L 1 (Σ j ) is the solution to the following problem (see [5, Thm. 2, Chap. XI, p. 120]):      ∆u = 0 in Ω j and R 3 \ Ω j [u] = 0 across Σ j [∂ n u] = 1 across Σ j (34)
The solution to problem (34) can be computed analytically by using the separation of variable method in spherical coordinates. It is found to be:

u = G(x, •) L 1 (Σ j ) =      ε 2 |x -c j | if x ∈ R 3 \ Ω j ε if x ∈ Ω j (35)
We split the sum S 1 into two parts. The first one, denoted S near 1 in the sequel, corresponds to the indexes j such that the center c j of the particle Ω j is at a distance of x smaller than 2δ and the second one, denoted S far 1 in the sequel, corresponds to the indexes j such that the center c j is at a distance of x strictly larger than 2δ. Accordingly,

S 1 = N i=1 |x-c i |<2δ G(x, •) L 1 (Σ i ) + N j=1 |x-c j | 2δ G(x, •) L 1 (Σ j ) := S near 1 + S far 1 . ( 36
)
The number of terms in the sum S near 1 is bounded independently of N . Since x ∈ Ω c , this number, denoted M x , corresponds to the number of disjointed balls of radius 1 that can be contained in a ball of radius 4. By volume consideration, one can show that M x 64. From (35), we deduce that S near 1 M x ε.

Using relation (35) and a strategy similar to the one used in Appendix B to show relation (90) or (91), we have

S far 1 = N j=1 |x-c j | 2δ ε 2 |x -c j | = 6 ε 2 πδ 3 N j=1 |x-c j | 2δ ˆB(c j , δ 2 ) dy |x -c j | 6 ε 2 πδ 3 N j=1 |x-c j | 2δ ˆB(c j , δ 2 ) dy |x -y| -δ 2 6 ε 2 πδ 3 ˆB(x,2+ δ 2 )\B(x, 3δ 2 ) dy |x -y| -δ 2 24 ε 2 δ 3 ˆ2+ δ 2 3δ 2 r 2 dr r -δ 2 = 24 ε 2 δ 3 ˆ2 δ (s + δ 2 ) 2 ds s B 1 2 δ 3 ε 2 ,
where we have introduced

B 1 = max δ>0 ˆ2 δ 3(s + δ 2 ) 2 ds s = 3 max δ>0 2 + 2δ - 3δ 2 2 + δ 2 4 ln 2 δ ≈ 8.385.
We finally deduce that

S 1 = S near 1 + S far 1 B 1 2 δ 3 ε 2 + M x ε, (37)
where the contribution M x ε is actually only present if x is at a distance less than 2δ from at least one particle.

The estimate of

S 2 = N j=1 ∂ n G(x, •) L 4 3 (Σ j )
is obtained in a very similar way. We split the sum over j into two parts according to the distance between the point x and the spheres Σ j :

S 2 = N i=1 |x-c i | 2δ ∂ n G(x, •) L 4 3 (Σ i ) + N j=1 |x-c j |>2δ ∂ n G(x, •) L 4 3 (Σ j ) := S near 2 + S far 2 . ( 38
)
The number of terms in the first sum S near 2 is, as before, bounded by M x , independently of N . In order to estimate the quantity ∂ n G(x, •)

L 4 3 (Σ i ) when |x -c i | < 2δ, we first introduce the projection x of x on Σ i . We note that (see Appendix E) ∀y ∈ Σ i |∂ n G(x, y)| |∂ n G(x , y)|,
and that the quantity

∂ n G(x , •) L 4 3 (Σ i )
can be explicitly evaluated (see Appendix E for details). We have

∂ n G(x, •) L 4 3 (Σ i ) A 3 ε 3 2 (39)
where

A 3 = ∂ n G((1, 0, 0), •) L 4 3 (S 2 )
≈ 0.3439. (The result is proven in Appendix E for a ball of radius 1 but by dilation one easily reveals the factor ε

2 ). Finally, we obtain

S near 2 M x A 3 ε 3 2 .
Let us now consider the sum S far 2 in (38) over the indexes j such that x is at a distance larger than 2δ from the center c j of the particle Ω j . For x such that |x -c i | > 2δ and for y ∈ Σ j , we have

|x -y| |x -c i | -ε 2δ -ε 2δ |x -c j | > 3 4 |x -c j | so that ∂ n G(x, •) L 4 3 (Σ j ) 1 4π ˛Σj dy |x -y| 8 3 3 4 4 9π ˛Σj dy |x -c j | 8 3 3 4 = 4 9π (4πε 2 )
Summing over j and using the same strategy as in the case of the sum S far 1 , we deduce that

S far 2 B 2 2 δ 3 ε 3 2
where

B 2 = √ 2 3 π 1 4
max δ>0 8 -6δ -δ 2 + 8δ ln( 2 δ ) ≈ 6.611. We finally deduce that

S 2 = S near 2 + S far 2 B 2 2 δ 3 ε 3 2 + M x A 3 ε 3 2 ,
where the contribution

M x A 3 ε 3
2 is only present if x is at a distance smaller than 2δ to any ball.

We are now in position to state the following point-wise estimate for ψ under the same assumptions as for theorem 1.

Theorem 2. Let Ω 1 , . . . , Ω N be N isolated spherical particles in B(0, 1) ⊂ R 3 , centered respectively at c j , j ∈ {1, . . . , N }, with the same radius ε. Let δ be the minimal distance between two particles as defined in (1). Let ϕ c j = ϕ ref • T c j where ϕ ref is given by (15) be the magnetic potential induced by particle Ω j taken in isolation from the others. For all x ∈ Ω c , we have

ϕ(x) - N j=1 ϕ c j (x) 2 (µ p -1) 2 µ p + 2 2ε δ 3 ε B 3 + 3 ln 2 δ B 1 2 δ 3 ε + M x + A 1 A 2 (µ p -1) N ε 2 B 2 2 δ 3 + M x A 3 , (40)
where M x corresponds to the number of particles whose distance at x is less than 2δ (we note that M x 64 and that M x N ). Note that all the constants have been explicitly estimated and are summerized in table 2 on p. 36.

Proof. From (32), ψ = ϕ -N j=1 ϕ c j satisfies for all x ∈ Ω c , |ψ(x)| (µ p -1) N j=1 g 0 L ∞ (Σ j ) G(x, •) L 1 (Σ j ) + ψ L 4 (Σ j ) ∂ n G(x, •) L 4 3 (Σ j ) . ( 41
)
From Lemma 1, we have the estimates

N j=1 G(x, •) L 1 (Σ j ) B 1 2 δ 3 ε 2 + M x ε, N j=1 ∂ n G(x, •) L 4 3 (Σ j ) B 2 2 δ 3 ε 3 2 + M x A 3 ε 3 2 .
It remains to estimate ψ L 4 (Σ j ) . From trace theory and Sobolev embeddings, it is possible to obtain the following inequality (see Appendix A):

ψ L 4 (Σ j ) A 2 ψ W 1 (R 3 ) , (42)
for all j ∈ {1, . . . , N } with

A 2 = 2 3 3 8 √ π
. Since ψ is solution to problem (21) which differs from problem [START_REF] Hecht | New development in FreeFem++[END_REF] only by the boundary data h, ψ satisfies the estimate [START_REF] Kim | Cerebral blood volume MRI with intravascular superparamagentic iron oxide nanoparticles[END_REF] with h replaced by g 0 , i.e. for all j ∈ {1, . . . , N }, we have

ψ L 4 (Σ j ) A 1 A 2 (µ p -1) N ε 3 2 g 0 L ∞ (Σ) . (43)
It follows from (41) that

|ψ(x)| (µ p -1) g 0 ∞ ε B 1 2 δ 3 ε + M x + A 1 A 2 (µ p -1) N ε 2 B 2 2 δ 3 + M x A 3 .
We conclude the proof using the estimate of g 0 ∞ provided in proposition 4.

Remark 3. Both Theorem 1 and Theorem 2 can be interpreted as follows. In the extreme case when we consider only two particles, the error made by adding the contribution to the magnetic potential of the two particles assuming no mutual influence, is roughly speaking proportional to the cube of the ratio of the particles radius to the distance separating the two particles. In the other extreme case when considering a very large number of particles able to fill the area of interest completely and evenly, the ratio ε 3 δ 3 is proportional to the volume fraction of particles. Thus, when the volume fraction is low, that is to say when the particles are spaced at a fairly large distance from each other, computing the magnetic potential by neglecting magnetic field reciprocal interactions between particles gives a quite accurate approximation.

Taking into account pairwise interactions

4.1. Reference problem for the pairwise interactions. When the approximation ϕ ≈ N j=1 ϕ c j is not accurate enough, one can compute the error term ψ in (20) using the same idea as the one introduced in the previous section for the magnetic potential ϕ, that is to say by expanding ψ as

ψ = N j=1 N k=1 k =j ψ c k ,c j + χ (44)
where

ψ c k ,c j , j, k = 1, . . . , N , k = j, are found to satisfiy    ∆ψ c k ,c j = 0 in Ω j and Ω j µ ∂ψ c k ,c j ∂n = [µ] (∇ϕ ref • T c k ) • n across Σ j (45)
It should be noted that each couple of functions (ψ c j ,c k , ψ c k ,c j ) is related to the potential created by the two particles Ω j and Ω k . Namely, the sum ϕ c j +ϕ c k +ψ c j ,c k +ψ c k ,c j corresponds to the potential created by the two particles Ω j and Ω k , without taking into account the other particles.

The solution ψ c k ,c j to problem (45) can be expressed as

ψ c k ,c j = ψ c k -c j ,0 • T c j (46)
where

ψ c k -c j ,0 ∈ W 1 0 (R 3 ) is solution to the following reference problem for u = c k -c j    ∆ψ u ref = 0 in B(0, ε) and B(0, ε) µ ∂ψ u ref ∂n = -[µ] g u across S(0, ε). ( 47
)
where

g u = -(∇ϕ ref • T u ) • n.
Note that when studying the non-interaction assumption in the last section, the reference problem was independent of the ball center: for symmetry reasons, by a translation one could consider a ball centered at origin. When considering pairwise interaction, there is no more symmetry and the reference problem depends on the relative position of the two balls.

The objective now is to obtain estimates for ψ u ref similar to the ones obtained for ϕ ref in proposition 2. The main difficulties are that we do not have an explicit expression for ψ u ref as it was the case for ϕ ref and that there is no longer invariance by rotation of the solution to problem (47) as it was the case for problem [START_REF] Kiselev | Analytical model of susceptibility-induced MR signal dephasing: Effect of diffusion in a microvascular network[END_REF].

4.2.

Explicit solution to the second reference problems. Since ψ u ref satisfies the Laplace equation, it can be expanded in spherical coordinates (r, θ, η) as

ψ u ref (r, θ, η) = +∞ =0 m=- α m r + β m 1 r +1 Y m (θ, η) (48)
where α m and β m are complex numbers to be determined and Y m denotes the Spherical Surface Harmonic function of degree and order m [START_REF] Sansone | Orthogonal Functions[END_REF].

The following lemma specifies the expression of the series expansion of ψ u ref . This expression is purely formal at this stage. Convergence of the series expansion will be proved in proposition 5.

Lemma 2. The Spherical Surface Harmonics series expansion of ψ u

ref reads

ψ u ref (r, θ, η) =              +∞ =0 m=- µ p -1 (µ p + 1) + 1 r ε -1 g m Y m (θ, η) if r < ε +∞ =0 m=- µ p -1 (µ p + 1) + 1 ε 2+ r 1+ g m Y m (θ, η) if r ε
where g m = ´2π 0 ´π 0 g u (θ, η) Y m (θ, η) sin(θ) dθ dη are the coefficients of the Spherical Surface

Harmonics series expansion of the source term g u introduced in (47).

Proof. We consider the formal series expansion (48) of the solution ψ u ref to problem (47). Because the magnetic potential must be bounded for r = 0, it reads in B(0, ε)

ψ u ref (r, θ, η) = +∞ =0 m=- α m r Y m (θ, η). (49)
Moreover, because the magnetic potential must tend to zero at infinity, we have in B(0, ε)

ψ u ref (r, θ, η) = +∞ =0 m=- β m 1 r +1 Y m (θ, η). ( 50
)
The constants α m and β m are determined by the conditions across S(0, ε). The continuity of the solution ψ u ref across the sphere S(0, ε) implies

∀ ∈ N ∀m ∈ {-, . . . , } α m = β m ε 2 +1 . ( 51
)
Let us examine the condition on the normal derivative in (47). From ( 16), we have for all x ∈ S(0, ε)

g u (x) = - µ p -1 µ p + 2 ε 3 |x -u| 5 3 e 3 • (x -u) (x -u) • n -e 3 • n (x -u) • (x -u) . (52)
This function belongs to C ∞ (S(0, ε)) and it can be expanded as a linear combination of Spherical Surface Harmonics as

g u (θ, η) = +∞ =0 m=- g m Y m (θ, η) (53)
where

g m = ˆ2π 0 ˆπ 0 g u (θ, η) Y m (θ, η) sin(θ) dθ dη. (54)
Thus, the interface condition in (45) is found to give rise to the relations

∀ ∈ N ∀m ∈ {-, . . . , } µ p α m ε -1 + ( + 1) β m ε +2 = (µ p -1) g m .
From (51), it follows that

∀ ∈ N ∀m ∈ {-, . . . , } β m = µ p -1 (µ p + 1) + 1 ε +2 g m . ( 55
)
According to Lemma 2, in order to estimate the behavior of the solution ψ u ref to the reference problem (47), we have to control the coefficients g m . This is the subject of the following lemma where the notation ∆ s refers to the Laplace-Beltrami operator on the sphere

∆ s = 1 sin(θ) ∂ θ sin(θ)∂ θ + 1 sin 2 (θ) ∂ 2 η .
For all integer p 1, we set

∆ p s = ∆ s • • • • • ∆ s (p times)
and by convention ∆ 0 s is the identity. For convenience, we also introduce the gradient operator on the sphere :

∇ s =    ∂ ∂θ 1 sin(θ) ∂ ∂η    Lemma 3. Coefficients g m for ∈ N, m ∈ {-, .
. . , } in the spherical harmonic expansion of g u defined in (53)-(54) satisfy g 0 0 = 0 and for any choice of p ∈ N

∀ ∈ N \ {0} ∀m ∈ {-, . . . , } |g m | √ 4π ( ( + 1)) p ∆ p s g u ∞ ,
with

g u ∞ 4 µ p -1 µ p + 2 ε 3 (|u| -ε) 3 , (56a)
∆ s g u ∞ µ p -1 µ p + 2 ε 3 (|u| -ε) 5 450ε 2 + 139ε|u| + 8|u| 2 . ( 56b
)
Proof. Let us first show that g 0 0 = 0. From (54), we have

g 0 0 = - 1 √ 2π ˛S(0,ε) ∂ϕ ref ∂n (x -u) dσ x = - 1 √ 2π ˛u+S(0,ε) ∂ϕ ref ∂n (y) dσ y .
It follows from Stokes' formula that

g 0 0 = - 1 2 √ π ˆu+B(0,ε) ∆ϕ ref (y) dy = 0.
Let us now consider the estimates of g m for 1. To highlight the decrease in coefficients with respect to , the idea is to use the fact that the elements Y m are eigenfunctions to the laplacian. More precisely, we have

∆ s Y m = -( + 1) Y m . It follows from (54) that g m = - 1 ( + 1) ˆ2π 0 ˆπ 0 g u (θ, η) ∆ s Y m (θ, η) sin(θ) dθ dη.
Since g u ∈ C ∞ (S(0, ε)), using Green's formula, we have

g m = - 1 ( + 1) ˆ2π 0 ˆπ 0 ∆ s g u (θ, η) Y m (θ, η) sin(θ) dθ dη.
Then, repeated use of Green's formula p -1 more times yields

g m = (-1) p ( ( + 1)) p ˆ2π 0 ˆπ 0 ∆ p s g u (θ, η) Y m (θ, η) sin(θ) dθ dη.
It follows that

|g m | 1 ( ( + 1)) p ˆ2π 0 ˆπ 0 ∆ p s g u (θ, η) Y m (θ, η) sin(θ) dθ dη 1 ( ( + 1)) p ∆ p s g u ∞ ˆ2π 0 ˆπ 0 Y m (θ, η) sin(θ) dθ dη √ 4π ( ( + 1)) p ∆ p s g u ∞ .
This complete the first part of the proof. It remains to show the two estimates (56). From (52), we have for all x ∈ S(0, ε)

|g u (x)| 4 µ p -1 µ p + 2 ε 3 |x -u| 3 .
The estimate (56a) follows from the inequality: |x -u| |u| -|x| = |u| -ε.

In order to obtain the estimate for ∆ s g u ∞ , we proceed as follows. We have

g u (x) = - µ p -1 µ p + 2 3 h u ( x) | x -u| 5 - e 3 • n | x -u| 3 , (57)
where we have set x = x/ε, u = u/ε and

h u ( x) = e 3 • ( x -u) ( x -u) • n . For p ∈ N, p 1, we set ∀ x ∈ S 2 and ∀ u / ∈ S 2 G p ( x, u) = 1 | x -u| p . A simple, yet cumbersome, direct calculation shows that ∆ s G p ( x, u) = p(p + 2)|J s ( x) u| 2 G p+4 ( x, u) -2p 1 + ( x • ( u -x)) G p+2 ( x, u)
where J s ( x) denotes the following 2 × 3 matrix

J s ( x) = cos(θ) cos(η) cos(θ) sin(η) -sin(η) -sin(η) cos(η) 0 = ∇ s x.
Moreover, we have

∇ s G p ( x, u) = p J s ( x) u G p+2 ( x, u). Note that J s ( x) u = J s ( x) ( u -x) and |J s ( x) u| √ 2 | x -u|. Therefore, |∆ s G p ( x, u)| 2p(p + 3) G p+2 ( x, u) + 2p G p+1 ( x, u) (58)
|∇ s G p ( x, u)| √ 2p G p+1 ( x, u). ( 59
)
Another direct calculation shows that

∆ s h u ( x) = 2 1 -n • ( x -u) ( x -u) • e 3 -2 x • e 3 ( x -u) • n + 2 ∂ θ x • e 3 ∂ θ x • ( x -u) and ∇ s h u ( x) = ( x -u) • n J s ( x) e 3 -( x -u) • e 3 J s ( x) u . Because | x| = |∂ θ x| = |∂ η x| = 1, we have the estimates |h u ( x)| | x -u| 2 (60) |∆ s h u ( x)| 2 | x -u| 2 + 6 | x -u| (61)
|∇ s h u ( x)| √ 2 | x -u| 2 + | x -u|. (62)
Now, using the differential calculus formula

∆ s (f × g) = f ∆ s g + g ∆ s f + 2 ∇ s f • ∇ s g we deduce that ∆ s g u (x) = - µ p -1 µ p + 2 3h u ( x) ∆ s G 5 ( x, u) + 3G 5 ( x, u) ∆ s h u ( x) + 6∇ s G 5 ( x, u) • ∇ s h u ( x) -e 3 • n ∆ s G 3 ( x, u) -G 3 ( x, u) ∆ s e 3 • n -2∇ s G 3 ( x, u) • ∇ s e 3 • n .
To get the estimate for ∆ s g u (x) we proceed by estimating each term involved in the preceding expression:

|∆ s g u (x)| µ p -1 µ p + 2 8G 3 ( x, u) + 6(19 + √ 2)G 4 ( x, u) + 6(46 + 5 √ 2)G 5 ( x, u) µ p -1 µ p + 2 1 | x -u| 5 8| x -u| 2 + 123| x -u| + 219 = µ p -1 µ p + 2 ε 3 |x -u| 5 8|x -u| 2 + 123ε|x -u| + 219ε 2
From the inequalities |x -u| |x| + |u| and |x -u| |u| -|x| = |u| -ε, we finally obtain

|∆ s g u (x)| µ p -1 µ p + 2 ε 3 (|u| -ε) 5 450ε 2 + 139ε|u| + 8|u| 2 . ( 63
)
We are now in position to prove convergence of the series expansion for ψ u ref introduced in lemma 2.

Proposition 5. The series (50) where β m , ∈ N, m ∈ {-, . . . , }, is defined in (55) converges uniformly on B(0, ε). The series (49) where α m , ∈ N, m ∈ {-, . . . , }, is defined in (51) converges uniformly on B(0, ε).

Proof. Let a (r, θ, η) = m=-α m r Y m (θ, η) be the coefficient of rank of the series (49). From (51) and (55), we have

a (r, θ, η) = (µ p -1)ε (µ p + 1) + 1 r ε m=- g m Y m (θ, η).
It follows from lemma 3 (considered with p = 1) that

|a (r, θ, η)| √ 4π ( + 1) ∆ s g u ∞ (µ p -1)ε (µ p + 1) + 1 r ε m=- |Y m (θ, η)| . (64)
Using the following estimate for Spherical Surface Harmonics (see Appendix D for a proof)

m=- |Y m (θ, η)| 3 √ 3 √ 4π 3 2 (65)
we conclude that ∈N a (r, θ, η) is uniformly convergent on B(0, ε). Similarly, let b (r, θ, η) = m=-β m 1 r +1 Y m (θ, η) be the coefficient of rank of the series (50). From (55), we have

b (r, θ, η) = (µ p -1)ε (µ p + 1) + 1 ε r +1 m=- g m Y m (θ, η). ( 66
)
From lemma 3 and from (65), we deduce that

|b (r, θ, η)| (µ p -1)ε (µ p + 1) + 1 ∆ s g u ∞,S(0,ε) 3 √ 3 3 2 ( + 1) ε r +1 (67)
and we conclude that the function series ∈N b (r, θ, η) converges uniformly on B(0, ε).

Finally, we are in position to prove the following result for ψ u ref which is the analogue of proposition 2 for ϕ ref .

Proposition 6. Let ψ u ref be the solution to problem (47). For all x ∈ B(0, ε), we have

|ψ u ref (x)| 12 √ 3 λ p ε 3 (|u| -ε) 3 ε 2 |x| Li -1 2 ε |x| |∇ψ u ref (x)| A 4 λ p ε 3 (|u| -ε) 3 ε 2 |x| 2 Li -3 2 ε |x| (68)
where λ p := (µ p -1) 2 (µ p + 1)(µ p + 2) ∈ [0, 1] and Li s (z) refers to the Polylogarithm function of order s (also called Jonquière's function). Moreover, the constant A 4 is such that A 4 ≈ 59.93.

Proof. The estimate of ψ u ref (x) is deduced from the Spherical Surface Harmonic series expansion given in Lemma 2 for r ε combined with the estimate of the coefficients g m given in lemma 3 for p = 0 and with the estimate (65) for the Spherical Surface Harmonics.

Let us consider the estimate of ∇ψ u ref (x). In the spherical coordinate basis, we have

|∇ψ u ref (x)| 2 = |∂ r ψ u ref (r, θ, η)| 2 + | 1 r ∂ θ ψ u ref (r, θ, η)| 2 + | 1 r sin(θ) ∂ η ψ u ref (r, θ, η)| 2 . ( 69
) Since ψ u ref (r, θ, η) = +∞ =1 m=-β m 1 r +1 Y m (θ, η) in B(0, ε), we deduce that ∂ r ψ u ref (r, θ, η) = - +∞ =1 m=- β m + 1 r +2 Y m (θ, η)
and an estimate for ∂ r ψ u ref can be obtained in exactly the same way as for ψ u ref .

Namely, using -estimate (55) for coefficients β m -estimate of the coefficients g m given in lemma 3 considered with p = 0

|g m | 4 √ 4π µ p -1 µ p + 2 ε 3 (|u| -ε) 3 (70)
-estimate (65) on the Spherical Surface Harmonic functions Y m we find that

|∂ r ψ u ref (r, θ, η)| 24 √ 3 λ p ε 3 (|u| -ε) 3 ε 2 r 2 Li -3 2 ε r . (71)
Let us consider the second component of ∇ψ u ref given by

1 r ∂ θ ψ u ref (r, θ, η) = +∞ =1 m=- β m 1 r +2 ∂ θ Y m (θ, η).
Here again, using estimate (55) for coefficients β m combined with the estimate g m given in lemma 3 considered with p = 0 and the estimate of m=-|∂ θ Y m (θ, η)| given by (117) in Appendix D, we get

1 r ∂ θ ψ u ref (r, θ, η) 4 √ 6(1 + π) λ p ε 3 (|u| -ε) 3 ε 2 r 2 Li -3 2 ε r . ( 72
)
The third component of ∇ψ u ref is given by

1 r sin(θ) ∂ η ψ u ref (r, θ, η) = +∞ =1 m=- β m 1 r +2 1 sin(θ) ∂ η Y m (θ, η).
Proceeding as before using now the estimate of m=-

1 sin(θ) ∂ η Y m (θ, η)
given by (118) in Appendix D, we get

1 r sin(θ) ∂ η ψ u ref (r, θ, η) 6 √ 6 λ p ε 3 (|u| -ε) 3 ε 2 r 2 Li -3 2 ε r . (73)
Finally, combining (71), ( 72) and (73) together with (69), we conclude that

|∇ψ u ref (x)| A 4 λ p ε 3 (|u| -ε) 3 ε 2 |x| 2 Li -3 2 ε |x| , (74)
where

A 4 = (24 √ 3) 2 + (4 √ 6(1 + π)) 2 + (6 √ 6) 2 = 2 6(85 + 8π + 4π 2 ) ≈ 59.93.
4.3. W 1 -error estimate when taking into account pairwise interactions. When the approximation ϕ ≈ N j=1 ϕ c j is not accurate enough, a more accurate approximate solution can be obtained by using the two terms expansion

ϕ = N j=1 ϕ c j + N j=1 N k=1 k =j ψ c k ,c j + χ ( 75
)
where χ denotes the approximation error.

Proposition 7. The remainder χ in (75) satisfies

   ∆χ = 0 in Ω 1 , . . . , Ω N and Ω c µ ∂χ ∂n = [µ] g 1 across Σ 1 , . . . , Σ N (76)
where the data g 1 is defined by

∀x ∈ Σ i g 1 (x) = - N j=1 j =i N k=1 k =j ∇ψ c k -c j ,0 (T c j (x)) • n(x).
Moreover, the source term g 1 satisfies

g 1 L ∞ A 4 λ p 2ε δ 5 B 3 + 3 ln 2 δ B 3 Li -3 2 2ε δ + 3 2ε δ ln 2 δ . ( 77
)
Proof. Since the functions ϕ, ϕ c j and ψ c k ,c j for j, k ∈ {1, . . . , N }, k = j, all satisfy the Laplace equation in each domain Ω i , i = 1, . . . , N and in Ω c , it follows from (75) by linearity that χ also satisfies the Laplace equation. Moreover, for all i ∈ {1, . . . , N } we have

µ ∂χ ∂n = µ ∂ψ ∂n - N j=1 N k=1 k =j µ ∂ψ c k ,c j ∂n across Σ i . (78)
From ( 21), ( 22) and (45), we can express the interface condition (78) across Σ i as

µ ∂χ ∂n = [µ] g 0 - N k=1 k =i µ ∂ψ c k ,c i ∂n = 0 - N j=1 j =i N k=1 k =j µ ∂ψ c k ,c j ∂n . ( 79
)
The potential ψ c k ,c j solution to problem (45) is continuous in R 3 and has regularity C ∞ everywhere except on the boundary Σ j . Therefore, as soon as |c i -c j | > 2ε (that corresponds to the fact that the particles Ω i and Ω j of radius ε do not intersect) the normal derivative ∂ψc k ,c j ∂n is continuous across Σ i and (79) reads

µ ∂χ ∂n = -[µ] N j=1 j =i N k=1 k =j ∂ψ c k ,c j ∂n = -[µ] N j=1 j =i N k=1 k =j ∇ψ c k -c j ,0 • T c j • n. (80)
The first result stated in the lemma is proven. In order to obtain the estimate on the source term g 1 , we start with

∀i ∈ {1, . . . , N } ∀x ∈ Σ i |g 1 (x)| N j=1 j =i N k=1 k =j |∇ψ c k -c j ref (x -c j )|.
We deduce from Proposition 6 that

|g 1 (x)| A 4 λ p N j=1 j =i N k=1 k =j ε 5 (|c k -c j | -ε) 3 |c j -x| 2 Li -3 2 ε |c j -x| .
From the definition of the Polylogarithm function and the inequality

|c k -c j | -ε 1 2 |c k -c j |, it follows that |g 1 (x)| A 4 λ p +∞ =1 3 2 ε +2 N j=1 j =i 1 |c j -x| +2 N k=1 k =j ε 3 |c k -c j | 3 .
Using relation (90) of Appendix B, we can bound the last sum to obtain

|g 1 (x)| A 4 λ p 2ε δ 3 B 3 + 3 ln 2 δ +∞ =1 3 2 ε +2 N j=1 j =i 1 |c j -x| +2 .
Finally, using relations (90) and (91) of Appendix B, we can bound the sums

N j=1 j =i 1 |c j -x| +2
and for all 1 to conclude that

|g 1 (x)| A 4 λ p 2ε δ 5 B 3 + 3 ln 2 δ B 3 Li -3 2 2ε δ + 3 2ε δ ln 2 δ .
Following the same steps as in the proof of Theorem 1, we deduce the following result.

Theorem 3. Let Ω 1 , . . . , Ω N be N isolated spherical particles in B(0, 1) ⊂ R 3 , centered respectively at c j , j ∈ {1, . . . , N } with the same radius ε. For all j, k ∈ {1, . . . , N }, i = k, let ϕ c j ∈ W 1 0 (R 3 ) be the solution to problem [START_REF] Pathak | A novel technique for modeling susceptibility-based contrast mechanisms for arbitrary microvascular geometries: The finite perturber method[END_REF] and ψ c k ,c j ∈ W 1 0 (R 3 ) be the solution to problem (45). We have

ϕ - N j=1 ϕ c j + N j=1 N k=1 k =j ψ c k ,c j W 1 (R 3 ) A 1 A 4 (µ p -1) λ p N ε 3 2 2ε δ 5 B 3 + 3 ln 2 δ × B 3 Li -3 2 2ε δ + 3 2ε δ ln 2 δ ( 81
)
where δ is the minimal distance between two particles as defined in (1). Note that all the constants have been explicitly estimated and are summerized in table 2 on p. 36.

ϕ c j + N j=1 N k=1 k =j ψ c k ,c j W 1 (R 3 ) (µ p -1) λ p N ε 3 2 2ε δ 6 1 + ln 2 δ 2
and this new bound is the product of a term N ε 3 2 of the order of ϕ, and the square of the first order error 2ε δ 3 1 + ln 2 δ .

Numerical investigations

Based on the mathematical investigations conducted in the previous sections, we have developed a computer program under Matlab to compute the magnetic potential ϕ solution to problem [START_REF] Deny | Les espaces du type de Beppo Levi[END_REF]. The computer program offers the possibility to compute the magnetic potential with either one term or two terms in the series expansion corresponding respectively to the results stated in theorem 1 and theorem 3. Computation of the first term in the expansion given by N j=1 ϕ c j is straightforward since ϕ c j = ϕ ref • T c j and the expression of ϕ ref is given by [START_REF] Kunis | Fast spherical fourier algorithms[END_REF]. Computation of the second term

N j=1 N k=1 k =j ψ c k ,c j , where ψ c k ,c j = ψ c k -c j ref • T c j ,
that takes into account pairwise interaction relies on the series expansion given by Lemma 2. Namely, for all x = (r, θ, η) ∈ R 3

ψ c k -cj ref (r, θ, η) =              +∞ =1 r ε m=- γ m Y m (θ, η) if |x -c j | ε +∞ =1 ε r +1 m=- γ m Y m (θ, η) if |x -c j | > ε (82)
where we have set

γ m = (µ p -1) ε (µ p + 1) + 1 ˆ2π 0 ˆπ 0 g(θ, η) Y m (θ, η) sin(θ) dθ dη. ( 83
)
Of course, in a numerical approach, the series expansions (82) have to be truncated to an integer L and the efficiency of the method depends on the convergence to zero of the series coefficients γ m as increases to infinity.

5.2.

Illustration of the mathematical results. We start with an illustration of the behavior of the bound in the estimate provided in theorem 2. We consider two spherical particles (N = 2) with radius ε = 1 mm aligned in the e 1 direction. For two particles, we have

E := max x∈R 3 ϕ(x) - 2 j=1 ϕ c j (x) = max x∈R 3 2 j=1 2 k=1 k =j ψ c k ,c j (84)
since with two particles, there is nothing else than pairwise interactions. We have computed E considering two particles at various distance δ form each other for µ p = 10 and H 0 = 1 µ 0 1 0 0 H.m -1 . We have depicted in Fig. 6 the quantity E as a function of the distance δ between the two spherical particles in log-scale. One can observe that E varies as 1/δ 3 as predicted by theorem 2.

We then consider two spherical particles (N = 2) with radius ε aligned in the e 1 direction with centers at a distance δ = 11 mm. We have computed E considering various radius ε for µ p = 10 and H 0 = 1 µ 0 1 0 0 H.m -1 . We have depicted in Fig. 7 the quantity E as a function of the radius ε of the two spherical particles in log-scale. One can observe that E varies as ε 4 as predicted by theorem 2. By linear least-square interpolation, the vertical intercept is found to be 1.053865 × 10 3 . In the estimate provided in theorem 2, for δ = 1, the first term in the series expansion in ε is 16

M x N = 2.
The bound provided in theorem 2 is therefore not a sharp bound in this particular case.

Numerical experiments.

As a first numerical experiment, we consider the case of two spherical particles with magnetic permeability µ p = 10 and radius ε = 1 mm subjected to a magnetic field aligned with the ball centers with intensity H 0 = 1 µ H.m -1 . The two ball centers are separated by δ = 3 mm, i.e. the two spheres are at a distance of one radius. The plane of reference where the magnetic potential is computed is a plane passing through the ball centers with size 9ε along the y-axis and 4ε along the x-axis. We have depicted in Fig. 8 the potential ϕ c 1 + ϕ c 2 corresponding to the first term in our series expansion, i.e. obtained by neglecting interaction between the balls, the potential ψ c 1 ,c 2 +ψ c 2 ,c 2 corresponding to the second term in the series expansion, i.e. corresponding to the pairwise interactions between the two balls and the second order series approximation of the magnetic potential In this first test case, one can see that the second term in the series expansion that takes into account pairwise interactions between the two balls amounts to around 7 % of the total magnetic potential. In order to validate our computer program based on the series expansion presented in the paper, we have also compared the second order series approximation of the magnetic potential ϕ ≈ ϕ c 1 + ϕ c 2 + ψ c 1 ,c 2 + ψ c 2 ,c 2 to the magnetic potential ϕ computed by the Finite Element Method (FEM) using the free software FreeFem++ [START_REF] Hecht | New development in FreeFem++[END_REF]. Note that the solution provided by the FEM is to be considered as a reference solution. It is not the exact solution because of the discretization error (the magnetic potential is approached by Lagrange P 2 Finite Element on a fine triangulation mesh of the computational domain) and the need to bound the computational domain by introducing an artificial boundary and a (approximate) boundary condition that takes account of the behavior of the magnetic potential at infinity. The relative maximal error over the computational grid on the plane of reference between the FEM solution and the solution computed by the two terms series expansion was found to be 0.67% whereas the relative quadratic error was 0.63%. For the sake of completeness, we should add that comparing the FEM solution to the solution ϕ c 1 + ϕ c 2 obtained under the non-interaction assumption we obtain a relative maximal error of 5.78% and a relative quadratic error of 5.18%. Thus, in this test example, we can conclude that taking into account pairwise interactions significantly improves the accuracy.

ϕ ≈ ϕ c 1 + ϕ c 2 + ψ c 1 ,c 2 + ψ c 2 ,
To go further in the numerical investigations, we have computed the relative maximal error and the relative quadratic error for various values of the parameter δ/ε when only one term in the potential series expansion is taken into account (this corresponds to the "non-interacting assumption" case) and when the first two terms are taken into account (this corresponds to the "pairwise interactions" case). One can see on the results given in Table 1 that when the two ball centers are at a distance lower than 5ε, the error made on the computation of the magnetic potential assuming that the particles have no magnetic interaction lead to a significant error. On the contrary, taking into account pairwise interactions provide good results with an error about one percent. One can also observe that when the particles become very close to each other, the error obtained with two terms in the potential series expansion tends to increase. In our test example, it is around 3% when the particles are at a distance of a quarter of their radius. This indicates that when the particles are very close to each other, one should compute a third term in the potential series expansion. 1. Maximal Error and quadratic error obtained when using one and two terms in the magnetic potential series expansion as a function of δ ε , i.e. as a function of the distance between the two ball centers. Note that for δ ε = 2, the balls are touching each other. As a second test example, we have considered a cube of 4 3 = 64 balls with radius ε = 1 mm regularly spaced with a distance between their centers δ = 3 mm. The balls are subjected to a magnetic field aligned with x-direction with intensity H 0 = 1 µ 0 H.m -1 . The plane of reference where the magnetic potential is computed is a plane passing through the center of the cube with basis vectors corresponding to the x and y directions. The magnetic potential is computed on 100 2 equidistant points on a square grid with size 15ε. Computations lasted 59 872 s. (approx. 16h30) to compute the 64 × 63 = 4032 particles pairwise interactions. In this test example, the second term of the series expansion is about 3.5 % of the first term. To temper the important simulation time observed, several remarks can be made. The first one is that the main computational cost is the computation, for each pairwise interaction, of the series coefficients γ m defined in (83) for all m = -, . . . , and for all = 1, . . . , L where L is the number of terms taken into account in the truncated series (we took L = 16). We have used a basic quadrature approach to evaluate γ m but actually the spherical harmonics integrals involved in the definition (83) of γ m could be computed by fast spherical Fourier algorithms [START_REF] Kunis | Fast spherical fourier algorithms[END_REF]. We can add that for the sake of graphical representation, the magnetic potential was computed on a bi-dimensional grid but computation on a three-dimensional grid would not significantly increase the overall computation time. The second remark is that simulation was achieved on a desktop computer, not on a server dedicated to scientific computing. In such a computer, computation time could be drastically reduced by parallel programming. Indeed, each particle pairwise interaction can be computed independently from each other.

Conclusion

In the first part of the paper, we have investigated the validity of the approach consisting in neglecting magnetic interactions between particles in the numerical computation of the magnetic field inhomogeneities induced by a cluster of metallic particles subjected to a uniform magnetic field. Such a situation is encountered in various contexts in Magnetic Resonance Imaging. We have obtained bounds for the approximation error in terms of the three geometrical key parameters that were found to be the number of particles N , the radius ε of the particles (assumed to have a spherical shape for simplicity), and the distance δ between the two nearest particles. Two other physical parameters, the magnetic permeability µ p of the particles and the strength of the inductive magnetic field H 0 are also involved in the error bound. The main conclusion of the analysis carried out on the "non-interacting assumption" is that the error is proportional to ε 4 /δ 3 . This behavior has also been observed in numerical experiments presented in the paper.

When the "non-interacting assumption" is deficient, we have proposed in the second part of the paper a method to compute a better approximation of the magnetic potential by taking into account pairwise magnetic field interactions between particles that enters in a general framework for computing the scalar magnetic potential as a series expansion. Numerical computation of pairwise magnetic field interactions relies on the evaluation of spherical harmonics series. Each pairwise interaction between two particles requires the evaluation of a spherical harmonics series. Fortunately, the series coefficients tend very quickly toward zero and in practise only a dozen of coefficients need to be evaluated. Moreover, computation of the series coefficient relies on the evaluation of spherical harmonics integrals for which fast spherical Fourier algorithms are available [START_REF] Kunis | Fast spherical fourier algorithms[END_REF]. And lastly, the method has good parallelisation properties which suggests that an efficient simulation software could be developped to deal with e.g. the applications in MRI quoted in the introduction. We have also obtained error bounds in terms of the above mentioned key parameters for the approximation of the magnetic potential taking into account pairewise interactions. It's well known (see [5, page 119]) that functions belonging to W 1 0 (R 3 ) admit a trace on the surfaces Σ j , j ∈ {1, ..., N }, and that the trace defines a continuous mapping from W 1 0 (R 3 ) onto H 1 2 (Σ j ). From Sobolev's embedding results, we also deduce that

W 1 0 (R 3 ) is continuously injected in L 1 (Σ j ): ∀j ∈ {1, ..., N } ∃C j > 0 ∀ϕ ∈ W 1 0 (R 3 ) ϕ L 1 (Σ j ) C j ϕ W 1 (R 3 ) .
The following proposition provides the constants C j when the surfaces Σ j are balls with radius ε. Proposition 8. For all j ∈ {1, ..., N } and for all ϕ ∈ W 1 0 (R 3 ) we have

ϕ L 1 (Σ j ) A 1 ε 3 2 ϕ W 1 (R 3 ) , (85)
where

A 1 = √ 4π.
Proof. We first consider the case where Σ j is the unit sphere S 2 centered at the origin and with radius 1. By a density argument, we can consider functions ψ in C ∞ 0 (R 3 ). In the spherical coordinates system (r, θ, η), we have

ˆS2 |ψ(1, θ, η)| sin θ dθdη ˆ∞ 1 ˆS2 |∂ r ψ(r, θ, η)| sin θ dθdηdr ˆ∞ 1 ˆS2 1 r 2 sin θ dθdηdr ˆ∞ 1 ˆS2 |∂ r ψ(r, θ, η)| 2 r 2 sin θ dθdηdr. It follows that ψ L 1 (S) √ 4π ψ W 1 (R 3 ) .
(86) In the case where Σ j is a sphere centered in c j and with radius ε > 0, we use the change of variable y ∈ S 2 → x = c j + εy ∈ Σ j . More precisely, for ϕ ∈ W 1 0 (R 3 ) we define ψ such that ψ(y) = ϕ(x). We have

ϕ L 1 (Σ j ) = ε 2 ψ L 1 (S) and ϕ W 1 (R 3 ) = √ ε ψ W 1 (R 3 ) . (87)
The estimate (86) together with (87) directly implies (85).

The Sobolev injection indicates that H 1 2 (Σ j ) ⊂ L 4 (Σ j ) and we also have the following finer result.

Proposition 9. For all j ∈ {1, ..., N } and for all ϕ ∈ W 1 0 (R 3 ) we have

ϕ L 4 (Σ j ) A 2 ϕ W 1 (R 3 ) , (88)
where

A 2 = 2 3 3 8 √ π .
Proof. As in the proof of Proposition 8, we only have to deal with the case where Σ j is the unit sphere S 2 . We also assume that ϕ belongs to C ∞ 0 (R 3 ) since the estimate for ϕ ∈ W 1 0 (R 3 ) will follow by density. Using the spherical coordinates (r, θ, η), we have

ϕ 4 L 4 (S 2 ) = ˆS2 |ϕ(1, θ, η)| 4 sin θ dθdη 4 ˆ∞ 1 ˆS2 |ϕ(r, θ, η)| 3 |∂ r ϕ(r, θ, η)| sin θ dθdηdr 4 ˆ∞ 1 ˆS2 |ϕ(r, θ, η)| 6 sin θ dθdηdr ˆ∞ 1 ˆS2 |∂ r ϕ(r, θ, η)| 2 sin θ dθdηdr 4 ϕ 3 L 6 (R 3 ) ϕ W 1 (R 3 )
. Using the following Sobolev embedding result in R 3 , see [7, page 26],

ϕ 3 L 6 (R 3 ) 4 3 3 2 1 2π 2 ∇ϕ 3 L 2 (R 3 ) . ( 89
)
we deduce that ϕ L 4 (S 2 ) 2 3

3 8 √ π ϕ W 1 (R 3 ) . 1 |c j -x| p 137 2 δ p .
(91)

Proof. Let i ∈ {1, . . . , N }, x ∈ R 3 such that dist(x, Σ i ) δ 2 .
For p 1, we split the sum

S (p) i (x) = N j=1 j =i 1 |c j -x| p
into two parts. The first one concerns the indices j such that the particles Ω j with center c j are at a distance lower than 2δ from Ω i . The second one concerns the indices j such that the particles Ω j are not in this neighborhood. Moreover, since the volume of the ball B(c j , δ

2 ) is π 6 δ 3 , we have

S (p) i (x) = N j=1, j =i |c j -c i |<2δ 1 |c j -x| p + N j=1, j =i |c j -c i | 2δ 6 πδ 3 ˆB(c j , δ
2 )

dy |c j -x| p . ( 92
)
The number of terms in the first sum in the RHS of (92), denoted S

i,1 in the sequel, is bounded independently of N . Actually, this number is lower than the cardinal number C of the set of points in B(0, 2) \ B(0, 1), such that the pairwise distance is larger than 1. By volume considerations, this number is found to be less than the quotient obtained by dividing the volume of B(0, 5 2 ) \ B(0, 1 2 ) by the volume of B(0, 1 2 ) that is to say 124. Each term in S i,1 is smaller than 1 (δ/2) p so that

S (p) i,1 124 2 δ p . (93)
Let us now obtain an estimate of the second sum in the RHS of (92), denoted S (p) i,2 in the sequel. First of all, from the triangular inequality, we note that for all y ∈ B(c j , δ

2 ),

|y -c i | |y -c j | + |c j -x| + |x -c i | δ 2 + |x -c j | + δ 2 = |x -c j | + δ. It follows that |x -c j | |y -c i | -δ and S (p) i,2 6 π δ 3 N j=1, j =i |c j -c i | 2δ ˆB(c j , δ 2 ) dy (|y -c i | -δ) p 6 π δ 3 ˆB(c i ,2)\B(c i , 3 2 δ) dy (|y -c i | -δ) p (94)
since under our assumptions for all j = 1, . . . , N , j = i, we have B(c j , δ

2 ) ⊂ B(c i , 2)\B(c i , 3 2 δ). The last integral in (94) can be expressed in spherical coordinates, and then it can be explicitly evaluated. We obtain

S (p) i,2 24 δ 3 ˆ2 3 2 δ r 2 dr (r -δ) p = 24 δ 3 ˆ2-δ δ 2 (s + δ) 2
s p ds := 24 δ 3 I p .

We now distinguish the different possible values of p envisaged in the proposition since the integral takes different closed form expressions depending on the value of p.

-Case p = 3. We have

I 3 = ˆ2-δ δ 2 (s + δ) 2 s 3 ds = 48 -56δ + 15δ 2 2(2 -δ) 2 + ln 2 δ (2 -δ) .
From ( 92) and (93), we deduce that

S (3) i (3I 3 + 124) 2 δ 3 .
One can show that the mapping δ → I 2 is bounded by 18+ln 2+ln 2 δ . This conclude the proof of estimate (90).

-Case p 4. We have

I p = ˆ2-δ δ 2 (s + δ) 2 s p ds = s 3-p 3 -p + 2δs 2-p 2 -p + δ 2 s 1-p 1 -p 2-δ δ 2 13 3 2 δ p-3
.

We deduce that

S (p) i,2 13 2 δ p .
From (92) and (93), we deduce the estimate (91).

One of the important feature of the estimates given in proposition [START_REF] Jin | The finite element method in electromagnetics[END_REF] is that they are optimal. More precisely, the behaviour of the sums with respect to the variable δ is correctly estimated. This remark can be understood by looking at the case of a large number of particles (N → +∞) uniformly distributed in the ball B(0, 1), with one particle located at the origin. To illustrate the situation, we consider the case p = 3. Let assume that c 1 = (0, 0, 0) for simplicity. Since δ measures the distance between the centers of two "adjacent" particles, these centers have for coordinates

(k 1 δ, k 2 δ, k 3 δ) where (k 1 , k 2 , k 3 ) ∈ Z 3 with δ 2 (k 2 1 + k 2 2 + k 2 3 ) 1.
It is clear that the "worst" case, i.e. the one that makes the largest sum S

(3) i , is the case i = 1. Thus, let us consider S

(3) 1 (0). For all j ∈ {2, . . . , N }, we have

S (3) 1 (0) = N j=1 j =1 1 |c j | 3 = (k 1 ,k 2 ,k 3 )∈Z 3 1 k 2 1 +k 2 2 +k 2 3 1 δ 2 1 (δ k 2 1 + k 2 2 + k 2 3 ) 3 = 1 δ 3 (k 1 ,k 2 ,k 3 )∈Z 3 1 k 2 1 +k 2 2 +k 2 3 1 δ 2 1 (k 2 1 + k 2 2 + k 2 3 ) 3 2
.

When N is large, the last sum is the approximation by the mid-ordinate quadrature rule of the integral

I = 1 δ 3 ˆB(0, 1 δ )\B(0,1) dx (x 2 1 + x 2 2 + x 2 3 ) 3 2
.

This integral can be easily evaluated using spherical coordinates and we obtain This shows that the estimate given in proposition 10 in the case p = 3 is optimal in the following sense: for δ = o(1) (that precisely corresponds to the case of many particles uniformly distributed in the ball B(0, 1)), estimate (90) reads

I = 4π δ 3 ˆ1 δ 1 dr r = 4 π δ 3 ln 1 δ . (95)
S (3) 1 (0) a δ 3 ln 1 δ ,
the constant a being independent of N , δ and ε, whereas in the present example, using (95), we have

S (3) 1 (0) b δ 3 ln 1 δ ,
where b is another constant independent of N , δ and ε.

Appendix C. Explicit solution to the reference problem for the first term of the expansion

In this appendix, we detail the way the explicit expression of the solution ϕ ref ∈ W 1 0 (R 3 ) to problem ( 14) is obtained. C.1. Variable separation method. In spherical coordinates, the Laplace equation is separable. We are looking for a solution ϕ ref to problem [START_REF] Kiselev | Analytical model of susceptibility-induced MR signal dephasing: Effect of diffusion in a microvascular network[END_REF] in the form

ϕ ref (r, θ, η) = u(r) × v(θ) × w(η) (96)
where u, v and w denote three complex valued functions of the real variable. The angles η and θ are the azimuthal angle and the polar angle respectively, see Fig. 10. In spherical coordinates, the Laplace operator reads

∆ϕ ref = 1 r 2 ∂ ∂r r 2 ∂ ∂r ϕ ref + 1 r 2 1 sin(θ) ∂ ∂θ sin(θ) ∂ ∂θ ϕ ref + 1 sin 2 θ ∂ 2 ∂η 2 ϕ ref . (97)
The Laplace equation ∆ϕ ref = 0 lends itself to being separated into a system of three ordinary differential equations: upon substituting ansatz (96) in the Laplace equation, one finds that for all r ∈ R *

+ , θ ∈]0, π[ and η ∈]0, 2π[, 1 u(r) d dr r 2 u (r) = 1 sin 2 (θ) 1 v(θ) sin(θ) d dθ sin(θ) v (θ) + w (η) w(η) . (98)
The left hand side of the equation depends only on the variable r whereas the right hand side depends only on the angular variables θ, η. It follows that the two terms are constant and equal to the same constant p 2 , where p denotes a complex number. Thus, from the left hand side of (98), we deduce that for all r ∈ R *

+ 1 u(r) d dr r 2 u (r) = p 2
and, from the right hand side of (98), we deduce that for all θ ∈]0, π[ and η ∈]0, 2π[,

1 v(θ) sin(θ) d dθ sin(θ) v (θ) -p 2 sin 2 (θ) = w (η) w(η) .
Once again, the left hand side of the equation depends only on the polar variable θ whereas the right hand side depends only on the azimuthal variables η. It follows that the two terms are constant and equal to the same constant q 2 , where q denotes a complex number. Finally, we conclude that the three unknown functions u, v and w satisfy for all r ∈ R * + , θ ∈]0, π[ and η ∈]0, 2π[,

r 2 u (r) + 2r u (r) -p 2 u(r) = 0 (99a) 1 sin(θ) d dθ (sin(θ) v (θ)) + (p 2 - q 2 sin 2 (θ) ) v(θ) = 0 (99b) w (η) + q 2 w(η) = 0 (99c)
where the real parameters p and q are the separation constants.

C.2. Angular dependency. Since the geometrical and physical properties of the ball are assumed to be independent of the azimuthal angle η, it is necessary for w to be a periodic function with period 2π. Then, the solution to equation (99c) reads

w(η) = C 1 e imη + C 2 e -imη (100)
where the separation constant q coincides with an integer m ∈ Z and C 1 and C 2 denote two complex constants. The change of variable τ = cos(θ) in equation (99b) leads to the equation

(1 -τ 2 ) v (τ ) -2τ v (τ ) + (p 2 - m 2 1 -τ 2 ) v(τ ) = 0 (101)
where the new unknown function v is defined by the relations v(θ) = v(cos(θ)). This equation is referred in the literature as the general Legendre equation. The solutions to the linear second order ordinary differential equation (ODE) (101) are the so-called hypergeometric functions, see [1] chp. 15. This equation has non-zero solutions that are nonsingular at ±1 if and only if p 2 = ( + 1) with being a non negative integer such that |m| . In such a case, the solutions to equation (101) are the associated Legendre functions, see [START_REF] Sansone | Orthogonal Functions[END_REF]1], defined for m 0 by where Y m is known as the Spherical Surface Harmonics of degree and order m, see [START_REF] Sansone | Orthogonal Functions[END_REF].

The normalization constant C m is taken to be

C m = (2 + 1) 4π ( -m)! ( + m)! . ( 104
)
With this normalization convention, the complex conjugate of the Spherical Surface Harmonics of degree and order m is Y m (θ, η) = (-1) m Y -m (θ, η). The Spherical Surface Harmonics form a complete set of orthonormal functions and thus they form an orthonormal basis of the Hilbert space of square-integrable functions on the unit sphere [START_REF] Sansone | Orthogonal Functions[END_REF]. Namely, on the unit sphere any square-integrable function can be expanded as a linear combination of Spherical Surface Harmonics as:

h(θ, η) = +∞ =0 m=- h m Y m (θ, η) (105)
where the equality holds in L Two linearly independent solutions to the linear second order ordinary differential equation (106) are r → r and r → r --1 . We deduce that the radial function u can be expressed as u(r) = α r + β 1 r +1 where α and β denote two complex constant numbers.

We conclude that the general solution of the Laplace equation in spherical coordinates reads

ϕ ref (r, θ, η) = +∞ =0 m=- α m r + β m 1 r +1 Y m (θ, η) (107)
where α m and β m are complex numbers.

Moreover, since the magnetic potential must tend to zero at infinity, the solution in B(0, ε) reads

ϕ ref (r, θ, η) = +∞ =0 m=- β m 1 r +1 Y m (θ, η). (109)
The Legendre polynomial P 0 = P is bounded by 1 and we have the following bounds for the Associated Legendre functions, see [START_REF] Lohöfer | Inequalities for the associated Legendre functions[END_REF],

∀m ∈ {1, . . . , } ( -m)! ( + m)! max For 1, for m ∈ {1, . . . , } and for all x ∈] -1, 1[, we have [1,[START_REF] Olver | NIST Handbook of Mathematical Functions[END_REF]:

∂ x P 0 (x) = ∂ x P (x) = - 1 √ 1 -x 2 P 1 (x) ∂ x P m (x) = 1 2 √ 1 -x 2
( + m)( -m + 1) P m-1 (x) -P m+1 (x) .

It follows that, for 1, for m ∈ {1, . . . , } and for all θ ∈ [0, π], we have We conclude the proof using 2 + 1 3 and ( + 1) + π 2 ( + 1) 2 2(1 + π) 2 . Proposition 13. For ∈ N, the family of spherical harmonics (Y m ) -m satisfies: Following the same steps as in the proof of Proposition 12, we conclude that m=- ≈ 0.3439.

∀(θ, η) ∈ [0, π] × [0, 2π]
1 sin(θ) ∂ η Y m (θ, η) √ 2 + 1 2 √ 2π m=- ( -m)( -m -1) √ 2 + 1 2 √ 2π m=- ( -m) (2 + 1)
Proof. We proceed in two steps. First, we will compute the expression of ∂ n G(x , •)

L 4 3 (S 2 )
for all x ∈ S 2 . Then we will show that for all x ∈ R 3 , denoting by x ∈ S 2 the projection of x onto S 2 , we have

∀y ∈ S 2 \ {x } |∂ n G(x, y)| |∂ n G(x , y)|. (121)
For all x ∈ S 2 and 1 < p < 2, let us consider the quantity

A = ˛S2 n • (x -y) 4π|x -y| 3 p dy
where n denotes the outward unit normal to S 2 at y ∈ S 2 . In order to compute A, let us consider a frame centered on the point x and the following parametrization of the sphere S 2 :

(s, t) ∈ B(0,

1) ⊂ R 2 -→   s t 1 ± √ 1 -s 2 -t 2   ∈ S 2 ,
where the sign ± is introduced to describe the two hemispheres. In this frame, the normal unit vector n to the sphere S 2 is given by

n =   s t ± √ 1 -s 2 -t 2   .

Figure 2 .

 2 Figure 2. Notation for the translation of the coordinates system.

Proposition 3 .

 3 The function ψ introduced in (20) belongs to W 1 0 (R 3 ) and satisfies    ∆ψ = 0 in Ω 1 , . . . , Ω N and Ω c µ ∂ψ ∂n = [µ] g 0 across Σ 1 , . . . , Σ N

Figure 3 .

 3 Figure 3. Effect of the translation of the coordinates system for several particles.

Remark 4 . 2 (

 42 When 2ε δ becomes small, since Li -3 z) ∼ 0 z, estimate (81) reads ϕ -N j=1

5. 1 .

 1 Convergence toward zero of the series coefficients. In order to illustrate the fast convergence to 0 of the coefficients γ m given by (83) as tends to infinity, we have depicted in Fig.4the values |g m | for m = -, . . . , and = 0, . . . , 15 for µ p = 10, ε = 1 and u = c k -c j = (1, 1, 1) . We have also depicted in Fig. 5 the decimal logarithm of max m=-,..., |g m | for = 1, . . . , 29. The slope of the straight line is approximately -0.2. This means that from to + 1 the maximal value max m=-,..., |g m | is divided by 1.6.

(Figure 4 .Figure 5 .

 45 Figure 4. Values of |g m | for m = -, . . . , and = 0, . . . , 16 for µ p = 10, ε = 1 and u = c k -c j = (1, 1, 1) .

21

 21

Figure 6 .

 6 Figure 6. Values of E as a function of the distance δ between the two spherical particles in a log-scale. The slope of the solid line is -3.

Figure 7 .

 7 Figure 7. Values of E as a function of the radius ε of the two spherical particles in a log-scale. The slope of the solid line is 4.

c 2 .

 2 These three quantities were computed over a grid of 50 × 100 points in the plane of reference. Computation time (CPU time) was 51.03 s under Matlab (R2018a) on an Intel i5 Quad Core desktop computer.

Figure 8 .

 8 Figure 8. Two particles. From top to bottom : first term in the series expansion of the magnetic potential, second term in the series expansion and approximation of the magnetic potential obtained by summing the first two terms of the series expansion. From left to right : three-dimensional shaded surface representation and two-dimensional representation of the three quantities.

Figure 9 .

 9 Figure 9. Cube of 64 equidistant particles. From top to bottom : first term in the series expansion of the magnetic potential, second term in the series expansion and approximation of the magnetic potential obtained by summing the first two terms of the series expansion. From left to right : three-dimensional shaded surface representation and two-dimensional representation of the three quantities.

Figure 10 .

 10 Figure 10. Notation for the spherical coordinates.

  v(τ ) = P m (τ ) := (-1) m (1 -τ 2 ) m 2 d m dτ m P (τ )(102)where P denotes the Legendre polynomial of degree defined for all x ∈] -The associated Legrendre functions are sometimes defined without the multiplicative constant (-1) m which is known as the Cordon-Shortley phase factor. Following the standard convention, we set for negative values of the orderP -m (x) = (-1) m ( -m)! ( + m)! P m (x) ∀x ∈] -1, 1[.For a fixed integer m, the associated Legrendre functions satisfy the following orthogonality conditions for all k, ∈ N such that 0 m < ˆ1 -1 P m k P m dx = 0 and ˆ1 -1P m P m dx = 2( + m)! (2 + 1)( -m)! .Finally, we find the functions product v(θ) w(η) in the form of Y m (θ, η) = C m P m (cos(θ)) e imη -

C. 3 .

 3 2 ([0, π] × [0, 2π], C) and h m = ˆ2π 0 ˆπ 0 h(θ, η) Y m (θ, η) sin(θ) dθ dη. Radial dependency. It remains to achieve the identification of the radial function u. It satisfies the following Euler differential equation deduced from (99a) r 2 u (r) + 2ru (r) -( + 1)u(r) = 0.(106)

C. 4 .

 4 Boundary condition. Since ϕ ref must be bounded for r = 0, the solution to problem[START_REF] Kiselev | Analytical model of susceptibility-induced MR signal dephasing: Effect of diffusion in a microvascular network[END_REF] in B(0, ε) readsϕ ref (r, θ, η) = +∞ =0 m=-α m r Y m (θ, η).

Proposition 12 .

 12 For ∈ N, the family of spherical harmonics (Y m ) -m satisfies:∀(θ, η) ∈ [0, π] × [0, 2π] m=-|∂ θ Y m (θ,Using the definition of the Spherical Harmonics, we get m=-|∂ θ Y m (θ, η)| 2 + 1 4π |∂ θ (P 0 (cos(θ)))| + 2 m=1 ( -m)! ( + m)! |∂ θ (P m (cos(θ)))| .

∂ 2 P 2 -m 2 .( + 1 ) 2 -m 2 = ( + 1) + 2 m=1( + 1 ) 2 -m 2 (

 2221222122 θ (P 0 (cos(θ))) = P 1 (cos(θ))∂ θ (P m (cos(θ))) = 1 m+1 (cos(θ)) -( + m)( -m + 1) P m-1 (cos(θ)) .Using relations (116), it follows thatm=-|∂ θ Y m (θ, η)|By comparison to the quadrature rectangle rule, we have m=-

For 1 , 2 √ 1 -x 2 P

 1212 ∂ η Y m (θ, η) Using the definition of the Spherical Harmonics, we express the sum in the LHS of (118) using the associated Legendre polynomials: for m ∈ {1, . . . , } and for all x ∈] -1, 1[, we have[1,[START_REF] Olver | NIST Handbook of Mathematical Functions[END_REF]:m m (x) = -P m+1 -1 (x) -( + m)( + m -1) P m-1 -1 (x).It follows that, for 1, for m ∈ {1, . . . , } and for all θ ∈ [0, π], we have m 2 sin(θ) P m (cos(θ)) = -P m+1 -1 (cos(θ)) -( + m)( + m -1) P m-1 -1 (cos(θ)).

2 .where A 3 =

 23 Appendix E. L p -norm of the normal derivative of Green's kernel Proposition 14. For all x ∈ R 3 , we have∂ n G(x, •) ∂ n G((1, 0, 0), •)

1 |x -c j | 2 .

Appendix B. Some results on the sum of powers of the inverse distance

The results given in the following proposition are used in the study of the behavior of the various terms in the expansion of the magnetic potential.

Proposition 10. Let Ω 1 , . . . , Ω N be N spherical particles in B(0, 1) ⊂ R 3 , centered respectively at c j , j ∈ {1, . . . , N }, with the same radius ε > 0. Let δ be the minimal distance between two particles as defined in (1). For all i ∈ {1, . . . , N } and for all x ∈ R 3 such that dist(x, Σ i ) δ 2 , we have

where B 3 = 142 + 3 ln 2 ≈ 144.1. Moreover, for all p 4, we have

The constants α m and β m are determined by the following two conditions across S(0, ε) in ( 14)

In the situation considered here, the source term h depends only on the polar angle θ. This implies that the solution ϕ ref will not depend on the azimuthal angle η and reads

The interface conditions (110) give rise to the following linear system satisfied by the unknown coefficients α 0 and β 0 , ∈ N

Solving this linear system for α 0 and β 0 gives

and α 0 = β 0 = 0 for all = 1. Finally, we conclude that

The gradient of ϕ ref is given in spherical coordinates by

Appendix D. Some results on the Spherical Surface Harmonics Proposition 11. For ∈ N, the family of spherical harmonics (Y m ) -m satisfies:

Proof. From the definition of Spherical Surface Harmonics, we have

Since |y| 2 = 2n • y, we have A = A + + A -where

Using the change of variable u = 1 ± √ 1 -r 2 , we deduce that

For p = 4 3 , by evaluating A by quadrature, we found A ≈ 0.24093027. Finally,

In order to prove (120), we consider x ∈ R 3 and we denote by x its projection onto S 2 .

Considering the frame previously introduced and centered on x , the coordinate of x reads (0, 0, -d) where d = dist(x, S 2 ) 0. For any point y ∈ S 2 , we deduce from the relation

.

Similarly, one can show that

Since d 0 and |x -y| |x -y| we directly deduce that

This last inequality implies that ∂ n G(x, •)

. We conclude the proof of Proposition