Photo-acoustic inversion using plasmonic contrast agents:

www.oeaw.ac.at

www.ricam.oeaw.ac.at

Photo-acoustic inversion using plasmonic contrast agents: The full Maxwell

model

A. Ghandriche, M. Sini

RICAM-Report 2021-40

Photo-acoustic inversion using plasmonic contrast agents:

The full Maxwell model

Ahcene Ghandriche^a,∗, Mourad Sini^b

aRICAM, Austrian Academy of Sciences,, Altenbergerstrasse 69, Linz, 4040, Upper Austria, Austria

bRICAM, Austrian Academy of Sciences,, Altenbergerstrasse 69, Linz, 4040, Upper Austria, Austria

Abstract

We analyze the inversion of the photo-acoustic imaging modality using electromagnetic plasmonic nano-particles as contrast agents. We prove that the generated pressure, be- fore and after injecting the plasmonic nano-particles, measured at a single point, located away from the targeted inhomogeneity to image, and at a given band of incident frequen- cies is enough to reconstruct the (eventually complex valued) permittivity. Indeed, from these measurements, we define an indicator function which depends on the used incident frequency and the time of measurement. This indicator function has differentiating be- haviors in terms of both time and frequency. First, from the behavior in terms of time, we can estimate the arrival time of the pressure from which we can localize the injected nano-particle. Second, we show that this indicator function has maximum picks at in- cident frequencies close to the plasmonic resonances. This allows us to estimate these resonances from which we construct the permittivity.

To justify these properties, we derive the dominant electric field created by the injected nano-particle when the incident frequency is close to plasmonic resonances. This is done for the full Maxwell system. To this end, we use a natural spectral decomposition of the vector space (L²(D))³ based on the spectra of the Newtonian and the Magnetization operators. As another key argument, we provide the singularity analysis of the Green’s tensor of the Maxwell problem with varying permittivity. Such singularity is unusual if compared to the ones of the elliptic case (as the acoustic or elastic models). In addition, we show how the derived approximation of the electric fields propagates, as a source, in the induced pressure with the time.

Keywords: photo-acoustic imaging, plasmonic nano-particles, surface plasmon resonance, inverse problems, Maxwell system.

2000 MSC: 35R30, 35C20

∗Corresponding author.

Email addresses: [email protected](Ahcene Ghandriche ), [email protected](Mourad Sini)

1. Introduction and statement of the results 1.1. Introduction

The photo-acoustic experiment, in the general setting, applies to targets that are electrically conducting, in other words the imaginary part of the ’permittivity’ is quite pronounce, and it goes as follows. Exciting the target, with laser, or by sending an incident electric field, will create heat in its surrounding. This heat, in its turn, creates fluctuations, i.e. a pressure field, that propagates along the body to image. This pressure can be collected in an accessible part of the boundary of the target. The photo-acoustic imaging is to trace back the pressure and reconstruct the permittivity that created it.

In our settings, the source of the heat is given by the injected electromagnetic nano- particles. To describe the mathematical model behind this experiment, let us set E, T and pto be respectively the electric field, the heat temperature and the acoustic pressure.

Then, as described above, the photo-acoustic experiment is based on the following model coupling these three equations:











Curl Curl E−ω² ε µ E = 0, E :=E^s+Eⁱ, in R³, ρ0cp

∂T

∂t − ∇ ·κ∇T =ω Im (ε)|E|² δ0(t), in R³×R⁺, 1

c²

∂²p

∂t² −∆p=ρ₀ β₀ ∂²T

∂t² , inR³×R+,

whereρ₀ is the mass density, c_p the heat capacity, κis the heat conductivity,cis the wave speed andβ₀ the thermal expansion coefficient. To the last two equations, we supplement the homogeneous initial conditions T = p = ^∂p_∂t = 0, at t = 0 and the Silver-M¨uller radiation condition to E^s. Under the condition that the heat conductivity is relatively small, the model above reduces to the following one:

(1.1)







∂_t²p(x, t)−c²_s(x)∆

xp(x, t) = 0 in R³×R⁺, p(x,0) = ^{ω β}_c ⁰

p Im (ε)(x)|E|²(x), in R³,

∂_tp(x,0) = 0 in R³,

here c_s is the velocity of sound in the medium that we assume to be a uniform constant.

The constants β₀ and c_p are known and ω is an incident frequency. The source E is solution of the scattering problem

(1.2)

Curl Curl E−ω² ε µ E = 0, E :=E^s+Eⁱ, inR³, E^s(x) satisfies the Silver-M¨uller radiation conditions,

whereε=ϵ_p insideD,ε=ϵ₀(x) outsideDandϵ₀(x) =ϵ∞outside a bounded and smooth domain Ω (D⊂Ω being the injected nano-particle with permittivity ϵ_p and permeability µ). More details on the actual derivation of this model can be found in [39, 45] and more references therein. The permittivity ϵ₀(·) is variable and it is supposed to be smooth inside Ω. The needed smoothness will be discussed later.

We have two classes of such nano-particles: dielectric and plasmonic nano-particles.

The dielectric nano-particles enjoy the following features. They are highly localized as they are nano-scaled and they have high contrast permittivity. Under these scales, we can choose the incident frequency so that we excite the dielectric resonances which are related to the eigenvalues of the Newtonian operator. The main feature of the plasmonic nano-particles is that they enjoy negative values of the real part of their permittivity if we choose incident frequencies close to the plasmonic frequencies of the nano-particle. With such negative permittivity, we can excite the plasmonic resonances which are related to the eigenvalues of the Magnetization operator. To describe this, we use the Lorentz model where the permeabilityµis kept constant as the one of the homogeneous background while the permittivity has the form:

(1.3) ϵ_p =ϵ∞

1 + ω_p² ω₀²−ω²+iγω

where ω_p is the electric plasma frequency, ω₀ is the undamped frequency and γ is the electric damping frequency. We observe that if we choose the incident frequency ω so that ω² is larger than ω²₀, then the real part becomes negative. For such choices of the incident frequency, the nano-particle behaves as a plasmonic nano-particle.

The goal of the photo-acoustic imaging using nano-particles is to recover ϵ₀(·) in Ω from the measure of the pressure p(x, t), x ∈ ∂Ω and t ∈ (0, T) for large enough T. The decoupling of the original photo-acoustic mathematical model (1.1) into (1.1)-(1.2) suggests that we split the inversion into the following two steps.

1. Acoustic Inversion: Recover the source term Im (ε)(x)|u|²(x), x ∈ Ω, from the measure of the pressure p(x, t), x∈∂Ω and t∈(0, T).

2. Electromagnetic Inversion: Recover the permittivity ϵ₀(x), x∈Ω from Im (ε)(x)|u|²(x), x∈Ω.

The pressure is collected on the boundary of Ω in the following situations:

ˆ Before injecting any particle. The measured data is the pressure p(x, t), x ∈ ∂Ω and t ∈ (0, T) without injecting any particle. There is a large literature based on such data. Without being exhaustive, we cite the following references [3, 4, 6, 11, 12, 8, 9, 10, 17, 23, 31, 33, 34, 36, 37, 42, 44] devoted to such inversions. The general approach is that, using the Radon transform, one can recover the initial pressure, i.e. Im (ε)(x)|E|²(x),x∈Ω. The next step is to use these internal values to recover the permittivity ϵ₀(·).

ˆ After injecting nano-particles. The measured data is the pressure p(x, t), x ∈ ∂Ω and t ∈ (0, T) after injecting a nano-particle. The first work in this direction is [45] where plasmonic nano-particles are used and an optimization method was proposed to invert the electric energy fields. There, the 2D-model is stated and the magnetic field was used. Assuming the initial pressure to be already given, via one of the inversion methods as the Radon transform for instance, the authors propose a reconstruction method to recover the permittivity from the modulus of the electric (or the magnetic) field given in and around the particle D. For this, they use the contrasting behavior of the magnetic field across the interface of the particle.

ˆ Before and after injecting nano-particles. The measured data is the pressure p(x, t), x∈ ∂Ω and t ∈ (0, T) before and then after injecting the nano-particle. In [29, 28], we considered the 2D-model using dielectric nano-particles. There, we did not split the problem into two steps. Rather, we derived direct formulas linking the measured pressure collected only on asingle pointxon the accessible surface, to the internal values of the modulus of the electric field. In addition, using dimers (two close nano-particles), we showed that we can reconstruct, not only the electric field, but also the values of the (real part of the ) Green’s functions on the centers of the dimer’s nano-particles. From this Green function, we recover the permittivity. The main argument there is that under critical scales, on the size and the high values of the permittivity, we can choose the incident frequency so that we excite the dielectric resonances which are related to the eigenvalues of the Newtonian operator. Here, we propose to use plasmonic nano-particles. Measuring the induced pressure before and after injecting such a nano-particle, on a single point of ∂Ω but a band of frequencies, and taking their difference, we show that the generated curve has picks on incidence frequencies close to singular frequencies related to the eigenvalues of the Magnetization operator (that are the plasmonic resonances). With such behavior, we can construct those resonances. From these resonances, we extract the values of the permittivity. More details are given in section 1.3.

1.2. Statement of the results

Let Ω be aC²-smooth and bounded domain. The nano-particleDis taken of the form D :=a B+z where z models its location and a its relative radius withB as C²-smooth domain of maximum radius 1. For later use, we introduce the integral operators of the volume potential N^k(·) and the Magnetization potential ∇M^k(·), both acting on vector fields:

(1.4) N^k(f)(x) :=

Z

B

Φ_k(x, y)f(y)dy and ∇M^k(f)(x) :=∇ Z

B

∇yΦ_k(x, y)·f(y)dy,

where Φ_k(x, y) := ^e_4π|x−y|^ik|x−y| is the fundamental solution for Helmholtz equation in the entire space. Particularly, for k = 0 we obtain:

(1.5)

N(f)(x) :=

Z

B

1

4π|x−y|f(y)dy and ∇M(f)(x) :=∇ Z

B

∇y

1 4π|x−y|

·f(y)dy.

We recall the decomposition

(L²(B))³ =H⁰(div = 0)(B)⊕H⁰(Curl = 0)(B)⊕ ∇Harmonic(B)

where H⁰(div = 0)(B) := {u∈(L²(B))³,div (u) = 0;ν·u= 0 on ∂B}, H⁰(Curl = 0) :=

{u∈(L²(B))³, Curl(u) = 0;ν×u= 0 on∂B}and ∇Harmonic(B) :={u=∇ϕ, ∆ϕ = 0 in B}.

We can show, see later, thatN|_{H0(div =0)} andN|_H0(Curl=0) generate complete orthonor- mal bases

λ⁽¹⁾n , e⁽¹⁾n

n∈N

and

λ⁽²⁾n , e⁽²⁾n

n∈N

of H⁰(div = 0) and H⁰(Curl= 0) respec- tively. In addition, it is known that ∇M :∇Harmonic → ∇Harmonic has a complete basis

λ⁽³⁾n , e⁽³⁾n

n∈N

.

The permittivityε(·) is defined as

(1.6) ε(x) :=







ϵ∞ in R³\Ω, ϵ₀(x) in Ω\D, ϵ_p in D, where ϵ_p is given by (1.3).

Related to this, we set the index of refractionn, in R³, given by¹ n:=

(√

ϵ_pµ in D

n₀ in R³\D and n₀ :=

(pϵ₀(·)µ in Ω

√ϵ∞µ inR³\Ω. (1.7)

We assumeϵ₀(·) to be of class C¹. Let z ∈Ω and define f_n(ω, γ) :=ϵ₀(z)−(ϵ₀(z)−ϵ_p)λ³_n. We show that in the square

ω₀;

q

ω_p²+ω²₀ =:ω_max

× 0; ω_max

Im (ϵ0(·)) Re (ϵ₀(·)) L^∞(Ω)

=:γ_max

! ,

the dispersion equation f_n(ω, γ) = 0 has one and only one solution. For any n₀ fixed, we set (ω_n0, γ_n0) to be the corresponding solution for n=n₀.

From now on, we use the notation u₀ :=E, respectively u₁ :=E, for the solution of (1.2) before, respectively after, injecting the nano-particle D inside Ω.

1For z ∈ C, given by z = r e^iϕ with −π < ϕ ≤ π, the principal square root of z is defined to be:

√z=√ r e^iϕ2.

Theorem 1.1. We assume Ω and B (and hence D) to be of class C². In addition, we assume ϵ₀(·) to be of class C¹ and satisfies the conditions²

Reϵ₀(·)> ϵ_∞ and ∥ϵ₀(·)−ϵ_∞∥_L^∞_(Ω) ≤C

where the positive constant C := C(Ω, ω², µ, ϵ∞) is given by (3.19) and depends on the used frequency ω², the permeability µ, the permittivity ϵ∞ and the domain Ω through the mapping property of the Newtonian operator.

Let the used incident and damping frequencies (ω, γ) be such that (1.8) ω²−ω_n²₀ ∼a^h and γ−γ_n0 ∼a^h for h∈(0,1), where a >0 is the radius of the nano-particle D.

1) We have the following approximation of the electric field³

(1.9)

Z

D

|u₁|²(x)dx=a³

|ϵ₀(z)|²

⟨u₀(z);R

Be⁽³⁾n0(x)dx⟩

2

ϵ0(z)−(ϵ0(z)−ϵp) λ⁽³⁾n0

2 +O amin(3,4−3h) .

2) Let x ∈ ∂Ω and s ≥ diam(D) + dist(x, D). We have the following approximation of the average pressure:

(1.10)

p^⋆(x, s)−p^⋆₀(x, s) = a³

4πIm (ϵ_p)

|ϵ₀(z)|²

⟨u₀(z);R

Be⁽³⁾n0(x)dx⟩

2

ϵ₀(z)−(ϵ₀(z)−ϵ_p)λ⁽³⁾n0

2 +O amin(3−h,4−3h) ,

where

p^⋆(x, s) :=

Z s 0

r Z r

0

p(x, t)dt dr and p^⋆₀(x, s) :=

Z s 0

r Z r

0

p₀(x, t)dt dr

with p₀(·,·) being the pressure generated by the medium in the absence of the nano- particle. In (1.9) and (1.10) the expression

⟨u₀(z);R

Be⁽³⁾n0(x)dx⟩

2

should be understood as P

m

⟨u₀(z);R

Be⁽³⁾n0,m(x)dx⟩

2

, where e⁽³⁾n0;m are such that ∇ e⁽³⁾n0,m

=λ⁽³⁾n0 e⁽³⁾n0,m.

2The first condition is a natural one in applications. The second condition is needed to derive and analyze (the singularity of the) the Green’s kernel for the inhomogeneous Maxwell system.

3For g(·) positive function, we say that f(δ) =O(g(δ)), asδ →0, if there exist positive numbersβ and C such that for allδ with 0< δ < β,|f(δ)| ≤C g(δ). In other wordsf does not grow much faster thang.

To justify these results, we need to derive the dominating fields in both the acoustic and the electromagnetic models that constitute the photo-acoustic model. The most difficult issue is in deriving the dominating electric field generated by the plasmonic resonances.

These are related to the eigenvalues of the Magnetization operator, at zero wave number, restricted to the sub-space of grad-harmonic functions. To do this, we first use the natu- ral decomposition of the L²-space as L² = H⁰(div = 0)⊕H⁰(Curl = 0)⊕ ∇Harmonic to which we correspond a spectral decomposition given by the eigenvalues-eigenfunctions of the vector Newtonian operator restricted to H0(div = 0) and H0(Curl = 0) and the above mentioned eigenvalues-eigenfunctions of the Magnetization operator. As the elec- tromagnetic background is inhomogeneous, unlike the elliptic models, as in acoustics or elasticity for instance, the corresponding Green’s kernel has unusual singularities. We derive two main decompositions of this kernel in terms of the singularities. The first one is valid inR³ as a sum of the explicit kernel, of the Maxwell system with constant electro- magnetic parameters, and a kernel⁴ having L

3

2− integrable singularity. This allows us to define and justify the equivalent Lippmann-Schwinger equation. The second one is valid locally, near any fixed source point, as a sum of three kernels. The first one is the kernel of the grad-harmonic operator, i.e. the kernel of the Magnetization operator, the second is of the form ∇K, with K having an L³⁻ integrable singularity and the third one with anL³⁻ integrable singularity as well. This last decomposition is key in deriving the dom- inating electric field generated by the plasmonic resonances. These decompositions of the Green’s kernel, see for instance (3.52), are more precise and general than the ones derived in the case of two half spaces with two different permittivities constants, see for instance [18] and [32]. Let us finally mention the work [5], where the spectral decomposition of the Magnetization operator on ∇Harmonic is used for a different purpose and in a homoge- neous medium. In our work, we need the spectral decomposition of the whole L²-space using the spectrum of the Newtonian operator on H0(div = 0) and H0(Curl = 0) and the one of the Magnetization operator on∇Harmonic. In addition, here we deal with an inhomogeneous background which makes the analysis more involved but useful for itself.

1.3. Inversion of the photo-acoustic imaging modality using plasmonic contrast agents Here, we discuss how to use the approximation formulas we derived in the previous section to the actual inversion procedure for the photo-acoustic modality using plasmonic nano-particles as contrast agents.

1. Let x ∈ ∂Ω be fixed. Let also (ω, γ) be any couple of incident and damping fre- quencies.

(a) If s < dist(x, D), then p^⋆(x, s) − p^⋆₀(x, s) = O a^3−h

for any (ω, γ). This property, which is related to the finiteness of the speed of propagation, can be shown by combining (2.37) and Lemma2.3.

4Fors >1 we denote byL^s−the space of functions belonging toL^s−δfor everyδsuch that 0< δ < s−1.

(b) If s ≥ diam(D) + dist(x, D), then, under the condition of the existence of n₀ ∈N such that R

Be⁽³⁾n0(x)dx̸= 0, we have p^⋆(x, s)−p^⋆₀(x, s)∼a^3−2h, for any (ω, γ) close to (ω_n0, γ_n0) as in (1.8). This comes from (1.10).

From these formulas, we can estimate dist(x, D) with an error of the order of diam(D) ∼ a. Therefore, measuring p^⋆(x, s) −p^⋆₀(x, s) for four different points x₁, x₂, x₃, x₄ on ∂Ω, such that two of these four points are contained in the same line, we can localize the injected nano-particle with an error of the order a. More details about the localization of the nano-particle Dcan be found in Appendix 5.4.

2. Let x∈∂Ω ands ≥diam(D) +dist(x, D) be fixed, we define Iz(ω, γ) := |p^⋆(x, s, ω)−p^⋆₀(x, s, ω)|

on the square

ω₀; q

ω_p²+ω₀² :=ω_max

×

0; ω_max

Im (ϵ0(·)) Re (ϵ0(·))

L^∞(Ω):=γ_max

.

According to (1.10), this functional has a sequence of picks (ω_n, γ_n), n = 1,2, ...

Observe that the index n is related to one of the eigenvalue λ⁽³⁾n of ∇M. From Lemma 5.7, we know that

(1.11) λ⁽³⁾_n < λ⁽³⁾_m ⇒ω²_n< ω²_m.

Therefore from the picks (ωn, γn) of the functional Iz(ω, γ) we can choose anyone of them, say (ωn0, γn0). From (1.11), to ωn0 we correspond a unique λ⁽³⁾n0, via the ordering of λ⁽³⁾n ’s. Fromf_n0(ω_n0, γ_n0) = 0, we obtain

ϵ₀(z) = − ϵ_pλ⁽³⁾n0

1−λ⁽³⁾n0

.

Observe that the validity of the imaging procedure works for nano-particles for which R

Be⁽³⁾n (x)dx̸= 0 for some n’s, see the condition mentioned in (1b). This condition can be clarified for particular shapes. For nano-particles B of ellipsoidal shape, with semi-axes given by r₁, r₂ and r₃ we use the fact that:

(1.12)

N(1)(x) = r1r2r3

4

Z ∞ 0

1−

3

X

j=1

x²_j s+r²_j

! 1

p(s+r²₁) (s+r₂²) (s+r₃²)ds, x∈B,

see for instance Theorem 1.1 of [43]. Therefore, by straightforward computations, using the relation ∇M(I) =−∇div (N(I)) =−∇div (N(1)I), we derive:

(1.13) ∇M(I) (x) = r₁r₂r₃ 2





I1(r₁, r₂, r₃) 0 0 0 I2(r₁, r₂, r₃) 0 0 0 I3(r₁, r₂, r₃)



, x∈B, where, for j = 1,2,3, we have:

(1.14) Ij(r₁, r₂, r₃) :=

Z +∞

0

1 s+r²_j

1

p(s+r²₁) (s+r²₂) (s+r²₃)ds.

Using the fact that λ⁽³⁾n e⁽³⁾n =∇M e⁽³⁾n

we get, λ⁽³⁾_n I·

Z

B

e⁽³⁾_n (x)dx = Z

B

λ⁽³⁾_n e⁽³⁾_n (x)dx= Z

B

∇M e⁽³⁾_n

(x)dx= Z

B

I· ∇M e⁽³⁾_n

(x)dx, which, by using the self-adjointness of the Magnetization operator, becomes

(1.15) λ⁽³⁾_n I· Z

B

e⁽³⁾_n (x)dx= Z

B

∇M(I) (x)·e⁽³⁾_n (x)dx.

Now, thanks to the constancy of ∇M(I) inside B, see (1.13), we deduce that:



λ⁽³⁾_n I− r₁r₂r₃ 2





I1(r₁, r₂, r₃) 0 0 0 I2(r₁, r₂, r₃) 0 0 0 I3(r₁, r₂, r₃)







· Z

B

e⁽³⁾_n (x)dx=



 0 0 0



. As the matrix on the left hand side is a diagonal one, we deduce that we can have at most three eigenvalues, λ⁽³⁾n , for which the corresponding eigenfunctions might have the property R

Be⁽³⁾n (x)dx ̸= 0. The computation of R

Be⁽³⁾n (x)dx in the case of ellipsoidal shape is a difficult task. We restrict our computations to the particular case of a unit ball, which corresponds to taker₁ =r₂ =r₃ = 1. By straightforward computations, using the definition of∇M(I), see (1.13), and the formula of Ij(r1, r2, r3), see (1.14), we obtain

(1.16) ∇M(I)(x) = 1

3I, x∈B.

Therefore, ¹₃ is the only eigenvalue for which the corresponding eigenfunctions might have non-zero average. These eigenfunctions are of the form∇SL(Re(Y₁^m)) or∇SL(Im(Y₁^m)), where Y₁^m, |m| ≤1, are the spherical harmonics. In addition, we show that

(1.17) Z

B

∇SL(Re(Y₁^m)) (x)dx= χ

m=−1

√6π 9



 1 0 0



+ χ

m=0

2 3

rπ 3



 0 0 1



− χ

m=1

√6π 9



 1 0 0





and (1.18)

Z

B

∇SL(Im(Y₁^m)) (x)dx=− χ

m=−1

√6π 9



 0 1 0



+ χ

m=0



 0 0 0



− χ

m=1

√6π 9



 0 1 0



, whereχis the discrete characteristic function. More details about the derivation of (1.17) and (1.18) can be found in Appendix 5.3.

Combining (1.16),(1.17),(1.18) and (1.15) we end up with λ⁽³⁾₁ = ¹₃, which is the first eigenvalue of the Magnetization operator ∇M(·) in the unit ball, and the corresponding eigenfunctions have non-zero average. In this case, by a straightforward computations we get:

⟨u0(z);

Z

B(0,1)

e⁽³⁾₁ (x)dx⟩

2

= X

|m|≤1

⟨u0(z);

Z

B(0,1)

∇SL(Re(Y₁^m)) (x)dx⟩

2

+ X

|m|≤1

⟨u₀(z);

Z

B(0,1)

∇SL(Im(Y₁^m)) (x)dx⟩

2

= 4π

27 |u₀(z)|². Hence, in the imaging procedure described above, using a nanoparticle with a spherical shape, the functional

I_z(ω, γ) := |p^⋆(x, s, ω)−p^⋆₀(x, s, ω)|

has one and only one pick in the square

ω₀; q

ω_p²+ω²₀ :=ω_max

× 0; ω_max

Im (ϵ₀(·)) Re (ϵ₀(·)) L^∞(Ω)

:=γ_max

! .

Remark 1.2. It is worth noticing that the results above can also be derived using the Drude model for the permittivity, see for instance [21] formula (1.5), instead of the Lorentz model.

Remark 1.3. In our analysis, as we have seen, it is mandatory to vary both the incident frequenciesωand the damping frequenciesγ. This is can turn out to be expensive from the point of view of applications as it would mean that one should change the nano-particles to change γ. In the following, we give two ways to overcome this eventual issue.

1. In the case where Im (ϵ₀(z)) is very small, or mathematically zero, then we can take the damping frequency γ small as well but fixed. In this case, we only need to vary the incident frequency. Observe that the traditional photo-acoustic experiment applies only to electrically highly conducting tissues. However, the photo-acoustic experiment based on using contrast agents can deal with non-conductive tissues as well (as for benign or early stage tumors).

2. Instead of varying both the incident and damping frequencies, we allow the frequen- cies ω to be in the complex plan and fix the damping frequency. In this case, we can derive similar reconstruction formulas.

More details can be found in Remark 5.8.

The remaining part of the manuscript is divided as follows. In Section 2, we give the proof of Theorem 1.1 postponing the construction of the Green’s tensor, the invertibility of the Lippmann-Schwinger equation and certain a priori estimates of the electric fields to the next sections. In Section 3, we construct the Green’s tensor for our Maxwell model and provide its singularity analysis. These properties are used then to derive, and give sense to, the Lippmann-Schwinger equation. In Section 4, we prove the a priori estimates used in the proof of Theorem 1.1. In Appendix 5.1, we provide the spectral decomposition of the vector space (L²(D))³ based on the eigenvalues of the Newtonian and Magnetization operators, in Appendix 5.2 we analyze the dispersion equations f_n(ω, γ) = 0 used also in the proof of Theorem 1.1, in Appendix 5.3 we provide detailed computations of the integral of the first eigenfunction of Magnetization operator and in Appendix 5.4 we justify the localization of the nano-particle D inside Ω.

2. Proof of Theorem 1.1

2.1. Approximation of the Lippmann Schwinger equation and proof of (1.9)

As shown in Section 3.2 the solution, in distributional sense, of (1.2) can be written as solution of the following Lippmann-Schwinger equation

(2.1) u₁(x) +ω² Z

D

G_k(x, y)·u₁(y) (n²₀(y)−n²(y))dy=u₀(x), x∈R³,

where G_k(·,·) is the Green kernel for Maxwell’s equation for the inhomogeneous back- ground, defined as solution of:

∇y × ∇

y ×G_k(x, y)−ω²n²₀(y)G_k(x, y) = δ

x(y)I, x, y∈R³,

such that each column of G_k(·,·) satisfies the outgoing radiation condition lim

|x|→+∞ |x|

∇y ×G_k(x, y)× x

|x| −i k G_k(x, y)

= 0.

Recalling the definition of the index of refraction n(·), see (1.7), we get (2.2) u₁(x) +ω²µ

Z

D

G_k(x, y)·u₁(y) (ϵ₀(y)−ϵ_p)dy =u₀(x), x∈R³,

The formal representation of the convolution part is justified in Section 3.2 as well. The following theorem, on the singularity analysis of the Green’s function, is of importance to reduce the complexities of the integral equation (2.2), and then invert it.

Theorem 2.1. The Green kernel G_k(·,·) admit the following decomposition:

(2.3) G_k(x, z) = Υ(x, z) + Γ(x, z), x̸=z, where Υ(·,·) is the kernel defined by

(2.4) Υ(x, z) := 1

ω²µ ϵ₀(z)∇

x div

x (Φ₀(x, z)I), x̸=z, and the remainder part Γ(·,·), for all element x near z, is given by:

(2.5) Γ(x, z) := −1

ω²µ(ϵ₀(z))² ∇

x ∇

xM(Φ0(·, z)∇ϵ0(z)) (x) +W4(x, z), x̸=z,

where, for an arbitrarily and sufficiently small positive δ, the term W4(·, z) is an element in L

3(3−2δ) (3+2δ) (D).

Proof. The proof and details about the decomposition of the kernel Gk(·,·) are given in Section 3.1.

We restrict the study of the equation (2.2) to the domain D and we use the decom- position (2.3) to rewrite it as:

(2.6) u₁(x) +ω²µ Z

D

Υ(x, y)·u₁(y) (ϵ₀(y)−ϵ_p)dy=u₀(x) +Err_Γ(x), where Err_Γ(x) is the vector field given by

(2.7) Err_Γ(x) :=−ω²µ Z

D

Γ(x, y)·u₁(y) (ϵ₀(y)−ϵ_p)dy, x∈D.

Now, using the expression (2.4), we reformulate (2.6) as u₁(x)− ∇

x

Z

D

∇yΦ₀(x, y)·u₁(y)(ϵ₀(y)−ϵ_p)

ϵ₀(y) dy =u₀(x) +Err_Γ(x).

We set

(2.8) η(·) := (ϵ₀(·)−ϵ_p)

ϵ₀(·)

and use the definition of the Magnetization operator, see (1.5), to rewrite the previous equation as:

u₁(x)− ∇M(u₁η) (x) = u₀(x) +Err_Γ(x),

and then, by Taylor expansion for the function η(·) near the center z, we get (2.9) u₁(x)−η(z)∇M(u₁)(x) =u₀(x) +Err₀(x) +Err_Γ(x), where

(2.10) Err₀(x) :=∇M

u₁(·) Z 1

0

∇η(z+t(· −z))·(· −z)dt

(x).

Set W(·) to be the scattering matrix defined by

(2.11) W(·) = h

I− η(z)∇Mi−1

(I) (·).

Then, successively, taking the inverse of [I− η(z)∇M], on both sides of (2.9), integrating over D the obtained equation and using the definition of the matrixW(·), we get

(2.12) Z

D

u₁(x)dx = Z

D

W(x)·[u₀(x) +Err₀(x) +Err_Γ(x)] dx= Z

D

W(x)dx·u₀(z) +Err₁, where

Err₁ :=

Z

D

W(x)· Z 1

0

∇u₀(z+t(x−z))·(x−z)dt+Err₀(x) +Err_Γ(x)

dx.

Next, we estimate Err1. For this, we split it as

(2.13) Err₁ :=S₁+S₂+S₃,

then, we define each term and estimate it. More precisely, we have:

1. Estimation of:

S₁ :=

Z

D

W(x)· Z 1

0

∇u₀(z+t(x−z))·(x−z)dt dx

|S1| ≤ ∥W∥

L²(D)

Z 1 0

∇u0(z+t(· −z))·(· −z)dt L²(D)

= O

a⁵² ∥W∥

L²(D)

. (2.14)

2. Estimation of:

S₂ :=

Z

D

W(x)·Err₀(x)dx S₂ ^(2.10)=

Z

D

W(x)· ∇M

u₁(·) Z 1

0

∇η(z+t(· −z))·(· −z)dt

(x)dx.

We apply the Cauchy-Schwartz inequality to obtain:

|S2| ≲ ∥W∥

L²(D)

∇M

u1(·) Z 1

0

∇η(z+t(· −z))·(· −z)dt

L²(D) (5.14)

≤ ∥W∥

L²(D)

u₁(·) Z 1

0

∇η(z+t(· −z))·(· −z)dt L²(D)

. Then,

(2.15) S2 =O

∥W∥

L²(D) a ∥u1∥

L²(D)

. 3. Estimation of:

S₃ :=

Z

D

W(x)·Err_Γ(x)dx^(2.7)≃ Z

D

W(x)· Z

D

Γ(x, y)·u₁(y) (ϵ₀(y)−ϵ_p)dy dx.

With the help of (2.5) we rewrite the previous formula as S₃ ≃

Z

D

W(x)· Z

D

∇x

h∇

xM(Φ₀(·, y)∇ϵ₀(y))i

(x)·u₁(y)(ϵ₀(y)−ϵ_p) ϵ²₀(y) dy dx +

Z

D

W(x)· Z

D

W₄(x, y)·u₁(y) (ϵ₀(y)−ϵ_p)dy dx.

We split the previous formula as S₃ = S_3,1+S_3,2, we define and we estimate each term.

(a) Estimation of:

S_3,1 :=

Z

D

W(x)· Z

D

∇x

h∇

xM(Φ₀(·, y)∇ϵ₀(y))i

(x)·u₁(y)(ϵ₀(y)−ϵ_p) ϵ²₀(y) dy dx

= Z

D

W(x)· Z

D

∇x∇

x

Z

D

∇tΦ₀(t, x)· ∇ϵ₀(y)Φ₀(t, y)dt·u₁(y)(ϵ₀(y)−ϵ_p) ϵ²₀(y) dydx

= Z

D

W(x)· ∇

x∇

x · Z

D

Z

D

∇tΦ0(t, x)· ∇ϵ0(y)Φ0(t, y)u1(y)(ϵ₀(y)−ϵ_p)

ϵ²₀(y) dtdydx

= Z

D

W(x)· ∇

x∇

x · Z

D

Z

D

Φ₀(t, y)∇ϵ₀(y)⊗u₁(y)(ϵ₀(y)−ϵ_p) ϵ²₀(y) dy· ∇

tΦ₀(t, x)dtdx.

For shortness, for every y∈D, we set

(2.16) u^⋆₁(y) :=∇ϵ₀(y)⊗u₁(y)(ϵ₀(y)−ϵ_p)

ϵ²₀(y) ∈L²(D), hence,

S_3,1 = Z

D

W(x)· ∇

x ∇

x · Z

D

Z

D

Φ₀(t, y)u^⋆₁(y)dy · ∇

tΦ₀(t, x)dt dx

= Z

D

W(x)· ∇

x ∇

x · Z

D

N(u^⋆₁) (t)· ∇

tΦ₀(t, x)dt dx.

By an integration by parts, we rewrite the last formula as:

S3,1 = − Z

D

W(x)· ∇ ∇ ·N(divN(u^⋆₁)) (x)dx +

Z

D

W(x)· ∇ ∇ ·SL(ν·N(u^⋆₁)) (x)dx.

Again, we splitS_3,1 asS_3,1 =S_3,1,1+S_3,1,2, we define each term and we estimate it.

As u^⋆₁ ∈L²(D) the function∇ ∇ ·N(divN(u^⋆₁))∈L²(D). This implies,

|S_3,1,1| :=

Z

D

W(x)· ∇ ∇ ·N(divN(u^⋆₁)) (x)dx

≤ ∥W∥

L²(D) ∥∇ ∇ ·N(divN(u^⋆₁))∥

L²(D),

and, from Calderon-Zygmund inequality, see [30], page 242, we reduce the previous inequality to

|S_3,1,1| ≤ ∥W∥

L²(D) ∥divN(u^⋆₁)∥

L²(D).

Just as divergence operator div and gradient operator ∇ are both differential operator of order one, we use (5.13) to finish the estimation ofS_3,1,1. We have,

(2.17) S_3,1,1 =O

a ∥u^⋆₁∥

L²(D) ∥W∥

L²(D)

. Also, we have,

|S_3,1,2| :=

Z

D

W(x)· ∇ ∇ ·SL(ν·N(u^⋆₁)) (x)dx

≤ ∥W∥

L²(D) ∥∇ ∇ ·SL(ν·N(u^⋆₁))∥

L²(D), and, after scaling, we obtain:

|S_3,1,2| ≤ ∥W∥

L²(D) a³² a²

∇ ∇ ·SL

ν·N ue^⋆₁

L²(B)

.

Using the continuity of the operator ∇ ∇ · SL : H

1

2(∂B) → L²(B), see for instance [35], corollary 6.14, page 210, we reduce the previous equation to

|S3,1,2|≲∥W∥

L²(D) a³² a² ν·N

ue^⋆₁

H^1/2(∂B). Now, from the continuity of the trace operator we deduce that

|S_3,1,2|≲∥W∥

L²(D) a³² a² N

ue^⋆₁ H¹(B).

In addition, using the continuity of the Newtonian operator and scaling back to end up with

(2.18) S3,1,2 =O

∥W∥

L²(D) a² ∥u^⋆₁∥

L²(D)

. By gathering (2.17) with (2.18) we deduce that:

S_3,1 =O

∥W∥

L²(D) a ∥u^⋆₁∥

L²(D)

.

Viewing the definition ofu^⋆₁(·), see for instance (2.16), and using the smoothness of η(·)

ϵ₀(·)∇ϵ₀(·) we deduce that

(2.19) S_3,1 =O

∥W∥

L²(D) a ∥u₁∥

L²(D)

. (b) Estimation of:

S_3,2 :=

Z

D

W(x)· Z

D

W₄(x, y)·u₁(y) (ϵ₀(y)−ϵ_p)dy dx.

Then, by Holder inequality, we obtain:

|S_3,2| ≤ ∥W∥

L²(D)

Z

D

W₄(·, y)·u₁(y) (ϵ₀(y)−ϵ_p)dy L²(D)

≤ ∥W∥

L²(D) ∥u₁∥

L²(D)

Z

D

Z

D

|W₄(x, y) (ϵ₀(y)−ϵ_p)|²dy dx ¹₂

≤ ∥W∥

L²(D) ∥u₁∥

L²(D) Sup

D

|ϵ₀(·)−ϵ_p| Z

D

Z

D

|W₄(x, y)|²dy dx ¹₂

. We know that (ϵ₀(·)−ϵ_p)∼1, with respect to the size a, we deduce that (2.20) |S3,2|≲∥W∥

L²(D) ∥u1∥

L²(D)

Z

D

Z

D

|W4(x, y)|²dydx ¹₂

.

Now, we recall from (2.5) thatW₄(·, y)∈L

3(3−2δ)

(3+2δ)(D) and, then, we approximate its singularity as

W₄(·, y)≃ |· −y|⁻(3+2δ)(3−δ) 3(3−2δ)

and, then, by scaling the double integral in (2.20) we get

(2.21) S_3,2 =O

∥W∥

L²(D) ∥u₁∥

L²(D) a^2δ

2−21δ+18 3(3−δ)

. Finally, by summing (2.19) and (2.21) we get:

(2.22) S₃ =O

∥W∥

L²(D) ∥u₁∥

L²(D) a .

Using (2.14),(2.15),(2.22) and the definition of Err₁, see (2.13), we get Err₁ =O

a⁵² ∥W∥

L²(D)+a ∥u₁∥

L²(D) ∥W∥

L²(D)

. Going back to (2.12) and plugging the expression ofErr₁ to obtain (2.23)

Z

D

u₁(x)dx= Z

D

W(x)dx·u₀(z) +O

a⁵² ∥W∥

L²(D)+a ∥u₁∥

L²(D) ∥W∥

L²(D)

. The next proposition show the estimates of the terms appearing on the right hand-side.

Proposition 2.2. We have:

(2.24) ∥u₁∥

L²(D)≤a^−h ∥u₀∥

L²(D), ∥W∥

L²(D)=O a^32−h and

Z

D

W(x)dx= a³ϵ₀(z)

ϵ₀(z)−(ϵ₀(z)−ϵ_p) λ⁽³⁾n0

Z

B

e⁽³⁾_n0(x)dx ⊗ Z

B

e⁽³⁾_n0(x)dx+O a³ . Proof. See Subsection 4.1 and Subsection 4.3.

Thanks to Proposition 2.2, the equation (2.23) takes the following form:

Z

D

u₁(x)dx= a³ ϵ₀(z)⟨u₀(z);R

Be⁽³⁾n0(x)dx⟩

ϵ₀(z)−(ϵ₀(z)−ϵ_p) λ⁽³⁾n0

Z

B

e⁽³⁾_n0(x)dx+O amin(3,4−2h) .

Consequently,

Z

D

u₁(x)dx

2

=

a⁶ |ϵ₀(z)|²

⟨u₀(z);R

Be⁽³⁾n0(x)dx⟩

2

ϵ₀(z)−(ϵ₀(z)−ϵ_p) λ⁽³⁾n0

2

Z

B

e⁽³⁾_n

0(x)dx

2

+O amin(6−h,7−3h) .

In the sequel, we derive a relation between the given data ∥u₁∥

L²(D) and R

Du₁(x)dx . For this, we recall from (2.9) that, in the domain D, we have

u₁ = [I−η(z)∇M]⁻¹[u₀+Err₀+Err_Γ], and after scaling to the domain B we obtain

˜

u₁ = [I−η(z)∇M]⁻¹h

˜

u₀+Errg₀+Errg_Γi . Recalling the decomposition of L²(B), see (5.1),

L²(B) =H0(div = 0)⊕^⊥ H0(Curl = 0)⊕ ∇^⊥ Harmonic, we project the previous equation into three subspaces.

1. Taking the inner product with respect toe⁽¹⁾n (·):

(2.25) ⟨u˜₁;e⁽¹⁾_n ⟩=⟨e⁽¹⁾_n ; [I−η(z)∇M]⁻¹h

˜

u₀+Errg₀+Errg_Γi

⟩, using the self-adjointness of ∇M(·) and the fact that ∇M

e⁽¹⁾n

= 0, see Lemma 5.5, we deduce that [I−η(z)∇M]⁻¹

e⁽¹⁾n

=e⁽¹⁾n and we reduce the equation (2.25) to:

⟨˜u₁;e⁽¹⁾_n ⟩=⟨e⁽¹⁾_n ;h

˜

u₀+Errg₀+Errg_Γi

⟩=⟨e⁽¹⁾_n ; ˜u₀⟩+Err_3,n, where Err_3,n is the term given by:

Err_3,n := ⟨e⁽¹⁾_n ;Errg₀⟩+⟨e⁽¹⁾_n ;Errg_Γ⟩

(2.10)

= −ω²µ a⟨e⁽¹⁾_n ;∇M

˜ u₁(·)

Z 1 0

∇η(zf +ta·)·(·)dt

⟩+⟨e⁽¹⁾_n ;Errg_Γ⟩.

We use the fact that ∇M e⁽¹⁾n

= 0 to reduce the previous expression to Err_3,n =

⟨e⁽¹⁾n ;Errg_Γ⟩.Hence,

⟨e⁽¹⁾_n ; ˜u₁⟩=⟨e⁽¹⁾_n ; ˜u₀⟩+⟨e⁽¹⁾_n ;Errg_Γ⟩.

By taking the square modulus of the previous equality and then the series with respect to n, we obtain:

X

n

⟨˜u₁;e⁽¹⁾_n ⟩

2 = X

n

⟨˜u₀;e⁽¹⁾_n ⟩

2+X

n

⟨Errg_Γ;e⁽¹⁾_n ⟩

2

+ O



 X

n

⟨˜u₀;e⁽¹⁾_n ⟩

2

!¹₂ X

n

⟨Errg_Γ;e⁽¹⁾_n ⟩

2!¹₂

,

and, using (4.6),(4.10) and (2.24), we get:

(2.26) X

n

⟨˜u₁;e⁽¹⁾_n ⟩

2 =X

n

⟨˜u₀;e⁽¹⁾_n ⟩

2+O

a^{(9−7δ−2δ2)(3−2δ) −h}

=O a^min

2;^{(9−7δ−2δ2)}_(3−2δ) −h

!

. 2. Taking the inner product with respect toe⁽²⁾n (·).

⟨u˜₁;e⁽²⁾_n ⟩=⟨e⁽²⁾_n ; [I−η(z)∇M]⁻¹h

˜

u₀+Errg₀+Errg_Γi

⟩, since ∇M

e⁽²⁾n

=e⁽²⁾n , see Lemma (5.5), then after taking the adjoint operator of [I−η(z)∇M] the previous equation will be reduced to

⟨˜u₁;e⁽²⁾_n ⟩ =

ϵ₀(z) h

⟨e⁽²⁾n ; ˜u₀⟩+⟨e⁽²⁾n ;Errg₀⟩+⟨e⁽²⁾n ;Errg_Γ⟩i ϵ_p

= ϵ₀(z)⟨e⁽²⁾n ; ˜u₀⟩ ϵp

+Err_4,n, (2.27)

where, obviously, the term Err_4,n is given by Err_4,n :=

ϵ₀(z) h

⟨e⁽²⁾n ;Errg₀⟩+⟨e⁽²⁾n ;Errg_Γ⟩i

ϵ_p .

Now, as ϵ₀(z) ϵp

∼1, with respect to the sizea, we approximate Err_4,n by:

Err4,n ≃ ⟨e⁽²⁾_n ;Errg0⟩+⟨e⁽²⁾_n ;ErrgΓ⟩.

Using the definition ofErr₀, see for instance (2.10), and the fact that ∇M e⁽²⁾n

= e⁽²⁾n to get:

Err_4,n ≃a⟨e⁽²⁾_n ; ˜u₁(·) Z 1

0

∇gη(z+ta·)·(·)dt⟩+⟨e⁽²⁾_n ;Errg_Γ⟩.

Consequently, X

n

|Err_4,n|² ≲a² ∥˜u₁∥²

L²(B)+X

n

⟨e⁽²⁾_n ;Errg_Γ⟩

2 (4.14)

= O

a² ∥˜u₁∥²

L²(B)

, and, using the a priori estimation (2.24), we obtain:

(2.28) X

n

|Err_4,n|² =O a^2−2h . Hence, using (2.28) in (2.27) we obtain:

(2.29) X

n

⟨˜u₁;e⁽²⁾_n ⟩

2 = |ϵ₀(z)|²

|ϵ_p|² X

n

⟨˜u₀;e⁽²⁾_n ⟩

2+O a^2−2h(4.6)

= O a^2−2h .