www.oeaw.ac.at

www.ricam.oeaw.ac.at

**Mathematical analysis of the** **acoustic imaging modality**

**using bubbles as contrast** **agents at nearly resonating**

**frequencies**

**A. Dabrowski, A. Ghandriche, M. Sini**

**RICAM-Report 2021-35**

MATHEMATICAL ANALYSIS OF THE ACOUSTIC IMAGING MODALITY USING BUBBLES AS CONTRAST AGENTS AT NEARLY RESONATING FREQUENCIES

ALEXANDER DABROWSKI^{∗}, AHCENE GHANDRICHE^{∗}AND MOURAD SINI^{‡}

Abstract. We analyze mathematically the acoustic imaging modality using bubbles as contrast agents. These bubbles are modeled by mass densities and bulk moduli enjoying contrasting scales. These contrasting scales allow them to resonate at certain incident frequencies. We consider two types of such contrasts. In the first one, the bubbles are light with small bulk modulus, as compared to the ones of the background, so that they generate the Minnaert resonance (corresponding to a local surface wave). In the second one, the bubbles have moderate mass density but still with small bulk modulus so that they generate a sequence of resonances (corresponding to local body waves).

We propose to use as measurements the far-fields collected before and after injecting a bubble, set at a given location point in the target domain, generated at a band of incident frequencies and at a fixedsingle backscattering direction. Then, we scan the target domain with such bubbles and collect the corresponding far-fields. The goal is to reconstruct both the, variable, mass density and bulk modulus of the background in the target region.

(1) We show that, for each fixed used bubble, the contrasted far-fields reach their maximum value at, incident, frequencies close to the Minnaert resonance (or the body-wave resonances depending on the types of bubbles we use). Hence, we can reconstruct this resonance from our data. The explicit dependence of these resonances in terms of the background mass density of the background allows us to recover it, i.e.

the mass density, in a straightforward way.

(2) In addition, this measured contrasted far-fields allow us to recover the total field at the location points of the bubbles (i.e. the total field in the absence of the bubbles). A numerical differentiation argument, for instance, allows us to recover the bulk modulus of the targeted region as well.

1. Introduction and statement of the results

Diffusion by highly contrasted small particles is of fundamental importance in several branches of applied sciences, as for example in material sciences and imaging. In this work, we focus on the acoustic imaging modality using microscaled bubbles as contrast agents, see [10, 17, 20, 21] for more details on related theoretical and experimental studies. We describe a modality using the contrasted scattered fields, by the targeted anomaly, measured before and after injecting microscaled bubbles. These bubbles are modeled by mass densities and bulk moduli enjoying contrasting scales. These contrasting scales allow them to resonate at certain incident frequencies. The main goal of this work is to analyze mathematically this contrasted scattered fields in terms of these scales with incident frequencies close to these resonances and derive explicit formulas linking the values of the unknown mass density and bulk modulus of the targeted region to the measured scattered fields.

To describe properly the mathematical model we are dealing with in this work, let us denote by D a small
particle inR^{3}of the formD:=εB+z, whereBis an open, bounded, simply connected set inR^{3}with Lipschitz
boundary, containing the origin, andz specifies the location of the particle. The parameterε >0 characterizes

2010Mathematics Subject Classification. 35R30, 35C20.

Key words and phrases. Acoustic imaging, bubbly media, Minnaert resonance, surface-wave resonances, Newtonian potential, body-wave resonances.

∗ RICAM, Austrian Academy of Sciences, Altenbergerstrasse 69, A-4040, Linz, Austria. Email: alexan- [email protected] This author is supported by the Austrian Science Fund (FWF): P 30756-NBL.

∗ RICAM, Austrian Academy of Sciences, Altenbergerstrasse 69, A-4040, Linz, Austria. Email:

[email protected] This author is supported by the Austrian Science Fund (FWF): P 30756-NBL.

‡RICAM, Austrian Academy of Sciences, Altenbergerstrasse 69, A-4040, Linz, Austria. Email: [email protected] This author is partially supported by the Austrian Science Fund (FWF): P 30756-NBL.

1

the smallness assumption on the particle. Let us consider a mass density (respectively, bulk modulus) that we note byρε(·) (respectively,kε(·)) of the form

ρε(x) :=

(ρ0(x), x∈R^{3}\D,

ρ1, x∈D, kε(x) :=

(k0(x), x∈R^{3}\D,
k1, x∈D,

whereρ1andk1are positive constants, whileρ0andk0are smooth enough functions which are constant outside of a bounded and smooth domain Ω. We denote respectively ¯ρ0 and ¯k0 to be the values of ρ0 andk0 outside Ω. Thusρ0andk0denote the density and bulk modulus of the background medium, and ρ1 andk1 denote the density and bulk modulus of the bubble respectively.

We are interested in the following problem describing the acoustic scattering by a bubble, see [11] and [12], given by the system

(1.1)

∇ · 1

ρ0

∇u

+ω^{2}
k0

u= 0 inR^{3}\D,

∇ · 1

ρ1

∇u

+ω^{2}
k1

u= 0 inD,
u|_{−}−u|_{+}= 0, on∂D,

1 ρ1

∂u

∂ν
_{−}

− 1 ρ0

∂u

∂ν
_{+}

= 0 on∂D,

where ω > 0 is a given frequency and ν denotes the external unit normal to ∂D. Here the total field is
u := u^{i}+u^{s}, where u^{i} denotes the incident field (we restrict to plane incident waves) and u^{s} denotes the
scattered waves which satisfy the following condition

(1.2) ∂u^{s}

∂|x|−iκ_{0}u^{s}=o
1

|x|

,|x| → ∞, (S.R.C).

We introduce the notation κ^{2}_{0} := ω^{2}ρ0/k0 and κ^{2}_{1} := ω^{2}ρ1/k1. The problem (1.1) is well posed, see [3, 4]

and [9]. In addition, the scattered fieldu^{s}can be expanded as
u^{s}(x, θ) = e^{iκ}^{0}^{|x|}

|x| u^{∞}(ˆx, θ) +O |x|^{−2}

, |x| →+∞,

where ˆx := x/|x| and u^{∞}(ˆx, θ) denotes the far-field pattern corresponding to the unit vectors ˆx, θ, i.e. the
incident and propagation directions respectively. We are interested in the regimes where the coefficients satisfy
the conditions:

ρ1=Cρε^{s}, s≥0 and k1=Ckε^{t}, t≥0,

with positive constants Cρ and Ck which are independent from ε, and real numbers s, t assumed to be non negative. The scattering problem described above models the acoustic wave diffracted in the presence of small bubbles. In this case, the parameterssandtfix the kind of medium we are considering, see [11, 12, 19] and [9].

Recall that the speed of propagation of the sound isc0 :=q

k0

ρ_{0} in the background and c1 :=q

k1

ρ_{1} and in the
bubble. We are interested in the following two regimes on the relative speed of propagation inside the bubble
as compared to the one in the background, i.e the ratio ^{c}_{c}^{1}

0.

(1) Moderate relative speed of propagation. In this case, we assume thats =t, then the relative speed of propagation is uniformly bounded from below and above or

c^{2}_{1}

c^{2}_{0} =ρ0k1

k0ρ1

= κ^{2}_{0}

κ^{2}_{1} ∼1, asε1.

(2) Small relative speed of propagation. In this case, we assume that s < t, then the relative speed of propagation is small or

c^{2}_{1}
c^{2}_{0} =κ^{2}_{0}

κ^{2}_{1} ∼ε^{t−s}, as ε1.

There is a major difference between these two regimes. In the first one, the contrast between the bubble and
the background comes only from the transmission coefficientsρ1/ρ0across the interface of the bubble. Hence,
there might be surface waves created by this contrast if it is pronounce enough. For certain scales ofρ1/ρ0, this
is indeed the case, as we will see it later. In the second regime, as the speed of propagation is very small inside
the bubble, then the sound wave (i.e. the fluctuation) slows down inside it and might create local spots, body
waves, even though the bubble has a small size. For certain scales of ^{c}_{c}^{2}^{1}2

0

this is indeed the case.

To highlight these differences, let us for the moment assume thatρ0 is constant everywhere in R^{3}. In this
case, the above problem can be equivalently formulated as

∆u+κ^{2}_{0}u= 0 inR^{3}\D,

∆u+κ^{2}_{1}u= 0 inD,
u|_{ư}ưu|_{+}= 0, on∂D,

1
ρ_{1}

∂u

∂ν ư

ư 1
ρ_{0}

∂u

∂ν
_{+}

= 0 on∂D,
uưu^{i} satisfies the SRC.

As we can see, the contrasts of the medium appear in the transmission conditions through the coefficient
1/ρ_{1} (or equivalently ρ_{0}/ρ_{1}), and through the speed of propagation, namely ρ_{0}/k_{0} and ρ_{1}/k_{1}. Based on the
Lippmann-Schwinger equation representation of the total fields, the second contrast appears on the (volumetric)
Newtonian potentialN_{D}^{ω} defined, fromL^{2}(D) toH^{2}(D), as

N_{D}^{ω}[f](x) :=

Z

D

Gω(x, y)f(y)dy,

while the first one appears on the (surface) Neumann-Poincar´e operator (K_{D}^{ω})^{?}defined, fromL^{2}(∂D) toL^{2}(∂D),
by

(K_{D}^{ω})^{?}[f](x) :=p.v
Z

∂D

∂Gω(x, y)

∂ν(x) f(y)dσ(y).

Precisely, the values of the fielduoutside the bubbleD is fully computable from the knowledge ofu(x), x∈D
and∂_{ν}u(x), x∈∂D. These last quantities are solutions of the following system of integral equations

(1.3) u(x)ưγω^{2}N_{D}^{ω}[u] (x) +αS_{D}^{ω}[∂_{ν}u] (x) =u^{i}(x), onD
and

(1.4) α

1
α+ρ_{0}

2 + (K_{D}^{ω})^{∗}

[∂νu]ưγω^{2}∂νưN_{D}^{ω}[u] (x) =∂νu^{i} on∂D,
whereS_{D}^{ω} is the single layer operator defined, fromL^{2}(∂D) toH^{3}^{2}(D), as

S_{D}^{ω}[f] (x) :=

Z

∂D

Gω(x, y)f(y)dσ(y),
andu^{i} is the incident field such that

(1.5) div

1 ρ0

∇u^{i}

(x) +ω^{2}
k0

u^{i}(x) = 0, x∈ D.

and we have adapted the succeeding notationsγ=βưαρ1/k1 andα:= 1/ρ1ư1/ρ0 withβ:= 1/k1ư1/k0. Here,Gω stands for the Green’s functions related to (1.5). More precisely,Gω is solution, in the distributional sense, of

(1.6) ∇

x ·h

ρ^{ư1}_{0} (x)∇

xGω(x, z)i
+ω^{2}

k_{0}Gω(x, z) =ưδ

z(x) for anyx, z∈R^{3}.
with the radiation conditions at infinity. For ˆxin the unit sphere, we note by

G^{∞}_{ω}(ˆx,·)

the far field associated toGω(·,·). Unless specified, we use the same notationGω(·,·) for the general setting whenρ0andk0 are locally variable.

Depending on the scales of the contrasts, we make the following observations.

(1) In the first regime, i.e s=t, we have γ ∼1, as ε <<1, and the Newtonian potential is negligible as
it scales asε^{2} as ε1. However, ifs=t= 2 then the contrasts on the mass densities, i.e. 1/α, can
approximate the spectrum of the Neumann-Poincar´e operatorK_{D}^{0}. For smooth domainD, this operator
defined onL^{2}(∂D) has a sequence of real eigenvalues accumulating at 0 in addition to the value−^{ρ}_{2}^{0}.
As the contrast is real, then we can only approximate the highest eigenvalue, which is−^{ρ}_{2}^{0}. This can
be done in this regime asα∼ε^{−2}. The frequencyωfor which this is possible is the Minnaert resonance
(corresponding to a surface wave type).

(2) In the second regime, if s < t, the high contrasts of the speed of propagation allow the Newtonian
operator to dominate the Neumann-Poincar´e operator. In addition, if we take t−s = 2, then the
contrast of the speed of propagation,γ∼ε^{−2}, will balance the scale of the Newtonian operator and we
might excite its eigenvalues. There is a discrete sequence of such eigenvalues (corresponding to local
body waves).

Microbubbles with scales fitting into the first regime are well known to exist in the nature. However, those related to the second regime, with small speed of propagation, are less known. Nevertheless, there are possibilities to artificially produce them, see the discussion in [24] and also in [23].

A first key observation in our analysis, which happens to be useful for the imaging later on, is that the Minnaert resonance and the sequence of body wave resonances are characterized by the bulk modulus of the bubble and the surrounding local mass density of the background. In addition, we show that the contrasted scattered fields reach their maximum values at, incident, frequencies close to the Minneart resonance (or the body-wave resonances depending on the types of bubbles we use). This allows us to recover these resonances by measuring the contrasted scattered waves at a band of incident frequencies but at a fixed single backscattering direction. A second key observation is that this measured contrasted scattered waves allows us to recover the total field at the location point of the bubble. Scanning the targeted region with such bubbles, we can recover the total field there up to a sign (i.e. the total field in the absence of the bubbles).

Based on these observations, we can reconstruct the density and the bulk modulus of the targeted region from the contrasted scattered waves (before and after injecting the bubbles) at a band of incident frequencies but at a fixed single backscattering direction. More details are given in section 2. Nevertheless, let us say it in short here that these contrasted scattered waves encodes the Minnaert and the body-waves resonance in its denominator and the total field in its numerator. From the first one, we extract the mass density while from the second one we derive the bulk modulus of the targeted region.

The following theorems translate these observations with more clear statements. We state the following conditions which are common to both the two results.

Conditions. Let Ωbe a bounded domain of diameterdiam(Ω)of order 1. Let also ρ_{0} andk_{0} be two functions
of classC^{1} and are constant outsideΩ. They are assumed to be positive functions. Let D:=z+ε B be a small
and Lipschitz smooth domain where z ∈ Ω away from its boundary. The relative diameter of D is small as
compared to the diameter ofΩ, i.e. ε

diam(Ω) << 1. The functions ρ_{0} and k_{0} are assumed to be independent
on the parameterε.

Theorem 1.1. Let the above Conditions be satisfied. In addition, let ρ_{1} and k_{1} be constants enjoying the
following scales

ρ1=ρ_{1}ε^{2}, k1=k1ε^{2} and k1

ρ_{1} ∼1, as ε <<1
andρ_{1} is large enough such thatmax

x∈Ωρ0(x)< ρ_{1}.

The solution of the corresponding problem (1.1), has the following expansions.

(1) The scattered field is approximated as

(1.7) u^{s}(·, θ, ω) =v^{s}(x, θ, ω)− ω^{2}ω_{M}^{2}

k1(ω^{2}−ω_{M}^{2} )|B|ε G_{ω}(x−z)v(z, θ, ω) +O ε^{2}
(ω^{2}−ω_{M}^{2} )^{2}

!

uniformly forxin a bounded domain away from D andθin the unit sphere.

(2) The farfield is approximated as

(1.8) u^{∞}(ˆx, θ, ω) =v^{∞}(ˆx, θ, ω)− ρ_{0}ω^{2}ω_{M}^{2}

k1(ω^{2}−ω^{2}_{M}) 4π|B|εv(z,−ˆx, ω)v(z, θ, ω) +O ε^{2}
(ω^{2}−ω_{M}^{2} )^{2}

! ,

uniformly forθ andxˆ in the unit sphere.

These expansions are valid under the condition that ε^{h}

(ω^{2}−ω^{2}_{M}) =O(1), h <1 asε <<1. Here, we have
(1.9) ωM :=ωM(z) :=

s 2k1

ρ0(z)µ∂B

, with µ∂B := 1

|∂B|

Z

∂B

Z

∂B

(x−y)·ν(x)

4π|x−y| dσ(x)dσ(y),

called the Minnaert frequency^{1}. In both the expansions v := v(x, θ, ω) = v^{∞}(ˆx, θ, ω) +u^{i}(x, θ) and v^{∞} :=

v^{∞}(ˆx, θ, ω)is the total field, and is the far field associated to the scattered field, solution of the problem (1.1) in
the absence of the bubbleD.

Note, from (1.9), that the Minnaert frequencyωM is such thatωM ∼1, in terms ofε, and sinceωapproaches ωM, we deduce that ω∼1, in terms ofε.

The first mathematical study of the Minnaert resonance was shown in [6] where it was estimated for bubbles injected in a homogeneous background. Later on, a series of works were devoted to its implications in different areas, see [3, 4, 7–9]. The approximations in (1.7) and (1.8) are extensions of those in [6] and [3] to the case when the background is heterogeneous (with variable mass density and bulk modulus). The surprising fact is that this resonance depends also on the surrounding background through its mass density.

To state the results related to the second regime, we first introduce, with some details, the Newtonian
operatorN^{ω}B defined fromL^{2}(B) to L^{2}(B) as^{2}

N^{ω}B[f](x) :=

Z

B

Γω(x, y)f(y)dy= Z

B

ρ0(z) e^{i ω}

q_{ρ}

0 (z) k0 |x−y|

4π|x−y| f(y)dy,
where Γ_{ω}(·,·) is the fundamental solution of

∆xΓω(x, z) +ω^{2}ρ0(z)

k_{0}(z) Γω(x, z) =−ρ0(z)δ

z(x), x, z ∈R^{3},

with radiations conditions at infinity. For ω = 0, the operator N^{0}B[·] is positive, compact, selfadjoint and by
ρ_{0}(z)λ^{B}_{n}, e^{B}_{n}

n∈Nwe design its sequence of eigenvalues with the corresponding eigenfunctions, where obviously
λ^{B}_{n}, e^{B}_{n}

n∈Nare related to the operator defined throughf(·)→R

B 1

4π|·−y|f(y)dy.

Theorem 1.2. Let the above Conditions be satisfied. In addition, let ρ1 and k1 be constants enjoying the following scales

(1.10) ρ1=ρ0(z) +O(ε^{j}), j >0, and k1=k1ε^{2} asε <<1.

In this regime, the solution of the problem (1.1), has the following expansions.

1Remark thatµ∂Bdepends only on the shape of the domainB.

2For convenience, we’ll omit the dependency notation ofN^{ω}B with respect toz.

(1) The scattered field has the approximation

(1.11) u^{s}(x, θ, ω) =v^{s}(x, θ, ω)− 1
k_{1}

ω^{2}ω^{2}_{n}

0

(ω^{2}−ω^{2}_{n}_{0})
Z

B

e^{B}_{n}_{0}^{2}

ε Gω(x;z)v(z, θ, ω) +O ε + ε^{1+min(1;j)}
ω^{2}−ω_{n}^{2}_{0}2

! ,

uniformly forxin a bounded domain away from D andθin the unit sphere.

(2) The farfield has the approximation (1.12)

u^{∞}(x, θ, ω) =v^{∞}(x, θ, ω)− ρ_{0}
4π k_{1}

ω^{2}ω^{2}_{n}

0

(ω^{2}−ω^{2}_{n}_{0})
Z

B

e^{B}_{n}

0

^{2}

ε v(z,−ˆx, ω)v(z, θ, ω) +O ε + ε^{1+min(1;j)}
ω^{2}−ω^{2}_{n}_{0}2

! ,

uniformly forθ andxˆ in the unit sphere.

These expansions are valid as soon as ε^{h}
ω^{2}−ω_{n}^{2}

0

=O(1), withh <min{1, j}, asε <<1, where

(1.13) ωn_{0}:=

s k1

ρ0(z)λ^{B}_{n}

0

.

Observe that R

Be^{B}_{n}_{0}2

means P

l

R

Be^{B}_{l} 2

, wherel is such that NB^{0}[e^{B}_{l} ] =λ^{B}_{n}_{0}e^{B}_{l} . ^{3}

Here again v := v(x, θ, ω) and v^{∞} := v^{∞}(ˆx, θ, ω) is the total field, and is the far field associated to the
scattered field, solution of the problem (1.1) in the absence of the bubbleD.

The body-wave resonances have been characterized already in [5, 22] in the framework of dielectric nanopar- ticles, in the scalar model related to the TM regime of the electromagnetic scattering, with a homogeneous background. There, the contrast comes from the dielectric nanoparticles with high permittivity and moderate permeability. In our context, the contrast comes from the fact that the density of the bubble is moderate while its bulk is still small. At the mathematical level, our formulas in (1.11) extend those in [5] to the case of the acoustic model, i.e. a divergence form model, with heterogeneous background. As we have said above, such bubble’s contrasts might not be available in nature but can be artificially designed, see [24].

We finish this section with the following observations.

(1) The approximations in Theorem 1.1 are similar to the ones in Theorem 1.2 up to the multiplicative factor appearing in the dominating term. The additional term O(ε) appearing in the error of the approximations (1.11) and (1.12) can be removed as follows:

(a) The scattered fields are approximated as
u^{s}(x, θ, ω) =v^{s}(x, θ, ω) +ω^{2}

k_{1} Gω(x, z)v(z, θ, ω)
Z

D

W(x)dx+O ε^{1+min(1;j)}
ω^{2}−ω^{2}_{n}_{0}^{2}

!

(b) The farfields are approximated as
u^{∞}(ˆx, θ, ω) =v^{∞}(ˆx, θ, ω) + ρ_{0}ω^{2}

4π k1

v(z,−ˆx, ω)v(z, θ, ω) Z

D

W(x)dx+O ε^{1+min(1;j)}
ω^{2}−ω^{2}_{n}

0

2

!

whereW := I−γ ω^{2}N^{0}B

−1

[1]. The termO(ε) appearing in (1.11) and and (1.12), is due to the fact that, see (4.20),

Z

D

W(x)dx=−ω_{n}^{2}_{0}
R

De^{D}_{n}_{0}(x)dx2

ω^{2}−ω_{n}^{2}_{0} +O ε^{3}
.

3This can be seen in (4.20).

(2) In both Theorem 1.1 and Theorem 1.2, the error terms contain the term _{(ω}_{2}_{−ω}^{ε}^{2}2

M)^{2} or _{(ω}2−ω^{ε}^{2}_{n}^{2}_{0})^{2}, respec-
tively. This is not optimal and we believe it can improved to allow us to takeω^{2}−ω^{2}_{M} and ω^{2}−ω^{2}_{n}_{0}
of the order of ε or less. This can be done at the expanse of improving the dominating term of the
corresponding approximation. This can be seen from the proofs when the background is homogeneous.

In the inhomogeneous background case, its justification makes the computations rather more involved and we prefer to skip this and stick to the results as stated above for clarity.

(3) Finally, we do believe that the condition max

x∈Ωρ0(x)< ρ_{1}used in Theorem 1.1 and the condition (1.10)
appearing in Theorem 1.2 might be removed.

2. An application to the acoustic imaging using resonating bubbles

Based on the expansions given in Theorem 1.1 and Theorem 1.2, in particular (1.8) and (1.12), we design
the following imaging procedure to reconstruct the mass density ρ_{0} and bulk modulus k_{0} inside the bounded
domain Ω where they are variable. This procedure is based on the following measured data.

Let [ωmin, ωmax] be an interval of possible incident frequencies under the following conditions
ω_{min}≤

v u u t

2k1

max

z∈Ωρ_{0}(z)µ_{∂B} ≤
v
u
u
t

2k1

min

z∈Ωρ_{0}(z)µ_{∂B} ≤ω_{max}.

This condition makes sense as soon as we know a priori a lower bound and an upper bound of the unknown mass densityρ0.

(1) Collect the farfields before injecting the bubble D, i.e. measure the backscattered farfield at a single
incident waveθ and a band of frequenciesω∈[ω_{min}, ω_{max}] : v^{∞}(−θ, θ, ω).

(2) Collect the farfield after injecting the bubbleD, centered at the pointz∈Ω, i.e. measure the backscat-
tered farfield at a single incident waveθand a band of frequenciesω∈[ωmin, ωmax] : u^{∞}(−θ, θ, ω, z).

The imaging procedure goes as follows. We set

(2.1) I(ω, z) :=u^{∞}(−θ, θ, ω, z)−v^{∞}(−θ, θ, ω)

as the imaging functional, remembering that the incident angle θ is fixed. We have the following properties from (1.8)

(2.2) I(ω, z)∼ − ρ_{0}

4π k_{1}

ω^{2}_{M}

(ω^{2}−ω_{M}^{2} (z))|B|ε [v(z, θ, ω)]^{2}.
We divide this procedure into two steps:

(1) Step 1. From this expansion, we recover ω^{2}_{M}(z) as the frequency for which the imaging functionω →
I(ω, z) gets its largest value. From the estimation of this resonance ω_{M}^{2} (z), we reconstruct the mass
density at the center of the injected bubblez, based on (1.9), as follows:

(2.3) ρ0(z) = 2k_{1}

ω_{M}^{2} (z)µ∂B

.

Scanning the domain Ω by such bubbles, we can estimate the mass density there.

(2) Step 2. To estimate now the bulk modulus, we go back to (2.2) or (1.8), and derive the values of the
totale field [v(z, θ, ω)]^{2}. This field corresponds to the model without the bubble. Hence, we have at
handv(z, θ, ω) forz∈Ω up to a sign (i.e. we know the modulus and the phase up to a multiple ofπ).

Use the equation ∇ ·ρ^{−1}_{0} ∇v+ω^{2}k_{0}^{−1}v = 0 to recover the values of k0 in the regions where v does
not change sign. This can be done by numerical differentiation for instance. Other ways are of course
possible to achieve this second step. In addition, we have at hand multiple frequency internal data.

The procedure described above uses the Minnaert resonance. The key point to recover the mass density is
the explicit dependance of this resonance on the value of the mass density on it’s ’center’, see (1.9). We can do
the same work using the sequence of resonances coming from the second regime, see (1.13). Therefore, we may
recoverρ_{0}and k_{0} for this regime as well. However, for technical reasons, we need to know the mass density as
we use the condition (1.10). But as we have said earlier, we believe that this condition might be removed.

3. Proof of Theorem 1.1

We divide the proof into two steps. In the first step, we provide the expansions in the case when the background is homogeneous. This allows to show the key parts in localizing the resonance and computing the scattered fields from incident frequencies close to these resonances. In the second step, we deal with the case when the background is heterogeneous and show how this perturbation influences the derivation of the expansions and the resonances as well.

3.1. Constant coefficients. We assume here that bothρ0 andk0 are constants everywhere inR^{3}. We recall
thatρ1= ¯ρ1ε^{2}, k1= ¯k1ε^{2}, where ¯ρ1,¯k1 do not depend onε. In this case, it is immediate to show that

Gω(x, y) =ρ0

e^{iκ}^{0}^{|x−y|}

4π|x−y|, x6=y, where κ0:=ωp ρ0/k0. Letube the solution of (1.1). From the Lippman-Schwinger representation we have

(3.1) u(x)−αdiv

x

Z

D

Gω(x, y)∇u(y)dy−β ω^{2}
Z

D

Gω(x, y)u(y)dy=u^{i}(x),
where α:= 1/ρ_{1}−1/ρ_{0} andβ:= 1/k_{1}−1/k_{0}.

Since∇

xG_{ω}(x, y) =−∇

yG_{ω}(x, y), by integration by parts and (1.6) we have
divx

Z

D

G_{ω}(x, y)∇u(y)dy=−ω^{2}ρ_{1}
k1

Z

D

G_{ω}(x, y)u(y)dy−
Z

∂D

G_{ω}(x, y)∂_{ν}u(y)dσ(y),
so (3.1) becomes

(3.2) u(x)−γω^{2}

Z

D

G_{ω}(x, y)u(y)dy+α
Z

∂D

G_{ω}(x, y)∂_{ν}u(y)dσ(y) =u^{i}(x),

whereγ=β−αρ_{1}/k_{1}.Taking the normal derivative asx→∂Dfrom insideD, from the jump relations of the
derivative of the single layer potential we obtain

(3.3)

1 + αρ_{0}
2

∂νu(x)−γω^{2}∂ν−

Z

D

Gω(x, y)u(y)dy+α(K_{D}^{ω})^{∗}[∂νu] (x) =∂νu^{i}(x).

Notice that due to the scaling of ρ_{1} and k_{1}, we have γ = O(1) as ε → 0. Expanding in z the fundamental
solution, we obtain forxaway fromD,

Z

D

Gω(x, y)u(y)dy=Gω(x, z) Z

D

u(y)dy+O

ε^{5}^{2}kuk_{L}2(D)

, as by the Cauchy-Schwartz inequality and the fact that|y−z|=O(ε) we have

Z

D

(y−z)u(y)dy

≤ k · −zk_{L}2(D)kuk_{L}2(D)=O

ε^{5}^{2}kuk_{L}2(D)

. In the same way, we have

Z

∂D

|y−z| ∂νu(y) dσ(y).ε^{2} k∂νuk_{L}2(∂D),
so that

Z

∂D

Gω(x, y)∂νu(y)dσ(y) =Gω(x, z) Z

∂D

∂νu(y)dσ(y) +O

ε^{2}k∂νuk_{L}2(∂D)

.

Therefore, we can rewrite (3.2) as
u^{s}(x) =γ ω^{2}Gω(x, z)

Z

D

u(y)dy −α Gω(x, z) Z

∂D

∂νu(y)dσ(y) + O

ε^{5}^{2}kuk_{L}2(D)+k∂νuk_{L}2(∂D)

. From the equation satisfied byu, see for instance (1.1), and the divergence theorem, we have

(3.4)

Z

D

u(y)dy=−k1

ω^{2}
Z

D

∇ · 1

ρ_{1}∇u

(y)dy=− k1

ω^{2}ρ_{1}
Z

∂D

∂νu(y)dσ(y), then

(3.5) u^{s}(x) =−

α+γk_{1}
ρ1

G_{ω}(x, z)
Z

∂D

∂_{ν}u(y)dσ(y) +O

ε^{5}^{2}kukL^{2}(D)+k∂νuk_{L}2(∂D)

.

Now, we derive the dominating term of R

∂D∂νu dσ and estimate kuk_{L}2(D) and k∂νuk_{L}2(∂D) in terms of ε.

Let us consider first the case whenγ= 0. In this case, the equation (3.3) becomes
(3.6) ((1/α+ρ0/2)I+ (K_{D}^{ω})^{∗}) [∂νu] =α^{−1}∂νu^{i},
and we can rewrite it as

(3.7) (1/α+ρ0/2)I+ (K_{D}^{0})^{∗}

[∂νu] + (K_{D}^{ω})^{∗}−(K_{D}^{0})^{∗}

[∂νu] =α^{−1}∂νu^{i}.
Let

(3.8)
A∂D:= ρ^{2}_{0}

k0

µ∂D, µ∂D:= 1

|∂D|

Z

∂D

Z

∂D

(x−y)·ν(x)

4π|x−y| dσ(x)dσ(y)and A(y) := ρ^{2}_{0}
k0

Z

∂D

(x−y)·ν(x) 4π|x−y| dσ(x).

By the divergence theorem we have A∂D >0, and it is immediate thatA−A∂D has average zero along∂D.

ExpandingGω(x, y) in terms of|x−y|, we obtain
(K_{D}^{ω})^{∗}[∂νu](x) :=

Z

∂D

∂ Gω(x, y)

∂ν(x) ∂νu(y)dσ(y)

= Z

∂D

ρ0

4π

"

(x−y)·ν(x)

|x−y|^{3} −1 +iκ_{0}|x−y|

∞

X

n=0

(iκ0|x−y|)^{n}
n!

#

∂_{ν}u(y)dσ(y)

= (K_{D}^{0})^{∗}[∂_{ν}u](x)−κ^{2}_{0}ρ_{0}
8π

Z

∂D

(x−y)·ν(x)

|x−y| ∂_{ν}u(y)dσ(y)

− iκ^{3}_{0}ρ_{0}
12π

Z

∂D

(x−y)·ν(x)∂νu(y)dσ(y) +O ε^{3}k∂νukL^{2}(∂D)

, (3.9)

and integrating (3.7) on∂D, as K_{D}^{0}[1] =−ρ0/2, see Appendix (5.3), we obtain
1

α−κ^{2}_{0}k0

2ρ0

A∂D

Z

∂D

∂νu(x)dσ(x) = 1 α

Z

∂D

∂νu^{i}(x)dx+iκ^{3}_{0}ρ0

12π Z

∂D

Z

∂D

(x−y)·ν(x)∂νu(y)dydx
+ κ^{2}_{0}k_{0}

2ρ0

Z

∂D

(A(y)−A_{∂D})∂_{ν}u(y)dy+O ε^{5}k∂_{ν}uk_{L}2(∂D)

. (3.10)

We can estimate the integral which containsA(·)−A_{∂D} by rewriting
Z

∂D

(A(y)−A_{∂D})∂_{ν}u(y)dσ(y) ^{(3.6)}= α^{−1}
Z

∂D

(A(y)−A_{∂D}) (ρ_{0}/2 + 1/α+ (K_{D}^{ω})^{∗})^{−1}[∂_{ν}u^{i}]

(y)dσ(y)

= α^{−1}
Z

∂D

(ρ_{0}/2 + 1/α+K_{D}^{ω})^{−1}[A(·)−A_{∂D}]

(y)∂_{ν}u^{i}(y)dσ(y)

≤ α^{−1}

(ρ0/2 + 1/α+K_{D}^{ω})^{−1}[A(·)−A∂D]
_{L}2(∂D)

∂νu^{i}

_{L}_{2}_{(∂D)}=O ε^{6}
,
(3.11)

the last equality being a consequence of the fact that (ρ0/2 + 1/α+K_{D}^{ω})^{−1}does not scale on L^{2}_{0}(∂D) :={f ∈
L^{2}(∂D) :R

∂Df dσ= 0}, and AandA∂D scale both asε^{2}. Then (3.10) becomes
1

α−iκ^{3}_{0}|D|ρ_{0}

4π −κ^{2}_{0}k_{0}
2ρ0

A_{∂D}
Z

∂D

∂_{ν}udσ= 1
α

Z

∂D

∂_{ν}u^{i}dσ+O ε^{5}k∂_{ν}uk_{L}2(∂D)+ε^{6}
,
where we have used the fact that R

∂D(x−y)·ν(x)dσ(x) = R

Ddiv [x−y]dx = 3|D|. Then, multiplying by α
(which scales likeε^{−2}), we obtain the expression of the following dominating term ofR

∂D∂νudσ,

1−i α κ^{3}_{0}|D|ρ_{0}

4π −α κ^{2}_{0}k_{0}
2ρ0

A_{∂D}
Z

∂D

∂_{ν}u dσ=
Z

∂D

∂_{ν}u^{i}dσ+O ε^{3}k∂_{ν}uk_{L}2(∂D)+ε^{4}
.

In the general case ofγ6= 0, instead of identity (3.6), we have 1

α+ρ0

2 + (K_{D}^{ω})^{∗}

[∂_{ν}u](x)−ω^{2}γ
α ∂_{ν−}

Z

D

G_{ω}(x, y)u(y)dy=α^{−1}∂_{ν}u^{i}(x).

Integrating on∂D, and integrating by parts the last integral, we obtain Z

∂D

1 α+ρ0

2 +(K_{D}^{ω})^{∗}

[∂νu](x)dσ(x)+ω^{2}γρ0

α
ω^{2}

k0

Z

D

Z

D

Gω(x, y)u(y)dy dx+ Z

D

u(x)dx

=α^{−1}
Z

∂D

∂νu^{i}(x)dσ(x).

Then, with the same estimates as in (3.10), we obtain 1

α−iκ^{3}_{0}|D|ρ_{0}

4π −κ^{2}_{0}k_{0}
2ρ0

A_{∂D}
Z

∂D

∂_{ν}u(x)dσ(x) + ω^{2}γρ_{0}
α

Z

D

u(x)dx=α^{−1}
Z

∂D

∂_{ν}u^{i}(x)dσ(x)

− ω^{2}γ κ^{2}_{0}
α

Z

D

Z

D

Gω(x, y)u(y)dydx+error, (3.12)

where

error:=O ε^{5}k∂_{ν}uk_{L}2(∂D)+ε^{6}
.

Next, with help of the Cauchy-Schwartz inequality, we estimate the double volume integral as

ω^{2}γ κ^{2}_{0}
α

Z

D

Z

D

Gω(x, y)u(y)dydx

.ε^{11}^{2} kuk_{L}2(D),
then, the equation (3.12) takes the following form

1

α−iκ^{3}_{0}|D|ρ0

4π −κ^{2}_{0}k_{0}
2ρ0

A∂D

Z

∂D

∂νu(x)dσ(x) +ω^{2}γρ_{0}
α

Z

D

u(x)dx= 1 α

Z

∂D

∂νu^{i}(x)dσ(x) +r,
where

(3.13) r:=O

ε^{11}^{2}kukL^{2}(D)+ε^{5}k∂νukL^{2}(∂D)+ε^{6}
.
We use (3.4) and the fact that ∆u^{i}=−κ^{2}_{0}u^{i} to obtain

(3.14)

1

α−iκ^{3}_{0}|D|ρ0

4π −ω^{2}

2 A∂D−γk1ρ0

αρ1

Z

∂D

∂νu=−ω^{2}ρ0

αk0

Z

D

u^{i}+r=−ω^{2}ρ0

α k0

|D|u^{i}(z) +r.

From the definition ofκ_{0}, recall thatκ_{0}:=ωp

ρ_{0}/k_{0}, we can see that the term between parenthesis on the left
hand side of the previous equation is cubic polynomial function onω. Now, we define the Minnaert frequency
ωM to be the dominating part of the zero of this cubic polynomial function. To find the dominant part of the
zero is equivalent to solve

ω^{2}= 2
A∂D

1

α−iκ^{3}_{0}|D|ρ0

4π −γk_{1}ρ_{0}
αρ1

= 2

A∂D

1 α

k_{1}ρ_{0}
ρ1k0

+O(ε) = 2k_{1}
µ∂Bρ0

+O(ε),

where the last two equalities are established using the definition of α, β, γ and the scales of |D| and A∂D. Setting,

(3.15) ω^{2}_{M} := 2k_{1}

µ∂Bρ0

,

we have (3.16)

1

α−iκ^{3}_{0}|D|ρ0

4π −ω^{2}

2 A∂D−γk_{1}ρ_{0}
αρ1

=ε^{2}ρ_{0}k_{1}
k0

ω_{M}^{2} −ω^{2}

ω^{2}_{M} +O ε^{3}
.
The equation (3.14), using (3.16), takes the following form

ε^{2}ρ0k1

k_{0}

ω^{2}_{M}−ω^{2}
ω^{2}_{M}

Z

∂D

∂_{ν}u(x)dσ(x) +O ε^{3}
Z

∂D

∂_{ν}u(x)dσ(x) =−ω^{2}ρ0

α k_{0} |D|u^{i}(z) +r.

By Cauchy-Schwartz inequality we estimate the second term on the left hand side equation asO

ε^{4} k∂νuk

L^{2}(∂D)

. Then, we obtain the following formula

Z

∂D

∂νu(x)dσ(x) = κ^{2}_{1} |D| u^{i}(z)

(ω^{2}−ω_{M}^{2} ) +O ε^{7}+r+ε^{4} k∂νuk

L^{2}(∂D)

ε^{2}(ω^{2}−ω_{M}^{2} )

!

(3.13)

= ω^{2}_{M}κ^{2}_{1} |D|u^{i}(z)

(ω^{2}−ω^{2}_{M}) +O ε^{4}+ε^{7}^{2} kuk

L^{2}(D)+ε^{2} k∂νuk

L^{2}(∂D)

(ω^{2}−ω_{M}^{2} )

! (3.17)

= O

ε^{3}
(ω^{2}−ω_{M}^{2} )

+O ε^{4}+ε^{7}^{2} kuk_{L}2(D)+ε^{2} k∂νuk_{L}2(∂D)

(ω^{2}−ω_{M}^{2} )

! . (3.18)

To estimate the error term, on the last expression, we need the following a priori estimates.

Proposition 3.1. Foru=u^{i}+u^{s}, solution of(1.1), it holds

(3.19) k∂νuk_{L}2(∂D)=O ε^{2}

ω^{2}−ω^{2}_{M}

, and

(3.20) kukL^{2}(D)=O ε^{3}^{2}

ω^{2}−ω^{2}_{M}

! ,

under the condition thatε/ ω^{2}−ω_{M}^{2}

is small enough.

Proof. Let us indicate asC a generic constant independent ofε. From (1.3) we have

(3.21) I−γω^{2}N_{D}^{ω}

[u] +α S_{D}^{ω}[∂νu] =u^{i}.

Sinceγ =O(1) and thus kN_{D}^{ω}k_{L} −−−→^{ε→0} 0, for ε small enough we have thatI−γω^{2}N_{D}^{ω} is invertible, so (3.21)
takes the following form

u=−α(I−γω^{2}N_{D}^{ω})^{−1}[S_{D}^{ω}[∂νu]] + (I−γω^{2}N_{D}^{ω})^{−1}[u^{i}].

Taking theL^{2}(D)-norm in both side of the last equation and using the fact that

I−γω^{2}N_{D}^{ω}−1

L(L^{2}(D))≤C
to obtain

(3.22) kuk_{L}2(D)≤α CkS_{D}^{ω}[∂_{ν}u]k_{L}2(D)+Cku^{i}k_{L}2(D).

In order to finish the last estimation we need to show precisely how the single layer potential scales. For this, by definition, we have

S_{D}^{ω}[f]

2

L^{2}(D) :=

Z

D

Z

∂D

Gω(x, y)f(y)dy

2

dx, ∀f ∈L^{2}(∂D)

= ε^{5}
Z

B

Z

∂B

G_{εω}(η−ξ) ˜f(ξ)dξ

2

dη:=ε^{5}

S_{B}^{εω} [ ˜f]

2
L^{2}(B)

(3.23)

and from the continuity ofS_{B}^{εω} fromL^{2}(∂B) toH^{3}^{2}(B) we have that

S_{D}^{ω}[f]

2

L^{2}(D)=ε^{5}
S_{B}^{εω}[ ˜f]

2

L^{2}(B)≤ε^{5}C
f˜

2

L^{2}(∂B)=ε^{3}C
f

2
L^{2}(∂D),

in particular

(3.24)

S_{D}^{ω}[∂νu]

_{L}_{2}_{(D)}≤ε^{3}^{2} C

∂νu
_{L}_{2}_{(∂D)}.
Combining (3.22) and (3.24), we obtain

(3.25) kuk_{L}2(D)≤ε^{−}^{1}^{2}Ck∂νuk_{L}2(∂D)+Cku^{i}k_{L}2(D).

To manage the termk∂νukL^{2}(∂D)we use the boundary integral equation given by (1.4), to write
(3.26) ∂νu=α^{−1}

1 α+ρ0

2 + (K_{D}^{ω})^{∗}
^{−1}

∂νu^{i}
+ω^{2}γ

α 1

α+ρ0

2 + (K_{D}^{ω})^{∗}
^{−1}

[∂νN_{D}^{ω}[u]] on ∂D.

Next, for convenience, we set

T :=

1
α+ρ_{0}

2 + (K_{D}^{ω})^{∗}
−1

and we rewrite (3.26) as

∂u

∂ν = 1 α T

∂u^{i}

∂ν − 1

|∂D|

Z

∂D

∂u^{i}

∂ν

+ 1

|∂D|

Z

∂D

∂u^{i}

∂ν 1 α T [1]

+ ω^{2}γ
α T

∂N_{D}^{ω}[u]

∂ν^{−} − 1

|∂D|

Z

∂D

∂N_{D}^{ω}[u]

∂ν^{−}

+ω^{2}γ
α

1

|∂D|

Z

∂D

∂N_{D}^{ω}[u]

∂ν^{−} T [1].
(3.27)

Since −^{ρ}_{2}^{0} is an eigenvalue of K_{D}^{0} with associated eigenspace consisting of constant functions, we have the
estimates

(3.28) kTk_{L}_{(L}2(∂D)) =

1 α+ρ0

2 + (K_{D}^{ω})^{∗}
−1

_{L}_{(L}_{2}_{(∂D))}

≤Cα,

and on the space of functions with zero average we have

(3.29) kTk_{L}_{(L}2

0(∂D))=

1 α+ρ0

2 + (K_{D}^{ω})^{∗}
−1

_{L}_{(L}2

0(∂D))

≤C.

Now, take theL^{2}(∂D)-norm in both side of (3.27), with the help of (3.28) and (3.29) we obtain

∂u

∂ν

_{L}_{2}_{(∂D)} . α^{−1}

∂u^{i}

∂ν − 1

|∂D|

Z

∂D

∂u^{i}

∂ν
_{L}2

0(∂D)

+ 1

|∂D|

Z

∂D

∂u^{i}

∂ν

k1k_{L}2(∂D)

+ α^{−1}

∂N_{D}^{ω}[u]

∂ν^{−} − 1

|∂D|

Z

∂D

∂N_{D}^{ω}[u]

∂ν^{−}
L^{2}_{0}(∂D)

+ 1

|∂D|

Z

∂D

∂N_{D}^{ω}[u]

∂ν^{−}

k1k_{L}2(∂D).
(3.30)

Obviously, we have (3.31)

Z

∂D

∂u^{i}

∂ν

= Z

D

∆u^{i}

=|κ^{2}_{0}|
Z

D

u^{i}

=O ε^{3}
and, by the triangular inequality and the smoothness of∂νu^{i}, we obtain
(3.32)

∂u^{i}

∂ν − 1

|∂D|

Z

∂D

∂u^{i}

∂ν
_{L}2

0(∂D)

.

∂u^{i}

∂ν
_{L}2

0(∂D)

=O(ε). We also have, recalling the definition of the Green function,

Z

∂D

∂N_{D}^{ω}[u]

∂ ν (x)dx = −ρ0

Z

D

u(x)dx−ω^{2}ρ0

k_{0}
Z

D

Z

D

Gω(x, y)u(y)dydx

= −ρ0

Z

D

u(x)dx+O

ε^{7}^{2} kuk_{L}2(D)

(3.4)

= k1ρ0

ω^{2}ρ_{1}
Z

∂D

∂νu(x)dσ(x) +O

ε^{7}^{2} kuk_{L}2(D)