Sub-wavelength focusing of acoustic waves in bubbly media

TL;DR

This paper investigates how bubbly media can be used to achieve super-focusing of acoustic waves by analyzing wave scattering near the Minnaert resonance frequency, combining theoretical derivations and numerical experiments.

Contribution

It introduces a new theoretical framework for modeling acoustic scattering in bubbly media and demonstrates super-focusing capabilities through numerical simulations.

Findings

Super-focusing of acoustic waves is achievable in bubbly media.

Derived scattering functions for single bubbles and multiple scattering approximations.

Numerical experiments confirm super-focusing potential near Minnaert resonance.

Abstract

The purpose of this paper is to investigate acoustic wave scattering by a large number of bubbles in a liquid at frequencies near the Minnaert resonance frequency. This bubbly media has been exploited in practice to obtain super-focusing of acoustic waves. Using layer potential techniques we derive the scattering function for a single spherical bubble excited by an incident wave in the low frequency regime. We then derive the point scatter approximation for multiple scattering by N bubbles. We describe several numerical experiments based on the point scatterer approximation that demonstrate the possibility of achieving super-focusing using bubbly media.

Equations

$$D:=\left\{x\in {\mathbb{R}}^3,|x|\le R\right\}.$$

$$\displaystyle\left\{\begin{array}[]{ll}&\displaystyle\nabla\cdot\frac{1}{\rho_{w}}\nabla u+\frac{\omega^{2}}{\kappa_{w}}u=0\quad\mbox{in }{\mathbb{R}}^{3}\backslash D,\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr&\displaystyle\nabla\cdot\frac{1}{\rho_{b}}\nabla u+\frac{\omega^{2}}{\kappa_{b}}u=0\quad\mbox{in }D,\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr&\displaystyle u_{+}-u_{-}=0\quad\mbox{on }\partial D,\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr&\displaystyle\frac{1}{\rho_{w}}\frac{\partial u}{\partial\nu}\bigg{|}_{+}-\frac{1}{\rho_{b}}\frac{\partial u}{\partial\nu}\bigg{|}_{-}=0\quad\mbox{on }\partial D,\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr&u^{s}:=u-u^{i}\,\,\,\mbox{satisfies the Sommerfeld radiation condition.}\end{array}\right.$$

$$c_w=\sqrt{\frac{{\kappa }_w}{{\rho }_w}},c_b=\sqrt{\frac{{\kappa }_b}{{\rho }_b}},k_w=\frac{\omega }{c_w},k_b=\frac{\omega }{c_b}.$$

$$\delta =\frac{{\rho }_b}{{\rho }_w},$$

$${\mathcal{S}}_D^k[\psi ](x)={\int }_{\partial D}G(x,y,k)\psi (y)d\sigma (y),x\in \partial D,$$

$$G(x,y,k)=G(x-y,k)=-\frac{e^{ik|x-y|}}{4\pi |x-y|},$$

$${\mathcal{K}}_D^{k,\ast }[\psi ](x)={\int }_{\partial D}\frac{\partial G(x,y,k)}{\partial \nu (x)}\psi (y)d\sigma (y),x\in \partial D.$$

$$u(x)=\left\{\begin{array}[]{lr}u^{in}+\mathcal{S}_{D}^{k_{w}}[\psi_{w}],&\quad x\in{\mathbb{R}}^{3}\backslash\bar{D},\\ \mathcal{S}_{D}^{k_{b}}[\psi_{b}],&\quad x\in{D},\end{array}\right.$$

$${\mathcal{A}}_D(\omega ,\delta )[\Psi ]=F,$$

$${\mathcal{A}}_D(\omega ,\delta )=\left(\begin{array}{cc}{\mathcal{S}}_D^{k_b}& -{\mathcal{S}}_D^{k_w}\\ \frac{1}{\delta }(-\frac{1}{2}Id+{\mathcal{K}}_D^{k_b,\ast })& -(\frac{1}{2}Id+{\mathcal{K}}_D^{k_w,\ast })\end{array}\right),\Psi =\left(\begin{array}{c}{\psi }_b\\ {\psi }_w\end{array}\right),F=\left(\begin{array}{c}{u}^{in}\\ \frac{\partial u^{in}}{\partial \nu }\end{array}\right).$$

$$X_1:=\mathrm{Vect}\left\{\left(\begin{array}{c}{\text{mathds}1}_{\partial D}\\ 0\end{array}\right),\left(\begin{array}{c}0\\ {\text{mathds}1}_{\partial D}\end{array}\right)\right\}.$$

$${\mathcal{A}}_D[\Psi ]=\left(\begin{array}{c}{\text{mathds}1}_{\partial D}\\ 0\end{array}\right),$$

$${\mathcal{S}}_D^k[{\text{mathds}1}_{\partial D}]=-\frac{{\mathrm{e}}^{\mathrm{i}kR}}{k}\sin (kR),$$

$${\left({\mathcal{K}}_D^k\right)}^{\ast }[{\text{mathds}1}_{\partial D}]=\frac{\sin (kR)}{kR}{\mathrm{e}}^{\mathrm{i}Rk}-\frac{{\mathrm{e}}^{2\mathrm{i}kR}}{2}.$$

$$A_D=\left(\begin{array}{cc}-\frac{{\mathrm{e}}^{\mathrm{i}{x}_b}}{x_b}\sin (x_b)& \frac{{\mathrm{e}}^{\mathrm{i}{x}_w}}{x_w}\sin (x_w)\\ & \\ \frac{{\mathrm{e}}^{\mathrm{i}{x}_b}}{\delta }\left(\frac{\sin (x_b)}{x_b}-\cos \left(x_b\right)\right)& -{\mathrm{e}}^{\mathrm{i}{x}_w}\left(\frac{\sin (x_w)}{x_w}-\mathrm{i}\sin \left(x_w\right)\right)\end{array}\right).$$

$$A_D\left(\begin{array}{c}{C}_b\\ C_w\end{array}\right)=\left(\begin{array}{c}\frac{u^{in}(0)}{R}\\ 0\end{array}\right).$$

$$\det (A_D)={\mathrm{e}}^{\mathrm{i}(x_b+x_w)}\sin (x_w)\left(\frac{\sin (x_b)}{x_b}\left(\frac{1}{x_w}-\mathrm{i}\right)-\frac{1}{\delta x_w}\left(\frac{\sin (x_b)}{x_b}-\cos (x_b)\right)\right),$$

$${\left(A_D\right)}^{-1}=\frac{-1}{\det (A_D)}\left(\begin{array}{cc}{\mathrm{e}}^{\mathrm{i}{x}_w}\left(\frac{\sin (x_w)}{x_w}-\mathrm{i}\sin \left(x_w\right)\right)& \frac{{\mathrm{e}}^{\mathrm{i}{x}_w}}{x_w}\sin (x_w)\\ & \\ \frac{{\mathrm{e}}^{\mathrm{i}{x}_b}}{\delta }\left(\frac{\sin (x_b)}{x_b}-\cos \left(x_b\right)\right)& \frac{{\mathrm{e}}^{\mathrm{i}{x}_b}}{x_b}\sin (x_b)\end{array}\right).$$

$$C_w=\left[\frac{{\mathrm{e}}^{-\mathrm{i}{x}_w}\left(\frac{1}{x_b}-\cot \left(x_b\right)\right)}{\sin (x_w)\left(\frac{1}{x_w}\left(\frac{1}{x_b}-\cot (x_b)\right)-\frac{\delta }{x_b}\left(\frac{1}{x_w}-\mathrm{i}\right)\right)}\right]\frac{u^{in}(0)}{R}.$$

$$\forall x\in {\mathbb{R}}^3\setminus D,{\mathcal{S}}_D^{k_w}[{\text{mathds}1}_{\partial D}](x)=-\frac{R\sin (k_{w}R)}{k_w}\frac{{\mathrm{e}}^{\mathrm{i}{k}_w|x|}}{|x|}.$$

$$u_{\mathrm{rad}}(|x|=R)=-C_{w}R\left(\frac{\sin (x_w)}{x_w}{\mathrm{e}}^{\mathrm{i}{x}_w}\right)=\left(\frac{x_b\cot \left(x_b\right)-1}{1-\delta -x_b\cot (x_b)+\mathrm{i}\delta \frac{c_b}{c_w}{x}_b}\right)u^{in}(0).$$

$$f_s(x_b)=\frac{R(1-x_b\cot \left(x_b\right))}{x_b\cot (x_b)-1+\delta -\mathrm{i}\delta \frac{c_b}{c_w}{x}_b},$$

$$u_{\mathrm{rad}}(x)=f_s(x_b)u^{in}(0)\frac{{\mathrm{e}}^{ik|x|}}{|x|}.$$

$$1-\delta -x_M\cot (x_M)=0\text{or}\tan (x_M)=\frac{x_M}{1-\delta }.$$

$$f_s(x_M)=\mathrm{i}\frac{c_w}{c_b{x}_M}\text{with}\frac{c_w}{c_b{x}_M}\approx 70.$$

$$\left|f_s(x_0)\right|=\underset{x\in [0,0.1]}{\sup }\left|f_s(x)\right|.$$

$$f_s(x_b)\approx \widetilde{f_s}(x_b):=\frac{1}{{\left(\frac{x_M}{x_b}\right)}^2-1-\mathrm{i}\left(\frac{c_b}{c_w}\right)x_b}.$$

$$D_j:=\left\{x\in {\mathbb{R}}^3,|x-x_j|\le R\right\},1\le j\le N.$$

$$u=\left\{\begin{array}[]{ll}\left(\sum\limits_{j=1}^{N}{\mathcal{S}}_{D_{j}}^{k_{w}}[\phi_{j,w}]\right)+u^{in},&x\in{\mathbb{R}}^{3}\setminus\overline{D},\\ &\\ {\mathcal{S}}_{D_{j}}^{k_{b}}[\phi_{j,b}]&x\in D_{j},\quad 1\leq j\leq N.\\ \end{array}\right.$$

$${\mathcal{A}}_{D_{1},\ldots D_{N}}\begin{pmatrix}\psi_{1,b}\\ \psi_{1,w}\\ \vdots\\ \psi_{N,b}\\ \psi_{N,w}\\ \end{pmatrix}=\begin{pmatrix}u^{in}\Big{|}_{\partial D_{1}}\\ \frac{\partial u^{in}}{\partial\nu_{1}}\Big{|}_{\partial D_{1}}\\ \vdots\\ u^{in}\Big{|}_{\partial D_{N}}\\ \frac{\partial u^{in}}{\partial\nu_{N}}\Big{|}_{\partial D_{N}}\\ \end{pmatrix},$$

Full text

Mathematics Subject Classification (MSC2000): 35R30, 35C20.

Keywords: Minnaert resonance, system of bubbles, sub-wavelength focusing.

1 Introduction

In this paper we are concerned with acoustic wave propagation in bubbly media. In particular, we are interested in the effect on the imaginary part of the Green’s function of exciting bubbles at frequencies near the Minnaert resonance frequency. It is well known that the focal spot size which limits the resolution limit in imaging is dependent on the imaginary part of the Green’s function due to the Helmholtz-Kirchoff Theorem [10]. Physical experiments have shown that it is possible to focus waves at the sub-wavelength scale by exploiting the strong scattering property of acoustic bubbles at their Minnaert resonance [14].

The Minnaert resonance is a quasi-static resonance in which the bubble size is much smaller than the wavelength of the incident wave. We use layer potential techniques to explicitly determine the acoustic surface potentials for a single bubble in the quasi-static regime, in which sound propagation can be neglected and the pressure field on the surface of the bubble may be considered a constant [4]. Knowledge of the acoustic potential corresponding to the exterior domain enables us to the derive the scattering function for a single bubble, which has a similar expression to the scattering function found in the physical literature. The Minnaert resonance corresponds to the value at which the scattering function is maximized.

In [3], a rigorous mathematical justification of the Minnaert resonance in the case of a single bubble in a homogeneous medium is established. The acoustic properties of the bubble are analyzed, and the approximate formula for the Minnaert resonance of an arbitrary shaped bubble is derived. In [1] the opening of a sub-wavelength phononic bandgap is demonstrated by considering a periodic arrangement of bubbles and exploiting their Minnaert resonance. As shown in [7], around the Minnaert resonant frequency, an effective medium theory can be derived in the dilute regime. Furthermore, above the Minnaert resonant frequency, the real part of the effective modulus is negative and consequently, the bubbly fluid behaves as a diffusive media for the acoustic waves. Meanwhile, below the Minnaert resonant frequency, with an appropriate bubble volume fraction, a high contrast effective medium can be obtained, making the sub-wavelength focusing or super-focusing of waves achievable [9]. These properties show that the bubbly fluid functions like an acoustic metamaterial and indicate that a sub-wavelength bandgap opening occurs at the Minneaert resonant frequency [15] . We remark that such behavior is rather analogous to the coupling of electromagnetic waves with plasmonic nanoparticles, which results in effective negative or high contrast dielectric constants for frequencies near the plasmonic resonance frequencies [2, 5, 6].

By generalizing the approach of [3] to the case of $N$ bubbles all having the same radius, we obtain an $N\times N$ block matrix in which the diagonal terms represent bubble self-interaction, while the off-diagonal terms represent the interaction between different bubbles. Letting the size of the bubbles go to zero we derive the point scatterer approximation which reduces the problem to an $N\times N$ linear system, which is recognizable as the classical Foldy-Lax formulation for a scattered field due to a collection of point scatterers [13].

We numerically simulate the system of point scatterers to analyze the super-focusing phenomenon observed in [14] and to simulate a time reversal experiment that demonstrates that it is possible to spatially localize a time reversed signal to a very high degree of accuracy in a bubbly media.

Due to the Helmholtz-Kirchhoff Theorem, the sharper the imaginary part of the Green’s function, the smaller the focal size and the higher the resolution achieved [10, 8, 9]. To better understand the effect of the bubbly media on the imaginary part of the Greens function, we excite the point scatterers over a wide range of frequencies. By averaging the imaginary part of the Green’s function function over a range of frequencies, we demonstrate that sub-wavelength focusing is achievable in the region slightly below the Minnaert resonance frequency.

The paper is organized as follows. In Section 2 we formulate the scattering problem for a single bubble and derive the scattering function. In Section 3 we derive the point scatterer approximation for a system of $N$ bubbles. In Section 4 we demonstrate the super-focusing phenomenon by numerically showing that the imaginary part of the Green’s function becomes sharper when the incident wave frequency is close to the Minnaert resonance frequency.

2 Single bubble

The purpose of this section is to justify the usual scattering function that one can find in the literature for bubbles. We explain the physical argument given in [13] using a layer potential approach [4]. In what follows we will only consider spherical bubbles. This is not a strong assumption for the phenomena we are interested in, and it already contains all the properties of a more complicated model; see [3]. We first consider a single bubble, which will be represented by a sphere $D$ of radius $R$, centered at the origin:

[TABLE]

We denote by $\rho_{b}$ and $\kappa_{b}$ the density and the bulk modulus of the air inside the bubble, respectively, and let $\rho_{w}$ and $\kappa_{w}$ be the corresponding parameters for the background media ${\mathbb{R}}^{3}\backslash D$. The scattering problem can be modeled by the following equations:

[TABLE]

We introduce four auxiliary parameters to facilitate our analysis:

[TABLE]

$c_{w}$ and $c_{b}$ represent the speed of sound in the background media and in the bubble, respectively. We also introduce the dimensionless contrast parameter

[TABLE]

and note that for water and air in regular conditions we have that $\delta\approx 1.2\times 10^{-3}$.

2.1 Layer potentials

We use layer potentials to represent the solution to the scattering problem (2.1). Let the single layer potential $\mathcal{S}_{D}^{k}$ associated with the domain $D$ and the wavenumber $k$ be defined by

[TABLE]

where

[TABLE]

is the Green function for the Helmholtz equation in ${\mathbb{R}}^{3}$, subject to the Sommerfeld radiation condition. We also define the boundary integral operator $\mathcal{K}_{D}^{k,*}$ by

[TABLE]

Then the solution $u$ can be written as

[TABLE]

for some surface potentials $\psi_{w},\psi_{b}\in L^{2}(\partial D)$. Using the jump relations for the single layer potentials, it is easy to derive that $\psi_{w}$ and $\psi_{b}$ satisfy the following system of boundary integral equations:

[TABLE]

where

[TABLE]

One can show that the scattering problem (2.1) is equivalent to the boundary integral equations (2.4).

2.2 The breathing approximation

2.2.1 The matrix $A_{D}$

Note that on the surface of the bubble, we have $u^{in}(x)=u^{in}(0)+O(k_{w}R)$, and $\frac{\partial u^{in}}{\partial\nu}=O(k_{w}R)$. If the wavelength is much larger than the size of the bubble, i.e., $k_{w}R\ll 1$, then we may approximate $u^{in}$ by $u^{in}(0)$ and $\frac{\partial u^{in}}{\partial\nu}$ by [math] on the surface of the bubble, respectively. In this case, we may look for solutions to (2.4) with the right-hand side given by $\begin{pmatrix}{\mathds{1}}_{\partial D}\\ 0\end{pmatrix}.$ This is the so-called breathing approximation [13].

Thanks to the spherical symmetry of $D$, the operator $\mathcal{A}_{D}$ can be diagonalized by spherical harmonic functions. Denote by $X=L^{2}(\partial D)\times L^{2}(\partial D)$. We define the finite-dimensional subspace $X_{1}$ of $X$ by

[TABLE]

The vector space $X_{1}$ is of dimension $2$. Using the breathing approximation, we only need to solve the following equation

[TABLE]

for $\Psi\in X_{1}$.

Let $A_{D}$ be the $2\times 2$ matrix representing the operator ${\mathcal{A}}_{D}$ on $X_{1}$. The values of the components of $A_{D}$ are deduced from the next lemma (see Section 5.1 for the proof).

Lemma 2.1**.**

It holds that

[TABLE]

and

[TABLE]

By introducing the dimensionless parameters $x_{b}=k_{b}R$ and $x_{w}=k_{w}R$, we obtain

[TABLE]

Denote $\psi_{b/w}=C_{b/w}{\mathds{1}}_{\partial D}$ where $(C_{b},C_{w})\in{\mathbb{C}}^{2}$. Then we have

[TABLE]

The determinant of $A_{D}$ is

[TABLE]

and the inverse of $A_{D}$ is

[TABLE]

Together with (2.5), we obtain that (we assume that $\sin(x_{b})\neq 0$ for simplicity)

[TABLE]

2.2.2 Scattering function

The pressure radiated by the bubble is ${u_{\rm rad}}:={\mathcal{S}}_{D}^{k_{w}}[\psi_{w}]$. Since $\psi_{w}=C_{w}{\mathds{1}}_{\partial D}$, we obtain ${u_{\rm rad}}=C_{w}{\mathcal{S}}_{D}^{k_{w}}[{\mathds{1}}_{\partial D}]$. The value of this last function is given by the next Lemma (see also (5.5)).

Lemma 2.2**.**

It holds that

[TABLE]

Together with (2.6), we deduce that the value of ${u_{\rm rad}}$ on the boundary $\partial D$ is (we assume that $\sin(x_{w})\neq 0$)

[TABLE]

Typical orders of magnitude are $\delta\approx 10^{-3}$ and $\delta\dfrac{c_{b}}{c_{w}}\approx 2\times 10^{-4}$. The dimensionless function

[TABLE]

is called the scattering function. It links the value of the incoming pressure on the boundary of the bubble to the value of the radiated pressure on the same boundary. Using the scattering function, we have the following monopole approximation for the scattered field

[TABLE]

Let us study the function $f_{s}$. The Minnaert resonance corresponds to the point $x_{b}=x_{M}$ such that

[TABLE]

This point corresponds to the frequency $x_{M}\approx 0.06$. At this particular frequency, it holds that

[TABLE]

Hence, at the Minnaert frequency, the pressure gains a phase of ${\mathrm{e}}^{{\mathrm{i}}\pi/2}$, and the amplitude is multiplied by a factor of $70$. For instance, for a bubble of radius $R=10^{-3}\ m$, we obtain a frequency $\omega_{M}/(2\pi)=\dfrac{c_{b}x_{M}}{2\pi R}\approx 3300\ \text{Hz}$, which is a perfectly audible sound. Note that $x_{M}$ does not correspond to any extremal point of $f_{s}$. This is however a vanishing point for the real part. We can also consider the frequency corresponding to $x=x_{0}$ such that

[TABLE]

Actually, it holds that $\left|x_{M}-x_{0}\right|\approx 6.10^{-4}$, so that $x_{M}$ is indeed a very good (and easily computable) approximation of $x_{0}$.

In the physical literature, the scattering function is usually approximated by the simpler function

[TABLE]

The functions $f_{s}(\omega)$ and $\tilde{f}_{s}(\omega)$ are shown in Figure 1. Numerically, we find that $\sup_{x\in[0,0.1]}\left|f_{s}-\widetilde{f_{s}}\right|\approx 0.47$, and hence $\widetilde{f_{s}}$ can be considered a highly accurate approximation of $f_{s}$.

3 System of N bubbles

We now consider a set of $N$ disjoint bubbles. In order to keep the notation simple, we will assume that all the bubbles have the same radius $R$. This is not a very restrictive assumption, and this case contains effects we are interested in. The volume occupied by bubble number $j$ is

[TABLE]

The distance between the centers of bubble $i$ and bubble $j$ is $d_{ij}=|x_{i}-x_{j}|$. The non-intersecting condition implies that $d_{ij}\geq 2R$ for all $i\neq j$. We finally denote by $D:=\bigcup_{1\leq j\leq N}D_{j}$ the volume occupied by all the bubbles.

We want to solve the Helmholtz equation $(\Delta+k^{2})u=0$ subject to boundary conditions analogous to the ones in (2.1). We make the ansatz

[TABLE]

By introducing $X:=\Pi_{j=1}^{N}L^{2}(\partial D_{j})\times L^{2}(\partial D_{j})$ , this leads to a system of equations of the form

[TABLE]

where ${\mathcal{A}}_{D_{1},\ldots,D_{N}}$ is acting on $X$ and where $\Phi:=(\psi_{1,b},\psi_{1,w},\ldots,\psi_{1,b},\psi_{1,w})^{T}\in X$ with $T$ being the transpose. The operator ${\mathcal{A}}_{D_{1},\ldots,D_{N}}$ has the block diagonal form

[TABLE]

where (see also (2.4))

[TABLE]

and, for $i\neq j$,

[TABLE]

Here, we introduced the operators ${\mathcal{S}}_{D_{i},D_{j}}^{k}:L^{2}(\partial D_{j})\to L^{2}(\partial D_{i})$ and ${\mathcal{L}}^{k}_{D_{i},D_{j}}:L^{2}(\partial D_{j})\to L^{2}(\partial D_{i})$, which describe the effect of bubble $j$ on bubble $i$. These operators are respectively defined by

[TABLE]

Such a system is usually called a system of $N$ scatterers for the Helmholtz equation, and has been the starting point for a large number of articles.

3.1 The point scatterer approximation for N-bubbles

We derive the point scatterer approximation for the N-bubble system. For this, we assume the breathing approximation for each of the bubbles, i.e., the pressure field applied to each bubble is constant. More precisely, for each $1\leq i\leq N$, we define $u^{in}_{i}$ to be the total field incident on the $i$-th bubble $D_{i}$, and $u^{s}_{i}$ to be the field scattered from $D_{i}$. It is clear that

[TABLE]

Due to our assumption, $u^{in}_{i}$ can be approximated by a constant function across the interface $\partial D_{i}$. This is a good approximation when the bubbles are well-separated, i.e., their distance is much greater than their size. Under this approximation, we can derive a good approximate solution for the system (3.2). Indeed, the constant $u^{in}_{i}(x_{i})$ can be used to approximate the pressure impinged on $D_{i}$. As a result, using (2.8), the scattered field $u^{s}_{i}$ can be approximated by

[TABLE]

By taking $x=x_{j}$ in the above equation and using (3.3), we obtain the following systems of linear equations for $u^{in}_{i}(x_{i})$:

[TABLE]

We define a matrix $M$ by setting

[TABLE]

Then the total field for the scattering problem by N bubbles can be approximated by

[TABLE]

where $(M^{-1})_{ij}$ denotes the $(ij)$-th component of the inverse matrix $M^{-1}$. This is called the point interaction approximation. See also [7] for a rigorous justification of this approximation.

4 Super focusing

We now perform a numerical investigation of the results that have been observed in [14] in order to assist us in developing a mathematical theory that explains the super-focusing phenomenon.

We take the density of water and air to be $\rho_{w}=1\times 10^{3}\ \text{kg m}^{-3}$, and $\rho_{b}=1.2\ \text{kg m}^{-3}$, respectively. We take the bulk modulus of water and air to be $\kappa_{b}=2.07\times 10^{9}\ \text{N m}^{-2}$, and $\kappa_{b}=127\times 10^{3}\ \text{N m}^{-2}$, respectively. This results in a speed of sound in water of $c_{w}\approx 1440\ \text{m s}^{-2}$.

We set the common radius of the bubbles to be $R=50\ \mu\ \text{m}$. In light of these parameters, the Minnaert resonant frequency is $\omega_{M}\approx 59\times 2\pi\ \text{kHz}\approx 390\ \text{kHz}$. We are interested in exciting the bubbles with a reference frequency $\omega_{R}$ which is slightly below the Minnaert resonance frequency. Hence, we choose $\omega_{R}=57.5\times 2\pi\ \text{kHz}$. We denote by by $\lambda_{w}=2\pi/k_{w}$ the wavelength in water. Typical values for the wavelength are $\lambda_{w}\approx 0.025\ m$. The scattering function $f_{s}$ is taken to be

[TABLE]

This expression is slightly different than (2.9), although it behaves in a similar manner. It is worth noting that the parameters we have chosen essentially mirror those used in the experiments described in [14].

Our experimental setup is also largely in agreement with [14]. We consider a source, located at the origin ${\mathbb{s}}=(0,0,0)^{T}$, that is surrounded by a cube containing uniformly distributed bubbles. The length of the cube is $L=0.01\ \text{m}\approx\lambda_{w}/2$. The gas volume fraction in the cube is $\Phi=2\times 10^{-4}$, which corresponds to a collection of approximately $380$ bubbles. We assume that the water outside this cube does not contain any bubbles. Finally, four receivers (time-reversal mirrors) are located at ${\mathbb{r}}_{1}=(0.02,0,0)^{T}$, ${\mathbb{r}}_{2}=(-0.02,0,0)^{T}$, ${\mathbb{r}}_{3}=(0,0.02,0)^{T}$ and ${\mathbb{r}}_{4}=(0,-0.02,0)^{T}$.

The experiment begins by sending a short pulse from the location ${\mathbb{r}}_{0}$. The signal $s(t)$ (shown in Figure 2) is given by

[TABLE]

To compute the forward problem, the signal is time-Fourier transformed with

[TABLE]

In practice, since we expect the recording to last much longer than the short impulse, the signal is recorded at each receiver for a time period of $[0,T]$, where $T=50\times 10^{-3}\ \text{s}$, with the numerical resolution of the signal being $\Delta t=1\times 10^{-6}\ \text{s}$. As a result, we are in possession of the Fourier transformed values for all $\omega\in[-\omega_{\rm max}/2,\omega_{\rm max}/2]$, with $\omega_{\rm max}\approx 3.1\times 10^{6}\ \text{Hz}$, and $\Delta\omega\approx 125\ \text{Hz}$. The plot of $\hat{s}(\omega)$ is given in Figure 3.

Due to the fact that the signal is real-valued, and oscillates near $\omega_{R}$, only the frequencies near $\omega_{R}$ are considered (in practice this mean all frequencies in the range $\omega\in[0,2\omega_{R}]$). This necessitates solving approximately 5750 Helmholtz problems. We use the point scatterer model described in the previous section to solve these problems, and obtain the frequency domain representations of the recorded signals $\hat{r}_{j}(\omega)$ for $1\leq j\leq 4$. Through application of the Fourier transform, we retrieve the recordings in the time domain $r_{j}(t)$ (see Figure 4). Note that while the initial signal is a short pulse, the recorded signal is much longer, which is due to the phenomenon of multiple echoing of waves among the bubbles. The non-null signal before $t_{0}$ is due to the Gibbs phenomenon.

We now simulate the time-reversal experiment. Each receiver now acts as a source, and emits the signal $r^{\sharp}_{j}(t)=r_{j}(T-t)$. Proceeding as in the forward problem, we time-Fourier transform the signal, and then solve the resulting Helmholtz problems. We record the time-reversed signal $s^{\sharp}(t)$ at the origin, and the result is shown in in Figure 5. In Figure 6 we plot the magnitudes of the normalized time-reversed pressure fields in the absence of bubbles, and in the presence of bubbles, respectively, for $x=(x_{1},0,0)^{\top}$ where $x_{1}\in[-L,L]$, and $t\in[0.035,0.05]$. While in both plots it is clear that the time-reversed signal arrives at the origin at the expected time, greatly enhanced spatial localization (super-focusing in space) is readily apparent in the presence of bubbles.

5 The imaginary part of the Green’s function

It is known in imaging that the level of focusing achievable is linked to the imaginary part of the Green’s function [10, 9]. In this section we concentrate on investigating the properties of the imaginary part of the Green’s function when a large number of small acoustic bubbles are introduced to the region under consideration.

As in the previous section, we consider a source at the origin surrounded by a cube of uniformly distributed bubbles, in which the gas volume fraction is $\Phi=2\times 10^{-4}$. We denote by $G_{\text{m}}(x,\omega)$ the Green’s function in the presence of the bubbles at frequency $\omega$. Using the point interaction approximation, we have

[TABLE]

where $f_{s}$ is the scattering function defined in (2.7), $c_{w}$ is the wave speed in water and

[TABLE]

We set $\omega_{\rm min}=15\times 2\pi\ \text{kHz}$ and $\omega_{\rm max}=155\times 2\pi\ \text{kHz}$. In Figure 7, we plot $|{\mathrm{Im}\,}(G(x,\omega))|$, and $|{\mathrm{Im}\,}(G_{\text{m}}(x,\omega))|$, for $x=(x_{1},0,0)^{\top}$ where $x\in[-L,L]$ and $\omega\in[\omega_{\rm min},\omega_{\rm max}]$ for three different random bubble configurations. We clearly observe that there are a lot of small oscillations with regard to the maps $\omega\mapsto{\mathrm{Im}\,}(G(x,\omega))$. This phenomenon was noticed and explained in the propagation of waves in random media, in the paraxial regime [11, 12]. We note that for each random configuration of bubbles, when the frequency is slightly less than the Minnaert resonance frequency (white dashed line), on average there is a very significant enhancement of the imaginary part of the effective Green’s function, as compared to the free-space Green’s function. Slightly above the Minnaert resonance frequency, on the other hand, the imaginary part of the Green’s function vanishes almost completely.

Finally, we consider the average of the imaginary part of the Green’s function over the frequency interval $[\omega_{\rm min},\omega_{\rm max}]$:

[TABLE]

In Figure 8 we compare the magnitude of this expression with the magnitude of the corresponding expression for the free-space Green’s function. We consider three separate cases. In the first case the averaging is performed over a range of frequencies just above the Minnaert resonance frequency, i.e., $\omega_{\rm min}=1.01\ \omega_{M}$ and $\omega_{\rm max}=2\ \omega_{M}$. In this case the presence of the bubbles can be seen to have a strong effect on propagation and $\langle\text{Im}(G_{\text{m}}(x))\rangle$ vanishes almost completely. Next we perform the averaging well away from and above the Minnaert resonance frequency, i.e., $\omega_{\rm min}=2\ \omega_{M}$ and $\omega_{\rm max}=4\ \omega_{M}$. It is clear that in this case the effect of the bubbles is greatly diminished and $\langle\text{Im}(G_{\text{m}}(x))\rangle$ appears almost the same as it would in a bubble free system.

Finally, we perform the averaging over a range of frequencies just below the Minnaert resonance frequency, i.e., we set $\omega_{\rm min}=0.8\ \omega_{M}$ and $\omega_{\rm max}=0.99\ \omega_{M}$. Here, we observe a strong peak near the origin, which is indicative of the super-focusing phenomenon. The full width half maximum (FWHM) of the free space Green’s function is $w\approx 0.0152$. In comparison, the FWHM of the Green’s function in the presence of bubbles is $w_{m}\approx 0.0016$, indicating that the presence of a large number of small bubbles, which have been excited near their Minnaert resonance frequency, leads to a substantial increase in focusing power.

5.1 Proof of Lemma 2.1

This is a classical result, which is a special case of a more general result. It is usually proven using the properties of spherical Bessel functions (see for instance [16]). We include here a simple proof for completeness.

Step 1: Proof for ${\mathcal{S}}_{D}^{k}$ with $x\in\partial D$.

Let us first prove the result for ${\mathcal{S}}_{D}^{k}$. We have that for $x\in\partial D$,

[TABLE]

Since the problem is invariant under rotation, it is sufficient to compute the integral for $x=(0,0,R)$. We use the spherical coordinates

[TABLE]

with $0\leq\theta\leq 2\pi$ and $0\leq\phi\leq\pi$. Then, for $y\in\partial D$, it holds that

[TABLE]

so that

[TABLE]

Together with the equality $\sin(\phi)=2\sin(\phi/2)\cos(\phi/2)$ and the change of variable $u=\sin(\phi/2)$, we end up with

[TABLE]

Step 2: Exact expression of ${\mathcal{S}}_{D}^{k}[{\mathds{1}}_{\partial D}]$ with $x\in\partial D$.

From the previous step, we can easily deduce an expression of ${\mathcal{S}}_{D}^{k}[{\mathds{1}}_{\partial D}]$ on ${\mathbb{R}}^{3}\setminus D$. Indeed, the function ${\mathcal{S}}_{D}^{k}[{\mathds{1}}_{\partial D}]$ outside $D$ is the (unique) solution to the Dirichlet equation

[TABLE]

together with the Sommerfeld radiation condition. The solution is therefore

[TABLE]

Step 3: Proof for $\left({\mathcal{K}}_{D}^{k}\right)^{*}$.

For all $x\in\partial D$,

[TABLE]

Together with (5.5), we obtain

[TABLE]