Mathematical imaging using electric or magnetic nanoparticles as contrast agents

We analyse mathematically the imaging modality using electromagnetic nanoparticles as contrast agent. This method uses the electromagnetic fields, collected before and after injecting electromagnetic nanoparticles, to reconstruct the electrical permittivity. The particularity here is that these nanoparticles have high contrast electric or magnetic properties compared to the background media. First, we introduce the concept of electric (or magnetic) nanoparticles to describe the particles, of relative diameter $\delta$ (relative to the size of the imaging domain), having relative electric permittivity (or relative magnetic permeability) of order $\delta^{-\alpha}$ with a certain $\alpha>0$, as $0<\delta<<1$. Examples of such material, used in the imaging community, are discussed. Second, we derive the asymptotic expansion of the electromagnetic fields due to such singular contrasts. We consider here the scalar electromagnetic model. Using these expansions, we extract the values of the total fields inside the domain of imaging from the scattered fields measured before and after injecting the nanoparticles. From these total fields, we derive the values of the electric permittivity at the expense of numerical differentiations.


1.
Introduction. Let D be a small body in R 3 of the form D := δB + z, where B is open, bounded with a Lipschitz boundary, simply connected set in R 3 containing the origin, and z specify the locations of D. The parameter δ > 0 is the B-relative radius of D. It characterises the smallness assumption on the body D as compared to B.
Let us consider the magnetic permeability µ δ and electric permittivity δ of the form (1) µ δ (x) = µ 0 , x ∈ R 3 D, µ 1 , x ∈ D, where µ 0 , µ 1 , 1 are positive constants while 0 := 0 (x) is variable in a bounded domain, i.e. there exists Ω bounded such that 0 is a constant in R 3 Ω. Here, Ω is of the same order as the reference body B. Thus µ 0 , 0 denote the permeability and permittivity of the background medium and µ 1 , 1 denote the permeability and permittivity of the scatterers respectively. We are then interested in the following scattering problem: (3) where ω > 0 is a given frequency. The total field V has the form V := V I + V s where V I denotes the incident field and V s denotes the scattered waves. The above set of equations have to be supplemented with the Sommerfeld radiation condition on V s which we shall henceforth refer to as (S.R.C). Here, V describes components of the electric field 1 . Keeping in mind the positivity of the permeabilities, the above problem can be equivalently formulated as where κ 2 0 = ω 2 0 µ 0 , κ 2 1 = ω 2 1 µ 1 and V s denotes the scattered field. We restrict to plane incident waves V I , i.e. V I (x) := e iκ0x·d , where x ∈ R 3 , d ∈ S 2 . Here S 2 denotes the unit sphere in R 3 .
The scattering problem (4) is well posed in appropriate spaces and the scattered field V s (x, d) has the following asymptotic expansion: , |x| → ∞, 1 In case we have invariance of the model in one direction, this Helmholtz model describes the propagation of the component of the electric field which is orthogonal to that axis. To simplify the exposition, we stated it in 3D instead of 2D to avoid the log-type singularities of the corresponding Green's functions. withx := x |x| , where the function V ∞ (x, d) for (x, d) ∈ S 2 × S 2 is called the far-field pattern.
Problem. Our motivation in this work is the use of electromagnetic nanoparticles as contrast agent to image the background, i.e. to reconstruct the coefficient 0 in Ω. The measured data in this case is described as follows: , for a single direction d} when no nanoparticle is injected.
, for a single direction d} when one single particle is injected to a given point z ∈ Ω. This measurement is needed for all the points (or a sample of points) z in Ω.
If, in (3) or (4), we exchange the roles of (µ 1 , µ 0 ) with ( 1 , 0 ), where µ 0 is now locally variable while 0 is a constant, then V describes components of the magnetic field. For such a model, the goal is to reconstruct the magnetic permeabilty µ 0 in Ω from the corresponding measured data generated by the nanoparticles 2 .
Imaging using electromagnetic nanoparticles as contrast agents has drawn a considerable attention in the very recent years, see for instance [9,12,19]. To motivate it, let us first recall that conventional imaging techniques, such as microwave imaging techniques, are known to be potentially capable of extracting features in breast cancer, for instance, in case of relatively high contrast of the permittivity and conductivity, between healthy tissues and malignant ones, [13]. However, it is observed that in the case of a benign tissue, the variation of the permittivity is quite low so that such conventional imaging modalities are limited to be used for early detection of such diseases. In such cases, creating such missing contrast is highly desirable. One way to do it is to use electromagnetic nanoparticles as contrast agents, [9,12]. In the literature, different types of nanoparticles have been proposed with this application in mind. Let us cite few of them: 1. To create contrast in the permittivity, carbon nanotubes, ferroelectric nanoparticles and Calcium copper titanate are used, see [9]. Such particles have diameter which is estimated around 10 nm, or 10 −8 m, and have relative electrical permittivity of the order 10 for the carbon nanotubes, 10 3 for ferroelectric nanoparticles and around 10 6 for Calcium copper titanate, see [22]. If the benign tumor is located at the cell level (which means that our Ω is that cell), with diameter of order 10 −5 m, then the Ω-relative radius of the particles are δ = 10 −3 . Hence, the relative permittivity of the types of nanoparticles are estimated of the order δ −α where α is 1 3 for the carbon nanotubes, 1 for ferroelectric nanoparticles and around 2 for Calcium copper titanate. 2. The human tissue is known to be nonmagnetic. To create magnetic contrasts, it was also proposed in [9] to use iron oxide magnetic nanoparticles for imaging early tumors. Such material has relative permeability of the order between 10 4 and 10 6 , see [21], and hence of the order δ −α with α between 4 3 and 2. 3. If the tumor is a relatively malignant one and occupy bodies having diameter of order of few millimeters or even centimeters, then we end up with values of α of order 1 5 , 3 5 and 1 respectively for the aforementioned electric nanoparticles and between 2 3 and 1 for the iron oxide magnetic nanoparticles.
This shows that for the detection of the tumors using such nanoparticles, we can model the ratio of the relative electric permittivity and relative magnetic permeability as follows: Definition 1.1. We shall call (D, 1 , µ 1 ) an electromagnetic nanoparticle of shape ∂D, diameter of order O(δ), δ << 1, and permittivity and permeability 1 , µ 1 respectively.
Let us stress here that these electromagnetic nanoparticles have large relative speeds of propagation (or indices of refraction).
To highlight the relevance of each type of data that we use in stating the imaging problem, we briefly analyse their importance: 1. Use of scattered data before injecting any nanoparticle. Assume that we have at hand the farfields V ∞ (x, d) in the case whenx and d are taken in the whole S 2 . The problem with such data is well studied and there are several algorithms to solve it , see [15,16,17]. It is also known that this problem is very unstable. Precisely, the modulus of continuity is in general of logarithmic type, see [18]. 2. Use of scattered data before and after injecting one single nanoparticle at time.
In [4,5], and the references therein, an approach of solving the problem using data collected before and after localized elastic perturbation has been proposed. However, their analysis goes also for the case of injecting a single nanoparticle at a time (precisely moderate nanoparticles, i.e. both the electrical permittivity and the magnetic permeability are moderate in terms of the size of the particles (taking α = 0 above) ). It is divided into two steps. In the first one, from these far-fields, we can extract the total energies due to the medium 0 in the interior of the support of ( 0 − 1). In the second step, we reconstruct 0 (x) from these interior values of the energies. In contrast to the instability of the classical inverse scattering problem, the reconstruction from internal measurements is stable, see [2,6,14,20]. However, one needs to handle the invertibility properties of the energy mappings. One can think of using measurements related to multiple frequencies ω 1 , ω 2 , etc., see [1], to invert such mappings.
Let us emphasize here that if we use electric (or magnetic) nanoparticles, instead of moderate nanoparticles, then we can extract, from the corresponding far-fields, the pointwise values of the total fields and not only the energies, see Theorem 2.1. This is one reason why we introduced the concept of electric (or magnetic) nanoparticles. The second step is then to extract the coefficients 0 from these data. This can be done at the expense of a numerical differentiation. It is known that numerical differentiation is an unstable step but it is only moderately unstable. Finally, observe that, we use the backscattered data, before and after injecting the particles, in only one and arbitrary direction d. We found that using electric nanoparticles (see the definition above), the electric far-field data are enough for the reconstruction. However, if we use magnetic nanoparticles, then, in addition to the electric far-fields, the electric near-fields are also required. Similarly, using magnetic nanoparticles, the magnetic far-field data are enough for the reconstruction. However, if we use electric nanoparticles, then, in addition to the magnetic far-fields, the magnetic near-fields are also needed. Hence, as a moral, if we use electric (resp. magnetic) nanoparticles, it is better to use electric (resp. magnetic) far-fields. In both cases, the error of the reconstruction is fixed by the parameter α which is related to the kind of (relative permittivity or permeability of the) nanoparticles used. In addition, the distance from the place where to collect the near-field data to the tumor is also fixed by the kind of the nanoparticles used. More details are provided in Remark 1.
The remaining part of the paper is organized as follows. In Section 2, we state the asymptotic expansion of the far-fields and near-fields generated by high contrast electric permittivity or magnetic permeability. Then, we discuss how these expansions can be used to derive the needed formulas to solve the imaging problem. In section 3, we give the details of the proof of the main theorem stated in section 2.
2. The asymptotic expansion of the fields and application to imaging.
2.1. The asymptotic expansion with singular coefficients. Let us recall our model: , |x| → ∞ (S.R.C), where κ 2 0 = ω 2 0 µ 0 , κ 2 1 = ω 2 1 µ 1 and V s denotes the scattered field. In the case when 0 is a constant every where in R 3 , we use the notation u instead of V in (6). The solution u can then be represented as where u I denotes the incident wave and φ, ψ are appropriate densities (see [7]).
When 0 is variable, we can represent V as and U t is the solution of the following scattering problem; From (7), we get where (G κ0 ) ∞ (x, y) is the far-field of the Green's function G κ0 (x, y) of the problem (8). The mixed reciprocity condition implies that From (9), we have (10) and from (7), we have where d(x, z) is the Euclidean distance from x to z.
In the following theorem, we state the asymptotic expansion of the far-fields and the near-fields in terms of δ at the first order.
Theorem 2.1. Assume that D is a Lipschitz domain, µ 0 , µ 1 and 1 are positive constants while 0 is a Lipschitz regular function and equals 1 outside a bounded and measurable set Ω. In addition, the coefficients 0 and µ 0 are assumed to be uniformly bounded in terms of the diameter δ of D. Then, at a fixed frequency ω, we have the following asymptotic expansions for the far-fields and the near-fields as δ → 0: 1. Far-fields: In the case of electric nanoparticles, where 0 < α < 1, and in the case of magnetic nanoparticles, where 0 < α < 1.

Near-fields:
in the case of electric nanoparticles, and in the case of magnetic nanoparticles.
is the polarization tensor and λ = 1 2 µ0+µ1 µ0−µ1 . 2.2. Application to imaging. For simplicity of exposition, we consider B to be a ball.
1. Using electric nanoparticles. In this case, we deduce from (12) that In particular forx = −d, we get Assume that we have at hand the farfield measured before using the nanoparticles, i.e. U ∞ (x, d), and after using the nanoparticles, i.e. V ∞ (x, d). Then from (17), we can compute ±U t (z, d) with an error of the order O(δ α ). By numerical differentiation we can compute ∆U t (z), in the regions where U t (z, d) does not vanish, and hence Observe that, we use the backscattering data for only one incident direction d.

2.
Using magnetic nanoparticles. In this case, we deduce from (13) where m 0 is related to the polarization tensor M as M = m 0 |B|Id + O(δ α ). Following [8] (Page 83) and using the fact that λ = 1 2 µ0+µ1 µ0−µ1 , we can derive using the facts that when x is close to z, see [3], Assume that we can measure the far-fields and the near fields before as well as after using the nanoparticles. (a) From (19), we can reconstruct with an error of the order δ α using one direction of incidence d and the corresponding backscattering directionx = −d.
Thus from (20), we can reconstruct ∇U t (z, d) with an error of the order Hence, we can reconstruct ( 1 − 0 ) (U t (z, d)) 2 with an error of the order δ min{α,1−α} if we use the far-field data and the near fields collected at a distance, from the tumor location, of the order d(x, z) ∼ O(δ α ). Let now K be a smooth region where we have sent the nanoparticles. Knowing ∇U t (z, d), for a sample of points z ∈ K, we can recover −ω Combining instance, in the case 1 > 0 (z), z ∈ Ω, we can consider the second-order polynomial equation It is easy to see that this equation has a positive root 0 since 1 > 0 . In both cases, the error (which is O(δ α ), where α ∈ (0, 1), using electric nanoparticles, or O(δ min{α,1−α} ), for α ∈ (0, 1), using magnetic nanoparticles) propagates through the numerical differentiation which is a moderately unstable step, i.e. with a loss of a polynomial rate of the error.

Remark 1.
We observe that to get the reconstruction of the coefficients 0 1. using electric nanoparticles, we need to measure the far-fields before and after the injection of these particles. However, we only need one single backscattered direction d ∈ S 2 . 2. using magnetic nanoparticles, we need to measure the far-fields before and after the injection of these particles in one single backscattered direction d ∈ S 2 and the near-fields in three different points x j , j = 1, 2, 3 such that the unit vectors d(xj ,z) are linearly independent. In both cases, the error of the reconstruction is fixed by the parameter α which is related to the kind of (relative permittivity or permeability of the) nanoparticles used. The distance from the place where to collect the near-field data to the tumor (which is of the order δ α for α ∈ (0, 1)) is also fixed by the kind of magnetic nanoparticles used.
3.1. The case when 0 is a constant. First, we recall that in the case where 0 is a constant every where in R 3 , we use the notation u instead of V and that the solution u can be represented as By the jump relations, the densities φ and ψ are solutions of the system of integral equations Using the fact that [− 1 2 Id + (K κ0 D ) * ] is invertible, from the second equation of (23), we have From the first identity, we have Notation. For φ ∈ L 2 (∂D), ψ ∈ L 2 (∂B), we shall denote byφ,ψ the functionŝ Proof. Let us first assume that To see this, we write is an isomorphism with a bounded inverse in terms of λ, |λ| ≥ 1 2 , see [8]. Hence (29) Combining the estimates from the above two cases, we can now prove the desired result. To do so, we first write [10,11]) and Remark 2. In case κ 0 is not constant, we shall see later that and not of order δ 2 . Therefore in the above lemma, we have to choose α accordingly. But surely this holds if β < 1.
Then from (27) it follows that Let us set W κ1 : Lemma 3.3. We have the following behavior of W κ1 : Proof. We note that H := W κ1 − u I satisfies the elliptic equation Multiplying the first equation of (35) byH and by applying the integration by parts formula, we obtain: Applying the Poincare inequality for the first term on L.H.S, we obtain (37) Now since the first Dirichlet eigenvalue is the same as the sharp Poincare constant C(D), by the Faber-Krahn inequality, we have where D * denotes the ball of same volume as D, δ 2 (diam(B)) 2 j 2 1 2 ,1 =: and j p,1 denotes the first positive zero of the Bessel function (of first kind) J p . Using the inequality (38) in (37) and by applying the Cauchy Schwartz inequality, we obtain: Now, rewrite (35) as, By scaling we can rewrite (40) on the reference body B as, By applying the elliptic regularity theorem, there exists a positive constant C 2 (B) such that Hence, we are done with proving (33). For any function Ξ ∈ H 2 (D), D ⊂ R N , one can prove the following by rescaling In our case, we have N = 3, Ξ = H and so we have the following inequality: By making use of the right inequality of (44) in (42), we get (45) Hence, we are done with proving (34).
Proof. We can write which implies As ∇ ξĤ = δ∇ x H, then From the estimate (see (42)) Using trace theorem, it then follows that Using this in (46), we have the required estimate.
We can rewrite the first identity in (23) as and the second identity in (23) as

Proof. Using Lemma 3.4, we can write
Let us also note that u I , ∇u and O(δ 2 ) in L 2 (∂D) respectively, and ignoring the terms of order O(δ 3 ) in (49), we obtain (51) 1 Keeping in mind the linearity of the identities, we split (50) and (51) into seven problems as following.
In the following, we shall treat each of these problems separately. We begin by noting the scaling property of the inverse of the single-layer potential.
Proof. Let f ∈ H 1 ♦ (∂D), that is, the meanf = 0. Using the scaling properties (see [10]) and surface Poincare inequality, we can write Bφ =f on ∂B, and we can write We rewrite the first identity in P1 as Let us also note that both H1 and H2 are harmonic in x, and hence satisfies Therefore by maximum principle, it follows that H1 − H2 = 0 in D, that is, Taking the normal derivative from the interior of the domain D, we then have This can also be written as Combining this with the second identity of P1, we then obtain µ0−µ1 . Now since µ 0 = µ 1 , the operator [λId + (K 0 D ) * ] is invertible and hence we can write (60) Next we note that since φ 1 satisfies the first identity in (52), we can also write Now using the fact ν i (x) = ∂ ν x i and the Green's formula, we have (61) 3.1.2. Treatment of Problem P2. As in the case of Problem P 1, we use the harmonicity of the right hand side of the first identity in (53) and proceed as following. Let H3(x) = 0 on ∂D.
Therefore by the maximum principle, we infer that H3(x) = 0 in D, and hence by taking the normal derivative from the interior we obtain (63) Using this in the second identity of (53), we have (64) 1 µ 1 Next we note that we can also write (65) Then using Green's formula, we have Next we estimate ∂D (x − z) i φ 2 (x) ds(x) for |j| = 1. Using Cauchy-Schwarz inequality, we have and from (64), using Minkowski's inequality, we have Thus we need to estimate ∂ ∂ν (x − z) j L 2 (∂D) . But and therefore we obtain that 3.1.3. Treatment of Problem P3. In this case, using the second identity of (54) we write and hence using Green's identity, we have Let us now estimate ∂D (x − z) i φ 3 (x) ds(x) for |i| = 1. Using Cauchy-Schwarz inequality, we have and also we have already seen that ( , that is, c is the average of f over the boundary ∂D of D. Then where C 1 is independent of δ. Nowf whereby we obtain that c = C 2 δ 2 , where the constant C 2 is now independent of δ. Now using the invertibility of (S 0 D ) −1 , we can rewrite the first identity in (54) as . Using this in the second identity of (54), we obtain whereby we obtain , and (S 0 D ) −1 is uniformly bounded with respect to δ when acting over the mean-zero function f −c (see Lemma 3.8). Therefore φ 3 L 2 (∂D) = O(δ 2 ) and hence 3.1.4. Treatment of Problem P4. In this case, we note that the right hand side of (55) is a constant, and hence independent of x. Therefore we can again use the harmonicity of both sides of (55) together with maximum principle to obtain the identity . Taking the normal derivative from the interior of the domain, we then obtain Using the second identity in (55), we then obtain 3.1.5. Treatment of Problem P5. In this case, using the second identity in (56), we can write Using Green's formula, we can then write Let us writef (x) := −(S κ1 D ψ(x) − S 0 D ψ(x)) and Then proceeding as in the case of Problem P 3, we can write Then where S κ1,δ Bψ (ξ) = ∂B e iκ 1 δ|ξ−η| 4π|ξ−η|ψ (η) ds(η). Using Cauchy-Schwarz inequality, we can estimate Next we note that and therefore S κ1,δ Bψ − S 0,δ Bψ L 2 (∂B) ≤ Cδ −1 δ 1− α 2 = Cδ − α 2 . Using this in (73), we get where C is independent of δ. We then have and hence φ 5 L 2 (∂D) ≤ C µ0 µ1 δ 2−α . Therefore, using Cauchy-Schwarz inequality, it follows that 3.1.6. Treatment of Problem P6. In this case, using the second identity in (57), we can write Using Green's formula, we can then write (76) Proof. From the first identity in (23), we have For ξ ∈ B, let us denote byĜ(ξ) := G(δξ + z). ThenĜ satisfies where S κ1,δ Bψ (ξ) := ∂B e iκ 1 δ|ξ−η| 4π|ξ−η|ψ (η) ds(η). Multiplying both sides of (78) byĜ and integrating by parts, we obtain Using Poincare's inequality, we get Using Holder's inequality on the right hand side, we then have For δ > 0 small enough, we have (C(B) − δ 2 κ 2 0 ) > 0, and hence For ξ ∈ B, Therefore, we have By using Holder's inequality, we can write But this implies We next estimate D S κ0 D φ(x) dx and D u I (x) dx. Using a Taylor series expansion around the centre z, we have the estimate Using (84) and (85) in (83), the required result follows.
Therefore from (76), we have From the first identity of (57), we have (using harmonicity and maximum principle and the fact that the left-hand side is a constant) Using this, from the second identity, we obtain This then implies that φ 6 L 2 (∂D) = O( µ0 µ1 δ 2−α ) and hence using Cauchy-schwarz inequality, we have 3.1.7. Treatment of Problem P7. Integrating the second identity of (58) over ∂D, using the fact that ∂D (K 0 D ) * ψ 7 (x) = − 1 2 ∂D ψ 7 (x), we can conclude that For functions F ∈ H 1 (∂D) and G ∈ L 2 (∂D) in the right-hand side of (58), we can estimate the L 2 -norm of φ 7 as follows. From the second identity in (58), we can write Again, from the first identity it follows that Combining the above two, we can write Now we can use the facts that − 1 2 are uniformly bounded with respect to δ and (S 0 D ) −1 scales as 1 δ to prove the estimate Using this in our setup, we obtain and hence 3.1.8. Asymptotics for the near-field and far-field pattern. In this section, we sum up the results obtained in the previous sections to find the asymptotics for the far field pattern. From the analysis in the previous sections, we have already seen that for α ∈ (0, 1), where |j| = 1.
We recall that the total field is of the form Thus the scattered field is given by 3.1.9. The near-fields. To estimate the near fields, we first derive the following estimate (91) where d := d(x, z). Then by (90), we deduce that 3.1.10. The far-fields. The far field pattern u ∞ (x) can be written as Therefore,

3.2.1.
Comparison of estimates to the homogeneous background case. We shall now estimate S κ0 D − S 0 D L(L 2 (∂D),H 1 (∂D)) and (K κ0 D ) * − (K 0 D ) * L(L 2 (∂D),L 2 (∂D)) . In contrast to the case of homogeneous background, we are able to derive slightly weaker estimates in this case.
The problem (95) can be dealt with just as in the previous case recalling that κ 0 is Lipschitz continuous. To study the problem (96), we proceed as follows.
Summing up the discussion above, we can therefore conclude that S κ0 D − S 0 D L(L 2 (∂D),H 1 (∂D)) = O(δ 2 |logδ|). The proof for the other case follows similarly.
It is clear that the new P1, P2, P3, P4 can be dealt in the same way as the ones before (only replacing U I (z) by U t (z)). Hence, we have the corresponding estimates as derived before.
Let us now discuss the way to handle P5, P6 and P7.