# On the detection of small moving disks in a fluid

Alexandre Munnier, Karim Ramdani

► To cite this version:

Alexandre Munnier, Karim Ramdani. On the detection of small moving disks in a fluid. SIAM Journal on Applied Mathematics, 2016, 76 (1), pp.159-177. 10.1137/141001226. hal-01098067v4

HAL Id: hal-01098067

https://inria.hal.science/hal-01098067v4

Submitted on 8 Mar 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# ON THE DETECTION OF SMALL MOVING DISKS IN A FLUID*

ALEXANDRE MUNNIER† AND KARIM RAMDANI‡

**Abstract.** We are interested in determining the positions and the velocities of moving rigid solids in a bounded cavity filled with a perfect fluid. We assume that the solids are small disks and that they move slowly. Using an integral formulation, we first derive the asymptotic expansion of the DtN map of the problem as the diameters of the disks tend to zero. Then, combining a suitable choice of exponential type data and the DORT method (which is usually used in inverse scattering for the detection of point-like scatterers), we propose a reconstruction method for the unknown positions and velocities.

**Key words.** inverse problems, perfect fluid, integral equations, asymptotic analysis, DtN operator, time reversal, DORT method

**AMS subject classifications.** 74F10, 45Q05, 35C20

**DOI.** 10.1137/141001226

**1. Introduction.** The geometric inverse problem that we consider in this paper is to determine the positions and the velocities of rigid solids moving in a perfect fluid. More precisely, let $\Omega$ be an open and simply connected bounded domain of $\mathbb{R}^2$ with smooth boundary $\Gamma := \partial\Omega$. The domain $\Omega$ is supposed to be filled with a perfect fluid and it contains $M$ rigid solids $D_m$, $m = 1, \dots, M$, where $D_m \subset \Omega$ is a closed disk with boundary $\gamma_m$. The domain occupied by the fluid is denoted by $\mathcal{F} := \Omega \setminus (\sum_{m=1}^M D_m)$. We also introduce the unit normal $n$ to $\partial\mathcal{F}$ directed towards the exterior of the fluid and $\tau$ the unit tangent vector to $\partial\mathcal{F}$ such that $\tau = n^{\perp}$ (throughout this paper, we set $x^{\perp} := (-x_2, x_1)$ for all $x = (x_1, x_2) \in \mathbb{R}^2$). We assume that every solid moves with a velocity $V_m \in \mathbb{R}^2$. At every given time, we can consider the Eulerian velocity field $\mathbf{U}(x)$, $x \in \mathcal{F}$, of the fluid. We assume that the flow is irrotational,

$$ (1.1a) \qquad \mathrm{rot}(\mathbf{U}) = 0 \quad \text{in } \mathcal{F}, $$

and circulation free, so that for every rectifiable Jordan curve $\gamma \subset \mathcal{F}$, there holds

$$ (1.1b) \qquad \int_{\gamma} \mathbf{U}(s) \cdot \tau(s) \, ds = 0. $$

The slip boundary conditions on $\partial\mathcal{F}$ read

$$ (1.2a) \qquad \mathbf{U} \cdot n = V_m \cdot n \quad \text{on } \gamma_m, \quad m = 1, \dots, M, $$

$$ (1.2b) \qquad \mathbf{U} \cdot n = 0 \quad \text{on } \Gamma. $$

We shall address the inverse problem of determining the positions and the velocities of the solids by using actuators and sensors located on the outer boundary $\Gamma$. We assume that we can prescribe the normal velocity of the fluid on $\Gamma$, leading to the modification of (1.2b) into

$$ (1.3) \qquad \mathbf{U} \cdot n = F \quad \text{on } \Gamma, $$

<footnote>*Received by the editors December 22, 2014; accepted for publication (in revised form) October 22, 2015; published electronically January 26, 2016.</footnote>

<footnote>http://www.siam.org/journals/siap/76-1/100122.html</footnote>

<footnote>†Université de Lorraine and CNRS, Institut Élie Cartan de Lorraine, UMR 7502, Vandœuvre-lès-Nancy, F-54506, France (alexandre.munnier@univ-lorraine.fr).</footnote>

<footnote>‡Inria, Université de Lorraine, Villers-lès-Nancy, F-54600, France (karim.ramdani@inria.fr).</footnote>
where $F$ is given. Regarding the available output, we assume that we can measure
the tangential velocity $\mathbf{U} \cdot \tau$ on $\Gamma$. According to conditions (1.1a) and (1.1b), we
classically introduce the stream function $\psi : \mathcal{F} \to \mathbb{R}$ such that $\mathbf{U} = -\nabla^\perp \psi$ in $\mathcal{F}$.
Then, (1.1a) and (1.2) read

$$
(1.4a) \qquad -\Delta\psi = 0 \qquad \text{in } \mathcal{F},
$$

$$
(1.4b) \qquad \psi = V_m^{\perp} \cdot x + c_m \quad \text{on } \gamma_m, m = 1, \dots, M,
$$

$$
(1.4c) \qquad \psi = f \qquad \text{on } \Gamma,
$$

where $\partial_{\tau} f = F$ on $\Gamma$ and the constants $c_m \in \mathbb{R}$, $k = 1, \dots, M$, are such that

$$
(1.4d) \qquad \int_{\gamma_m} \partial_n \psi(s) \, ds = 0.
$$

With these settings, the measurement reads $\mathbf{U} \cdot \tau = -\partial_n \psi$. The detection problem under consideration is then to recover the centers, the diameters, and the velocities of the solids from the DtN map $\Lambda : f \in H^{1/2}(\Gamma) \to \partial_n \psi \in H^{-1/2}(\Gamma)$.

Remark 1. Instead of using the stream function, one can equivalently formu-
late the problem in terms of the potential function and obtain Neumann boundary
conditions instead of Dirichlet in (1.4).

Although such geometric inverse problems of detecting moving solids in a fluid
appear in many applications, the associated literature is quite limited, as most con-
tributions deal with the case of motionless solids (i.e., obstacles). In particular, the
case of a fluid described by the Stokes equations has been investigated using opti-
mization methods by Caubet et al. in [12, 5, 11, 10] and more recently by Bourgeois
and Dardé in [6] using the quasi-reversibility method combined to a level set method.
The case of small obstacles has been studied in Caubet and Dambrine [9]. Regarding
moving obstacles, Conca et al. show in [15, 14] that the position and the velocity for
a single disk moving in a perfect fluid can be recovered from one measurement of the
velocity on part of the boundary. Linear stability estimates are also provided. In
Conca, Malik, and Munnier [16], the authors consider a moving rigid solid immersed
in a potential fluid and provide examples of detectable (ellipses for instance) and un-
detectable shapes. Conca, Schwindt, and Takahashi obtained in [17] an identifiability
result in the case of a rigid solid immersed in a viscous fluid.

<figure><img src="image_2.png" /><figcaption>FIG. 1. The domain $\Omega$ containing $M$ small rigid disks $D_m^\epsilon$, $m = 1, \dots, M$, and filled with a perfect fluid occupying the domain $\mathcal{F}^\epsilon = \Omega \setminus (\cup_{m=1}^M D_m^\epsilon)$.</figcaption></figure>

In this work, we restrict the analysis to the case where the solids are small disks
(see Figure 1) of typical size $\epsilon$, so that rotation plays no role. For such configurations,
we show an identifiability result and we provide a reconstruction method. More precisely, we denote each closed disk $D_m^\epsilon$, $m = 1, \dots, M$, and we denote by $\gamma_m^\epsilon$ its boundary. We assume that $D_m^\epsilon$ is centered at $r_m$ and is of radius $\epsilon R_m$, where the parameter $\epsilon$ is meant to tend to 0. Denoting by $\mathcal{F}^\epsilon = \Omega \setminus (\cup_{m=1}^M D_m^\epsilon)$ the domain occupied by the fluid and letting $\psi^\epsilon$ be the corresponding solution of (1.4), our inverse problem can be formulated as follows: *Knowing the DtN map*

$$ (1.5) \qquad \Lambda^{\epsilon} : f \in H^{1/2}(\Gamma) \mapsto \partial_{n} \psi^{\epsilon} \in H^{-1/2}(\Gamma) $$

is it possible, and if so how, to recover the number M of disks, their positions $r_m$, $m = 1, \dots, M$, their (rescaled) radii $R_m$, and their velocities $V_m$?

We answer this question in several steps. In section 2, we first derive the asymptotic expansion of the DtN operator $\Lambda^\epsilon$ as $\epsilon \to 0^+$ (Theorem 2.1) using a boundary integral formulation of the forward problem. In section 3, we combine this expansion with the DORT method<fnref n="1" /> to recover the number of disks and their positions, provided they are distant enough. Initially introduced by Fink and Prada [18], this method has been justified mathematically and used in the framework of wave systems for the detection of distant point-like scatterers in acoustics [21, 26, 7] and electromagnetics [4]. Once the positions have been determined, the velocities and rescaled radii can be easily recovered using suitable data. Finally, we collect in section 4 some numerical examples to illustrate the efficiency of the proposed reconstruction method.

**2. Asymptotic expansion for small disks.** Although the literature dealing with small inhomogeneities or small inclusions is quite rich (see, for instance, the review paper [2] and the book [3] by Ammari and Kang), the needed asymptotics is not, as far as we know, directly available for our problem (1.4), nor for the equivalent Neumann conjugate problem satisfied by the velocity potential. Quoting only the tightly related works, let us mention Il'in [22] who studied an elliptic boundary value problem set in a three dimensional (3D) domain containing a small hole, Maz'ya, Nazarov, and Plamenevskij who studied the Laplace problem with only one small hole (see [24, p. 59] for the Dirichlet case and [24, p. 291] for the Neumann case) and Friedman and Vogelius [19, Lemma 3.3] who obtained the first order asymptotics for several infinitely conducting small inhomogeneities (i.e., constant Dirichlet condition on the boundaries of the small holes and Neumann condition on the exterior boundary).

The next result provides the asymptotic expansion of the DtN map $\Lambda^\epsilon$ as $\epsilon \to 0^+$.

**THEOREM 2.1.** For every $f \in H^{1/2}(\Gamma)$, we denote by $U^f \in H^1(\Omega)$ the solution of the boundary value problem

$$ (2.1) \qquad \begin{array}{ll} -\Delta U^f = 0 & \text{in } \Omega, \\ U^f = f & \text{on } \Gamma. \end{array} $$

Let $\Lambda^0 \in \mathcal{L}(H^{1/2}(\Gamma), H^{-1/2}(\Gamma))$ the DtN map $\Lambda^0 : f \in H^{1/2}(\Gamma) \mapsto \partial_n U^f \in H^{-1/2}(\Gamma)$. Then, the DtN map $\Lambda^\epsilon \in \mathcal{L}(H^{1/2}(\Gamma), H^{-1/2}(\Gamma))$ defined by (1.5) admits the following expansion as $\epsilon \to 0^+$:

$$ \Lambda^{\epsilon} = \Lambda^{0} + \epsilon^{2}\Lambda^{2} + O(\epsilon^{3}), $$

where for every $f, g \in H^{1/2}(\Gamma)$,

$$ \langle \Lambda^2 f, g \rangle_{H^{-1/2}(\Gamma), H^{1/2}(\Gamma)} = 2\pi \sum_{\ell=1}^{M} R_{\ell}^2 \nabla U^f(r_{\ell}) \cdot \nabla U^g(r_{\ell}) + 2\pi \sum_{\ell=1}^{M} R_{\ell}^2 \nabla U^g(r_{\ell}) \cdot V_{\ell}^{\perp}. $$

<footnote n="1">¹DORT is the French acronym for Decomposition of the Time-Reversal Operator.</footnote>
The rest of this section is devoted to the proof of this result. The first ingredient is the following reciprocity identity, which can be easily obtained using Green's formula:

$$ (2.2) \qquad \langle (\Lambda^{\epsilon} - \Lambda^{0}) f, g \rangle_{H^{-1/2}(\Gamma), H^{1/2}(\Gamma)} = \sum_{\ell=1}^{M} [\psi^{\epsilon}, U^{g}]_{\gamma_{\ell}^{\epsilon}} $$

where we have set

$$ (2.3) \qquad [\phi_1, \phi_2]_{\gamma_{\ell}^{\epsilon}} := \langle \partial_n \phi_2, \phi_1 \rangle_{H^{-1/2}(\gamma_{\ell}^{\epsilon})} - \langle \partial_n \phi_1, \phi_2 \rangle_{H^{-1/2}(\gamma_{\ell}^{\epsilon})} $$

This formula shows in particular that the asymptotics of the bilinear form associated with $\Lambda^{\epsilon}$ can be obtained from the asymptotics of $\psi^{\epsilon}$ (which is given) and $\partial_n \psi^{\epsilon}$ on $\gamma_{\ell}^{\epsilon}$.

In order to obtain these asymptotics, we use a superposition result, usually referred to as Kirchhoff principle in the context of fluid dynamics. Thus, we have

$$ \psi^{\epsilon} = u^{\epsilon} + \sum_{m=1}^{M} V_{m}^{\perp} \cdot \begin{pmatrix} v_{m}^{\epsilon} \\ w_{m}^{\epsilon} \end{pmatrix}, $$

where $u^{\epsilon}$ satisfies

$$ \begin{array}{l} -\Delta u^{\epsilon} = 0 \quad \text{in } \mathcal{F}^{\epsilon}, \\ (2.4) \qquad u^{\epsilon} = a_{m}^{\epsilon} \quad \text{on } \gamma_{m}^{\epsilon}, (m = 1, \dots, M), \\ \qquad u^{\epsilon} = f \quad \text{on } \Gamma, \end{array} $$

and $v_m^\epsilon$ and $w_m^\epsilon$ solve for every $m = 1, \dots, M$ the following boundary value problems:

$$ \left| \begin{array}{ll} -\Delta v_m^\epsilon = 0 & \text{in } \mathcal{F}^\epsilon, \\ v_m^\epsilon = x_1 + b_{m,m}^\epsilon & \text{on } \gamma_m^\epsilon, \\ v_m^\epsilon = b_{m,k}^\epsilon & \text{on } \gamma_k^\epsilon, k \neq m, \\ v_m^\epsilon = 0 & \text{on } \Gamma, \end{array} \right. \qquad \left| \begin{array}{ll} -\Delta w_m^\epsilon = 0 & \text{in } \mathcal{F}^\epsilon, \\ w_m^\epsilon = x_2 + c_{m,m}^\epsilon & \text{on } \gamma_m^\epsilon, \\ w_m^\epsilon = c_{m,k}^\epsilon & \text{on } \gamma_k^\epsilon, k \neq m, \\ w_m^\epsilon = 0 & \text{on } \Gamma, \end{array} \right. $$

where the constants $a_m^\epsilon$, $b_{m,k}^\epsilon$, and $c_{m,k}^\epsilon$ are such that

$$ (2.5) \qquad \int_{\gamma_m^\epsilon} \partial_n u^\epsilon(s) \, ds = \int_{\gamma_m^\epsilon} \partial_n v^\epsilon(s) \, ds = \int_{\gamma_m^\epsilon} \partial_n w^\epsilon(s) \, ds = 0, \quad m = 1, \dots, M, $$

where these integrals should be understood in the sense of duality

$$ \int_{\gamma_m^\epsilon} \partial_n u^\epsilon(s) \, ds = \langle \partial_n u^\epsilon, 1 \rangle_{H^{-1/2}(\gamma_m^\epsilon), H^{1/2}(\gamma_m^\epsilon)}. $$

We detail in section 2.1 the derivation of the asymptotics of $u^\epsilon$, based on a boundary integral formulation of problem (2.4)–(2.5). Then, as the proofs are quite similar, we only sketch in section 2.2 the proofs for $v_m^\epsilon$ and $w_m^\epsilon$. These three asymptotic expansions are used to obtain the asymptotics of

$$ (2.6) \qquad [\psi^{\epsilon}, U^{g}]_{\gamma_{\ell}^{\epsilon}} = [u^{\epsilon}, U^{g}]_{\gamma_{\ell}^{\epsilon}} + \sum_{m=1}^{M} V_{m}^{\perp} \cdot \begin{pmatrix} [v_{m}^{\epsilon}, U^{g}]_{\gamma_{\ell}^{\epsilon}} \\ [w_{m}^{\epsilon}, U^{g}]_{\gamma_{\ell}^{\epsilon}} \end{pmatrix}, $$

and the result follows then from (2.2). For the sake of clarity and for reader's convenience, we preferred to give a constructive proof of these expansions.
**2.1. Asymptotic expansion for the function $u^\epsilon$.** The next result collects well-known properties of the single layer potential that we use in what follows (see, e.g., McLean [25], Rjasanow and Steinbach [27], or Steinbach [28, Chapter 6]).

**PROPOSITION 2.2.** For $x \in \mathbb{R}^2 \setminus \{0\}$, let $G(x) := -\frac{1}{2\pi} \log|x|$ denote the Green function of the operator $-\Delta$ in $\mathbb{R}^2$. Let $\gamma$ denote the smooth boundary of a bounded (possibly multiply connected) domain of $\mathbb{R}^2$. We introduce the single layer potential associated to a given density $p$:

$$ (2.7) \qquad (S_\gamma p)(x) := \int_{\gamma} G(x-y) p(y) \, d\sigma_y, \quad x \in \mathbb{R}^2 \setminus \gamma. $$

Then, the following assertions hold true:

1. $S_\gamma$ defines a continuous operator from $H^{-1/2}(\gamma)$ onto $H^1(\mathbb{R}^2)$.

2. The trace of $S_\gamma p$ on $\gamma$ is given by the boundary integral operator

$$ (2.8) \qquad (S_\gamma p)(x) := \int_{\gamma} G(x-y) p(y) \, d\sigma_y, \quad x \in \gamma. $$

Moreover, $S_\gamma$ defines an isomorphism from $H^{-1/2}(\gamma)$ to $H^{1/2}(\gamma)$ and the quantity $\|\cdot\|_{-\frac{1}{2}} := \langle \cdot, S_\gamma \rangle_{H^{-1/2}(\gamma), H^{1/2}(\gamma)}^{1/2}$ defines on $H^{-1/2}(\gamma)$ a norm equivalent to the classical norm.

3. The normal derivative of $S_\gamma p$ on $\gamma$ is given by the so-called jump condition

$$ (2.9) \qquad p = \left[ \frac{\partial(S_\gamma p)}{\partial n} \right] := \left. \frac{\partial(S_\gamma p)}{\partial n} \right|_{-} - \left. \frac{\partial(S_\gamma p)}{\partial n} \right|_{+} $$

in which $n$ denotes the exterior unit normal to $\gamma$ (pointing from the interior bounded domain inside $\gamma$ to the unbounded exterior) and the $+$ and $-$ signs refer, respectively, to the normal derivatives coming from the exterior or the interior of $\gamma$.

4. Assume that $\gamma$ has only one connected component. Then, there exists a unique density $\psi^{\text{eq}} \in H^{-1/2}(\gamma)$, called the equilibrium density of $\gamma$, such that $S_\gamma \psi^{\text{eq}}$ is constant on $\gamma$ and satisfies the normalization condition $\int_{\gamma} \psi^{\text{eq}}(y) \, d\sigma_y = 1$.

We seek $u^\epsilon$ as a single layer potential:

$$ (2.10) \qquad u^\epsilon(x) = \sum_{m=1}^{M} S_{\gamma_m^\epsilon} p_m^\epsilon(x) + S_\Gamma q^\epsilon(x), \quad x \in \mathcal{F}^\epsilon, $$

where $p_m^\epsilon$ and $q^\epsilon$ are densities in $H^{-1/2}(\Gamma)$ and $H^{-1/2}(\gamma_m^\epsilon)$ to be determined.

**Remark 2.** Since by Green's formula $\int_{\gamma_\ell^\epsilon} \frac{\partial}{\partial n}(S_{\gamma_m^\epsilon} p_m^\epsilon)|_- = 0$, relation (2.9) shows that $\int_{\gamma_m^\epsilon} p_m^\epsilon(y) \, d\sigma_y = - \int_{\gamma_m^\epsilon} \partial_n u^\epsilon(y) \, d\sigma_y = 0$ for all $m = 1, \dots, M$.

Taking into account Remark 2, system (2.4)–(2.5) satisfied by $u^\epsilon$ is equivalent to the following system of $(M+1)$ coupled integral equations on $\gamma_1^\epsilon, \dots, \gamma_M^\epsilon$ and $\Gamma$:

$$ (2.11\text{a}) \qquad S_{\gamma_\ell^\epsilon} p_\ell^\epsilon + \sum_{\substack{m=1 \\ m \neq \ell}}^{M} (S_{\gamma_m^\epsilon} p_m^\epsilon)|_{\gamma_\ell^\epsilon} + (S_\Gamma q^\epsilon)|_{\gamma_\ell^\epsilon} = a_\ell^\epsilon \quad \text{on } \gamma_\ell^\epsilon, \ell = 1, \dots, M, $$

$$ (2.11\text{b}) \qquad \sum_{m=1}^{M} (S_{\gamma_m^\epsilon} p_m^\epsilon)|_{\Gamma} + S_\Gamma q^\epsilon = f \quad \text{on } \Gamma, $$

$$ (2.11\text{c}) \qquad \int_{\gamma_\ell^\epsilon} p_\ell^\epsilon(y) \, d\sigma_y = 0, \quad \ell = 1, \dots, M. $$
Defining the curve $\gamma_{\ell} = r_{\ell} + \epsilon^{-1}(\gamma_{\ell}^{\epsilon} - r_{\ell})$ and setting $r_{\ell m} = r_{\ell} - r_{m}$, a simple rescaling
leads us to rewrite (2.11a)–(2.11b) in the equivalent form:

$$
\begin{equation} \tag{2.12a}
\begin{split}
& G(\epsilon) \int_{\gamma_{\ell}} \epsilon p_{\ell}^{\epsilon}(r_{\ell} + \epsilon(y - r_{\ell})) \, d\sigma_y + \int_{\gamma_{\ell}} G(x-y) \epsilon p_{\ell}^{\epsilon}(r_{\ell} + \epsilon(y - r_{\ell})) \, d\sigma_y \\
& \qquad + \sum_{\substack{m=1 \\ m \neq \ell}}^{M} \int_{\gamma_m} G(r_{\ell m} + \epsilon(x-y-r_{\ell m})) \epsilon p_m^{\epsilon}(r_m + \epsilon(y-r_m)) \, d\sigma_y \\
& \qquad + \int_{\Gamma} G(r_{\ell} - y + \epsilon(x-r_{\ell})) q^{\epsilon}(y) \, d\sigma_y = a_{\ell}^{\epsilon}, \quad x \in \gamma_{\ell},
\end{split}
\end{equation}
$$

$$
\begin{equation} \tag{2.12b}
\begin{split}
& \sum_{m=1}^{M} \int_{\gamma_m} G(x - r_m - \epsilon(y - r_m)) \epsilon p_m^{\epsilon}(r_m + \epsilon(y - r_m)) \, d\sigma_y \\
& \qquad + \int_{\Gamma} G(x - y) q^{\epsilon}(y) \, d\sigma_y = f(x), \quad x \in \Gamma.
\end{split}
\end{equation}
$$

The circulation-free condition (2.11c), once the moving disks are rescaled, becomes

$$
(2.13) \quad \int_{\gamma_m} \epsilon p_m^{\epsilon} (r_m + \epsilon(y - r_m)) \, d\sigma_y = 0 \quad (m = 1, \dots, M).
$$

For every $\epsilon > 0$, we introduce the following a priori decomposition of the densities:

$$
(2.14a) \quad \epsilon p_{\ell}^{\epsilon}(r_{\ell} + \epsilon(y - r_{\ell})) = p_{\ell,0}(y) + \epsilon p_{\ell,1}(y) + \epsilon^2 \mathcal{P}_{\ell}^{\epsilon}(y), \quad y \in \gamma_{\ell},
$$

$$
(2.14b) \qquad q^{\epsilon}(y) = q_0(y) + \epsilon q_1(y) + \epsilon^2 q_2(y) + \epsilon^3 Q^{\epsilon}(y), \quad y \in \Gamma.
$$

Likewise, we write the constants $a_{\ell}^{\epsilon}$ as

$$
(2.15) \qquad a_{\ell}^{\epsilon} = a_{\ell,0} + \epsilon a_{\ell,1} + \epsilon^2 \mathcal{A}_{\ell}^{\epsilon},
$$

where the terms appearing in the right-hand sides of (2.14)–(2.15) have to be deter-
mined. Classically, this can be achieved by plugging the above a priori expansions
in (2.12) and identifying the terms of same order in the expansion. More precisely,
(2.12) read

$$
\begin{align*}
& (2.16\text{a}) \quad G(\epsilon) \int_{\gamma_{\ell}} [p_{\ell,0}(y) + \epsilon p_{\ell,1}(y) + \epsilon^2 \mathcal{P}_{\ell}^{\epsilon}(y)] \, d\sigma_y \\
& \qquad + \int_{\gamma_{\ell}} G(x-y) [p_{\ell,0}(y) + \epsilon p_{\ell,1}(y) + \epsilon^2 \mathcal{P}_{\ell}^{\epsilon}(y)] \, d\sigma_y \\
& \quad + \sum_{\substack{m=1 \\ m \neq \ell}}^{M} \int_{\gamma_m} G(r_{\ell m} + \epsilon(x-y-r_{\ell m})) [p_{m,0}(y) + \epsilon p_{m,1}(y) + \epsilon^2 \mathcal{P}_m^{\epsilon}(y)] \, d\sigma_y \\
& \qquad + \int_{\Gamma} G(r_{\ell} - y + \epsilon(x-r_{\ell})) [q_0(y) + \epsilon q_1(y) + \epsilon^2 q_2(y) + \epsilon^3 Q^{\epsilon}(y)] \, d\sigma_y \\
& \qquad = a_{\ell,0} + \epsilon a_{\ell,1} + \epsilon^2 \mathcal{A}_{\ell}^{\epsilon}, \quad x \in \gamma_{\ell}, \ell = 1, \dots, M,
\end{align*}
$$
$$ (2.16\text{b}) \quad \sum_{m=1}^{M} \int_{\gamma_m} G(x-r_m - \epsilon(y-r_m)) [p_{m,0}(y) + \epsilon p_{m,1}(y) + \epsilon^2 \mathcal{P}_m^\epsilon(y)] \, d\sigma_y \\ \qquad + \int_{\Gamma} G(x-y) [q_0(y) + \epsilon q_1(y) + \epsilon^2 q_2(y) + \epsilon^3 Q^\epsilon(y)] \, d\sigma_y = f(x), \quad x \in \Gamma, $$

while the circulation free condition (2.13) implies the following conditions:

$$ (2.17) \quad \int_{\gamma_{\ell}} p_{\ell,0}(y) \, d\sigma_y = \int_{\gamma_{\ell}} p_{\ell,1}(y) \, d\sigma_y = \int_{\gamma_{\ell}} \mathcal{P}_{\ell}^{\epsilon}(y) \, d\sigma_y = 0, \quad \ell = 1, \dots, M. $$

For every $X, Y \in \mathbb{R}^2$, $X \neq 0$, we have

$$ G(X + \epsilon Y) = G(X) + \epsilon \eta_1^\epsilon(X, Y) = G(X) + \epsilon \nabla G(X) \cdot Y + \epsilon^2 \eta_2^\epsilon(X, Y), $$

where

$$ \eta_1^\epsilon(X, Y) = \int_0^1 \nabla G(X + \epsilon sY) \cdot Y \, ds, $$

$$ \eta_2^\epsilon(X, Y) = \int_0^1 (1-s) \langle D^2 G(X + \epsilon sY), Y, Y \rangle \, ds. $$

Using the above Green's function asymptotics in (2.16) and identifying the terms of order 0 in $\epsilon$ yields

$$ (2.18\text{a}) \quad S_{\gamma_{\ell}} p_{\ell,0}(x) + S_{\Gamma} q_0(r_{\ell}) = a_{\ell,0}, \quad x \in \gamma_{\ell}, \ell = 1, \dots, M, $$

$$ (2.18\text{b}) \quad S_{\Gamma} q_0(x) = f(x), \quad x \in \Gamma. $$

The last equation uniquely determines $q_0 = S_{\Gamma}^{-1} f \in H^{-1/2}(\Gamma)$. In particular, recalling the notation $U^f$ introduced in (2.1), we have $S_{\Gamma} q_0 = U^f$ in $\Omega$. Choosing $p_{\ell,0} = 0$ and $a_{\ell,0} = U^f(r_{\ell})$ for every $\ell = 1, \dots, M$, the first equation in (2.18) is then fulfilled.

Now, we identify the terms of order 1 in $\epsilon$ of system (2.12), and get

$$ (2.19\text{a}) \quad S_{\gamma_{\ell}} p_{\ell,1}(x) + S_{\Gamma} q_1(r_{\ell}) + \nabla U^f(r_{\ell}) \cdot (x - r_{\ell}) = a_{\ell,1}, \quad x \in \gamma_{\ell}, \ell = 1, \dots, M $$

$$ (2.19\text{b}) \quad S_{\Gamma} q_1(x) = 0, \quad x \in \Gamma. $$

From the last equation, we have $q_1 = 0$. Taking the scalar product of the first equation of (2.19) with the equilibrium density $\psi_{\ell}^{\text{eq}}$ of $\gamma_{\ell}$ (see Proposition 2.2), we get that $\nabla U^f(r_{\ell}) \cdot (\int_{\gamma_{\ell}} (x-r_{\ell})\psi_{\ell}^{\text{eq}}(x) d\sigma_x) = a_{\ell,1}$. Since the moment of order 1 of the equilibrium density of a circle coincides with its center (this follows in particular from relation (B.4) in Lemma B.1 of Appendix B, with $V = 0$ and $c = (2\pi R_{\ell})^{-1}$), we have $a_{\ell,1} = 0$. The first equation of (2.19) then yields $S_{\gamma_{\ell}} p_{\ell,1}(x) = \nabla U^f(r_{\ell}) \cdot (r_{\ell} - x)$, and according to (B.4) in Lemma B.1, we thus have $p_{\ell,1}(x) = \frac{2}{R_{\ell}} \nabla U^f(r_{\ell}) \cdot (r_{\ell} - x)$.

Let us now consider the terms of order 2 in (2.12b) (this will be enough to derive the expected asymptotic expansion of $\Lambda^{\epsilon}$). We easily obtain that

$$ (2.20) \quad S_{\Gamma} q_2(x) = \sum_{m=1}^{M} \nabla G(x - r_m) \cdot \left( \int_{\gamma_m} (y - r_m) p_{m,1}(y) \, d\sigma_y \right), \quad x \in \Gamma. $$

Applying once again Lemma B.1, with $V = -\nabla U^f(r_{\ell})$ and $c = \nabla U^f(r_{\ell}) \cdot r_{\ell}$, it follows easily from (B.4) that for every $m \in \{1, \dots, M\}$,

$$ (2.21) \quad \int_{\gamma_m} (y - r_{\ell}) p_{m,1}(y) \, d\sigma_y = -2\pi R_m^2 \nabla U^f(r_{\ell}). $$
Consequently, (2.20) shows that $q_2$ solves the boundary integral equation

$$ (2.22) \qquad S_{\Gamma} q_2(x) = - \sum_{m=1}^{M} 2\pi R_m^2 \nabla U^f(r_m) \cdot \nabla G(x-r_m), \quad x \in \Gamma. $$

Let us focus now on the remainders $P_{\ell}^{\epsilon} \in H^{-1/2}(\gamma_{\ell})$, $Q^{\epsilon} \in H^{-1/2}(\Gamma)$, and $A^{\epsilon} \in \mathbb{R}$ appearing in (2.14) and (2.15), and prove that they are bounded. This will provide a justification of the formal a priori expansion of $p_{\ell}^{\epsilon}$, $q^{\epsilon}$ and $a_{\ell}^{\epsilon}$. To do so, we need the following general result, proved in Appendix A.

**LEMMA 2.3.** Let $\Gamma_1$ and $\Gamma_2$ be two rectifiable Jordan curves contained in two open, bounded, and nonintersecting sets $V_1$ and $V_2$. Let $K: (x, y) \in V_1 \times V_2 \mapsto K(x, y) \in \mathbb{R}$ be a function of class $C^{\infty}$. Then there exists a constant $C > 0$ depending only on $V_1$, $V_2$, and $\|K\|_{H^2(V_1 \times V_2)}$ such that for all $p \in H^{-1/2}(\Gamma_1)$ and $q \in H^{-1/2}(\Gamma_2)$,

$$ \left| \int_{\Gamma_1} \int_{\Gamma_2} K(x, y) p(x) q(y) \, d\sigma_y d\sigma_x \right| \le C \|p\|_{H^{-1/2}(\Gamma_1)} \|q\|_{H^{-1/2}(\Gamma_2)}. $$

We are now in position to prove the following result.

**LEMMA 2.4.** The remainders $P_{\ell}^{\epsilon} \in H^{-1/2}(\gamma_{\ell})$, $Q^{\epsilon} \in H^{-1/2}(\Gamma)$, and $A^{\epsilon} \in \mathbb{R}$ appearing in (2.14) and (2.15) are uniformly bounded in $\epsilon$.

*Proof.* On $\Gamma$, we have

$$ S_{\Gamma} Q^{\epsilon}(x) = - \sum_{m=1}^{M} \int_{\gamma_m} \eta_{1}^{\epsilon}(x-r_m, r_m-y) P_{m}^{\epsilon}(y) d\sigma_y + \int_{\gamma_m} \eta_{2}^{\epsilon}(x-r_m, r_m-y) p_{m,1}(y) d\sigma_y. $$

Multiplying by $Q^{\epsilon}$, we get (recall that $\|Q^{\epsilon}\|_{-\frac{1}{2}}^2 = \langle Q^{\epsilon}, S_{\Gamma} Q^{\epsilon} \rangle_{H^{-1/2}(\Gamma), H^{1/2}(\Gamma)}$)

$$ \|Q^{\epsilon}\|_{-\frac{1}{2}}^2 = - \sum_{m=1}^{M} \int_{\Gamma} \int_{\gamma_m} \eta_{1}^{\epsilon}(x-r_m, r_m-y) P_{m}^{\epsilon}(y) Q^{\epsilon}(x) d\sigma_y d\sigma_x \\ + \int_{\Gamma} \int_{\gamma_m} \eta_{2}^{\epsilon}(x-r_m, r_m-y) p_{m,1}(y) Q^{\epsilon}(x) d\sigma_y d\sigma_x. $$

Lemma 2.3 implies the existence of a constant $C > 0$ independent of $\epsilon > 0$ such that

$$ (2.23) \qquad \|Q^{\epsilon}\|_{-\frac{1}{2}} \le C \left( 1 + \sum_{m=1}^{M} \|P_{m}^{\epsilon}\|_{-\frac{1}{2}} \right). $$

On $\gamma_{\ell}$, we have

$$ (2.24) \quad S_{\gamma_{\ell}} P_{\ell}^{\epsilon}(x) + \sum_{\substack{m=1 \\ m \neq \ell}}^{M} \int_{\gamma_m} \eta_1^{\epsilon}(r_{\ell m}, x-y-r_{\ell m}) p_{m,1}(y) \, d\sigma_y \\ + \epsilon \int_{\gamma_m} \eta_1^{\epsilon}(r_{\ell m}, x-y-r_{\ell m}) P_m^{\epsilon}(y) \, d\sigma_y + \int_{\Gamma} \eta_2^{\epsilon}(r_{\ell} - y, x-r_{\ell}) q_0(y) \, d\sigma_y \\ + \int_{\Gamma} G(r_{\ell} - y + \epsilon(x-r_{\ell})) q_2(y) \, d\sigma_y + \epsilon \int_{\Gamma} G(r_{\ell} - y + \epsilon(x-r_{\ell})) Q^{\epsilon}(y) \, d\sigma_y = A_{\ell}^{\epsilon}. $$
Multiplying by the equilibrium density $\psi_{\ell}^{\text{eq}}$ of $\gamma_{\ell}$, integrating on $\gamma_{\ell}$, and applying Lemma 2.3, we get, since $\int_{\gamma} \mathcal{P}_{\ell}^{\varepsilon}(y) d\sigma_{y} = 0$,

$$ (2.25) \qquad |\mathcal{A}_{\ell}^{\varepsilon}| \le C \left( 1 + \varepsilon \sum_{\substack{m=1 \\ m \ne \ell}}^{M} \|\mathcal{P}_{m}^{\varepsilon}\|_{-\frac{1}{2}} + \varepsilon \|\mathcal{Q}^{\varepsilon}\|_{-\frac{1}{2}} \right), $$

where the constant $C > 0$ does not depend on $\varepsilon$. Multiplying now (2.24) by $\mathcal{P}_{\ell}^{\varepsilon}$, integrating on $\gamma_{\ell}$, and using once again Lemma 2.3, we get

$$ (2.26) \qquad \|\mathcal{P}_{\ell}^{\varepsilon}\|_{-\frac{1}{2}} \le C \left( 1 + \varepsilon \sum_{\substack{m=1 \\ m \ne \ell}}^{M} \|\mathcal{P}_{m}^{\varepsilon}\|_{-\frac{1}{2}} + \varepsilon \|\mathcal{Q}^{\varepsilon}\|_{-\frac{1}{2}} \right). $$

Combining the estimates (2.23) and (2.26), we immediately obtain that

$$ \|\mathcal{P}_{\ell}^{\varepsilon}\|_{-\frac{1}{2}} \le C \left( 1 + \varepsilon \sum_{\substack{m=1 \\ m \ne \ell}}^{M} \|\mathcal{P}_{m}^{\varepsilon}\|_{-\frac{1}{2}} \right), $$

which clearly implies that $\|\mathcal{P}_{\ell}^{\varepsilon}\|_{-\frac{1}{2}}$ is uniformly bounded, for every $\ell = 1, \dots, M$. The conclusion of the lemma follows then from (2.23) and (2.25). $\square$

Summing up, we have proved the following result.

**PROPOSITION 2.5.** With the above notation, we have for every $\ell = 1, \dots, M$,

$$ (2.27a) \qquad \varepsilon p_{\ell}^{\varepsilon}(r_{\ell} + \varepsilon(y - r_{\ell})) = \varepsilon p_{\ell,1}(y) + \varepsilon^2 \mathcal{P}_{\ell}^{\varepsilon}(y), \qquad y \in \gamma_{\ell}, $$

$$ (2.27b) \qquad q^{\varepsilon}(y) = q_{0}(y) + \varepsilon^{2}q_{2}(y) + \varepsilon^{3}Q^{\varepsilon}(y), \qquad y \in \Gamma, $$

$$ (2.27c) \qquad a_{\ell}^{\varepsilon} = a_{\ell,0} + \varepsilon^2 \mathcal{A}_{\ell}^{\varepsilon}, $$

where

$$ a_{\ell,0} = U^{f}(r_{\ell}), \qquad p_{\ell,1} = \frac{2}{R_{\ell}} \nabla U^{f}(r_{\ell}) \cdot (r_{\ell} - x), $$

$$ q_0 = S_{\Gamma}^{-1} f, \qquad q_2 = -S_{\Gamma}^{-1} \left( \sum_{m=1}^{M} 2\pi R_m^2 \nabla U^f(r_m) \cdot \nabla G(-r_m) \right), $$

and where $\|\mathcal{P}_{\ell}^{\varepsilon}\|_{H^{-1/2}(\gamma_{\ell})}$, $\|\mathcal{Q}^{\varepsilon}\|_{H^{-1/2}(\Gamma)}$, and $\mathcal{A}_{\ell}^{\varepsilon}$ are uniformly bounded in $\varepsilon$.

Now, we can compute the contribution $[u^{\varepsilon}, U^{g}]_{\gamma_{\ell}^{\varepsilon}}$ to the expansion of $\Lambda^{\varepsilon}$ (see (2.2) and (2.6)). On $\gamma_{\ell}^{\varepsilon}$, we have on the one hand $u_{\ell}^{\varepsilon} = a_{\ell}^{\varepsilon}$ and, on the other hand, $p_{\ell}^{\varepsilon} = [\partial_{n}u^{\varepsilon}] = \partial_{n}u^{\varepsilon}$ (the first equality follows from the jump relation (2.9) and the second one from the fact that the continuous harmonic extension of $u^{\varepsilon}$ inside $D_{\ell}^{\varepsilon}$ is the constant function $a_{\ell}^{\varepsilon}$). Using the asymptotic expansions of Proposition 2.5, simple computations lead to the following result.

**THEOREM 2.6.** Let $u^{\varepsilon} \in H^1(\mathcal{F}^{\varepsilon})$ be the solution of (2.4)-(2.5). Using the notation (2.3), we have

$$ [u^{\varepsilon}, U^{g}]_{\gamma_{\ell}^{\varepsilon}} = \varepsilon^2 (2\pi R_{\ell}^2 \nabla U^{f}(r_{\ell}) \cdot \nabla U^{g}(r_{\ell})) + O(\varepsilon^3). $$
## 2.2. Asymptotic expansion for the functions $v_m^\epsilon$ and $w_m^\epsilon$.
Let us deal with the function $v_1^\epsilon$, as the other cases ($v_m^\epsilon, m \ge 2$ and $w_m^\epsilon, m = 1, \dots, M$) can be treated in a similar way. For the sake of simplicity, we drop the subscript referring to this disk's number in the proof and we merely denote this function by $v^\epsilon$. Similarly, we rename $D^\epsilon$ the disk $D_1^\epsilon$, and $\gamma^\epsilon$ its boundary, while the other disks are renamed $D_m^\epsilon, m = 1, \dots, N := M-1$, with boundaries $\gamma_m^\epsilon$. With this notation, $v^\epsilon$ solves

$$ (2.28a) \qquad -\Delta v^{\epsilon} = 0 \qquad \text{in } \mathcal{F}^{\epsilon}, $$

$$ (2.28b) \qquad v^{\epsilon} = x_1 + b^{\epsilon} \qquad \text{on } \gamma^{\epsilon}, $$

$$ (2.28c) \qquad v^{\epsilon} = b_{m}^{\epsilon} \qquad \text{on } \gamma_{m}^{\epsilon}, m = 1, \dots, N, $$

$$ (2.28d) \qquad v^{\epsilon} = 0 \qquad \text{on } \Gamma, $$

where the constants $b^\epsilon$ and $b_m^\epsilon$, $m = 1, \dots, N$, are such that

$$ (2.28e) \qquad \int_{\gamma^{\epsilon}} \partial_n v^{\epsilon} \, d\sigma_y = 0, \qquad \int_{\gamma_m^{\epsilon}} \partial_n v^{\epsilon} \, d\sigma_y = 0. $$

We seek $v^\epsilon$ as a single layer potential: $v^\epsilon(x) = S_{\gamma^\epsilon} p^\epsilon(x) + \sum_{m=1}^N S_{\gamma_m^\epsilon} p_m^\epsilon(x) + S_\Gamma q^\epsilon(x)$, for $x \in \mathcal{F}^\epsilon$, where $p^\epsilon \in H^{-1/2}(\gamma^\epsilon)$, $p_m^\epsilon \in H^{-1/2}(\gamma_m^\epsilon)$, and $q^\epsilon \in H^{-1/2}(\Gamma)$ are to be determined. Problem (2.28) is equivalent to the system of integral equations:

$$ S_{\gamma^{\epsilon}} p^{\epsilon}(x) + \sum_{m=1}^{N} (S_{\gamma_m^{\epsilon}} p_m^{\epsilon})_{|\gamma^{\epsilon}}(x) + (S_{\Gamma} q^{\epsilon})_{|\gamma^{\epsilon}}(x) = x_1 + b^{\epsilon}, \quad x \in \gamma^{\epsilon}, $$

$$ (S_{\gamma^{\epsilon}} p^{\epsilon})_{|\gamma_{\ell}^{\epsilon}}(x) + S_{\gamma_{\ell}^{\epsilon}} p_{\ell}^{\epsilon}(x) + \sum_{\substack{m=1 \\ m \neq \ell}}^{N} (S_{\gamma_m^{\epsilon}} p_m^{\epsilon})_{|\gamma_{\ell}^{\epsilon}}(x) + (S_{\Gamma} q^{\epsilon})_{|\gamma_{\ell}^{\epsilon}}(x) = b_{\ell}^{\epsilon}, \quad x \in \gamma_{\ell}^{\epsilon}, $$

$$ (S_{\gamma^{\epsilon}} p^{\epsilon})_{|\Gamma}(x) + \sum_{m=1}^{N} (S_{\gamma_m^{\epsilon}} p_m^{\epsilon})_{|\Gamma}(x) + S_{\Gamma} q^{\epsilon}(x) = 0, \quad x \in \Gamma, $$

with the constraints (see Remark 2) $\int_{\gamma^\epsilon} p^\epsilon(y) \, d\sigma_y = \int_{\gamma_m^\epsilon} p_m^\epsilon(y) \, d\sigma_y = 0$ for all $m = 1, \dots, N$. Rescaling the above equations, we obtain with obvious notation,

$$ (2.29a) \qquad \int_{\gamma} G(x-y) \epsilon p^{\epsilon}(r + \epsilon(y-r)) \, d\sigma_y \\ \qquad + \sum_{m=1}^{N} \int_{\gamma_m} G(r-r_m + \epsilon(x-r-y+r_m)) \epsilon p_m^{\epsilon}(r_m + \epsilon(y-r_m)) \, d\sigma_y \\ \qquad + \int_{\Gamma} G(r+\epsilon(x-r)-y) q^{\epsilon}(y) \, d\sigma_y = \epsilon x_1 + b^{\epsilon}, \quad x \in \gamma, $$

$$ (2.29b) \qquad \int_{\gamma} G(r_{\ell}-r+\epsilon(x-r_{\ell}-y+r)) \epsilon p^{\epsilon}(r+\epsilon(y-r)) \, d\sigma_y + \int_{\gamma_{\ell}} G(x-y) \epsilon p_{\ell}^{\epsilon}(r_{\ell}+\epsilon(y-r_{\ell})) \, d\sigma_y \\ \qquad + \sum_{\substack{m=1 \\ m \neq \ell}}^{N} \int_{\gamma} G(r_{\ell}-r_m + \epsilon(x-r_{\ell}-y+r_m)) \epsilon p_m^{\epsilon}(r_m + \epsilon(y-r_m)) \, d\sigma_y \\ \qquad + \int_{\Gamma} G(r_{\ell}+\epsilon(x-r_{\ell})-y) q^{\epsilon}(y) \, d\sigma_y = b_{\ell}^{\epsilon}, \quad x \in \gamma_{\ell}, $$
$$ (2.29c) \quad \int_{\gamma} G(x-r-\epsilon(y-r))\epsilon p^{\epsilon}(r+\epsilon(y-r)) d\sigma_y \\ + \sum_{m=1}^{N} \int_{\gamma_m} G(x-r_m-\epsilon(y-r_m))\epsilon p^{\epsilon}(r_m+\epsilon(y-r_m)) d\sigma_y \\ + \int_{\Gamma} G(x-y)q^{\epsilon}(y) d\sigma_y = 0, \quad x \in \Gamma, $$

with the constraints

$$ (2.29d) \quad \int_{\gamma} p_{\epsilon}(r + \epsilon(y - r)) d\sigma_y = \int_{\gamma_m} p_{\epsilon}^{m}(r + \epsilon(y - r)) d\sigma_y = 0, \quad m = 1, \dots, N. $$

For all $\ell = 1, \dots, N$, we introduce the formal asymptotic expansions of the quantities:

$$ (2.30a) \quad \epsilon p^{\epsilon}(r + \epsilon(y - r)) = p_0(y) + \epsilon p_1(y) + \epsilon^2 \mathcal{P}^{\epsilon}(y), \quad y \in \gamma, $$

$$ (2.30b) \quad \epsilon p_{\ell}^{\epsilon}(r + \epsilon(y - r)) = p_{\ell,0}(y) + \epsilon p_{\ell,1}(y) + \epsilon^2 \mathcal{P}_{\ell}^{\epsilon}(y), \quad y \in \gamma_{\ell}, $$

$$ (2.30c) \quad q^{\epsilon}(y) = q_0(y) + \epsilon q_1(y) + \epsilon^2 q_2(y) + \epsilon^3 \mathcal{Q}^{\epsilon}(y), \quad y \in \Gamma, $$

$$ (2.30d) \quad b^{\epsilon} = b_0 + \epsilon b_1 + \epsilon^2 \mathcal{B}^{\epsilon}, $$

$$ (2.30e) \quad b_{\ell}^{\epsilon} = b_{\ell,0} + \epsilon b_{\ell,1} + \epsilon^2 \mathcal{B}_{\ell}^{\epsilon}. $$

Using these expressions in system (2.29), the order 0 (in $\epsilon$) reads

$$ \begin{aligned} S_{\gamma} p_0(x) + S_{\Gamma} q_0(r) &= b_0, & x \in \gamma, \\ S_{\gamma_{\ell}} p_{\ell,0}(x) + S_{\Gamma} q_0(r_{\ell}) &= b_{\ell,0}, & x \in \gamma_{\ell}, \quad \ell = 1, \dots, N, \\ S_{\Gamma} q_0(x) &= 0, & x \in \Gamma, \end{aligned} $$

which, thanks to (2.29d), leads to $q_0 = p_0 = p_{\ell,0} = 0$ and $b_0 = b_{\ell,0} = 0$. Identifying the terms of order 1 in $\epsilon$ of system (2.29) leads to

$$ \begin{aligned} S_{\gamma} p_1(x) + S_{\Gamma} q_1(r) &= x_1 + b_1, & x \in \gamma, \\ S_{\gamma_{\ell}} p_{\ell,1}(x) + S_{\Gamma} q_1(r_{\ell}) &= b_{\ell,1}, & x \in \gamma_{\ell}, \quad \ell = 1, \dots, N, \\ S_{\Gamma} q_1(x) &= 0, & x \in \Gamma. \end{aligned} $$

Multiplying the first equation by the $\psi_{\gamma}^{\text{eq}}$ of $\gamma$ and integrating over $\gamma$, we find that $b_1 = -\int_{\gamma} \psi_{\gamma}^{\text{eq}}(x) x_1 d\sigma_x = -r \cdot e_1$, since the first momentum of the equilibrium density of a circle is its center. From the other equations we deduce that $q_1 = 0$, $p_{\ell,1} = 0$, and $b_{\ell,1} = 0$. Finally, $S_{\gamma} p_1(x) = (x-r) \cdot e_1$ for $x \in \gamma$, and according to Lemma B.1, we have $p_1(x) = (2/R)(x-r) \cdot e_1$. Finally, the identification of the second order terms in $\epsilon$ in (2.29c) yields

$$ \nabla G(x - r) \cdot \left( \int_{\gamma} (r - y) p_1(y) d\sigma_y \right) + S_{\Gamma} q_2(x) = 0, \quad x \in \Gamma. $$

Identity (B.4) of Lemma B.1 shows that $\int_{\gamma} (r-y)p_1(y) d\sigma_y = -2\pi R^2 e_1$, and therefore, $q_2$ is the unique solution of the integral equation

$$ S_{\Gamma} q_2(x) = 2\pi R^2 \nabla G(x - r) \cdot e_1, \quad x \in \Gamma. $$
The remainders in (2.30) can be shown to be bounded using the same arguments as in section 2.1. Summing up, we have proved the following result.

PROPOSITION 2.7. With the above notation, we have for every $\ell = 1, \dots, N$:

$$ (2.31\text{a}) \qquad \varepsilon p^{\varepsilon}(r + \varepsilon(y - r)) = \varepsilon p_1(y) + \varepsilon^2 \mathcal{P}^{\varepsilon}(y), \qquad y \in \gamma, $$

$$ (2.31\text{b}) \qquad \varepsilon p_{\ell}^{\varepsilon}(r + \varepsilon(y - r)) = \varepsilon^2 \mathcal{P}_{\ell}^{\varepsilon}(y), \qquad y \in \gamma_{\ell}, $$

$$ (2.31\text{c}) \qquad q^{\varepsilon}(y) = \varepsilon^2 q_2(y) + \varepsilon^3 \mathcal{Q}^{\varepsilon}(y), \qquad y \in \Gamma, $$

$$ (2.31\text{d}) \qquad b^{\varepsilon} = \varepsilon b_1 + \varepsilon^2 \mathcal{B}^{\varepsilon}, $$

$$ (2.31\text{e}) \qquad b_{\ell}^{\varepsilon} = \varepsilon^2 \mathcal{B}_{\ell}^{\varepsilon}, $$

where

$$ b_1 = -r \cdot e_1, \qquad p_1(x) = \frac{2(x-r) \cdot e_1}{R}, \qquad q_2 = S_{\Gamma}^{-1} (2\pi R^2 \nabla G(\cdot - r) \cdot e_1), $$

and where the quantities $\|\mathcal{P}^{\varepsilon}\|_{H^{-1/2}(\gamma)}$, $\|\mathcal{P}_{\ell}^{\varepsilon}\|_{H^{-1/2}(\gamma_{\ell})}$, $\|\mathcal{Q}^{\varepsilon}\|_{H^{-1/2}(\Gamma)}$, $\mathcal{B}^{\varepsilon}$, and $\mathcal{B}_{\ell}^{\varepsilon}$ are uniformly bounded in $\varepsilon$.

Using these expansions and similar formulae for $w_m^{\varepsilon}$, one can easily check that the contributions $[v_m^{\varepsilon}, U^g]_{\gamma_{\ell}^{\varepsilon}}$ and $[w_m^{\varepsilon}, U^g]_{\gamma_{\ell}^{\varepsilon}}$ to the expansion of $\Lambda^{\varepsilon}$ (see (2.2) and (2.6)) read as follows.

THEOREM 2.8. Let $\ell \in \{1, \dots, M\}$ and let $\mathcal{I}_{\ell}$ be defined in (2.3). For all $m \in \{1, \dots, M\}$, we have

$$ V_m^{\perp} \cdot \begin{pmatrix} [v_m^{\varepsilon}, U^g]_{\gamma_{\ell}^{\varepsilon}} \\ [w_m^{\varepsilon}, U^g]_{\gamma_{\ell}^{\varepsilon}} \end{pmatrix} = \begin{cases} O(\varepsilon^3) & \text{if } m \neq \ell, \\ \varepsilon^2 (2\pi R_{\ell}^2 \nabla U^g(r_{\ell}) \cdot V_{\ell}^{\perp}) + O(\varepsilon^3) & \text{if } m = \ell. \end{cases} $$

According to (2.2) and to the splitting (2.6), Theorem 2.1 follows then immediately from Theorems 2.6 and 2.8.

### 3. Target identification using exterior sources.
Our goal is to determine the number $M$, the positions $r_m$, the rescaled radii $R_m$, and the velocities $V_m$ of the small moving disks from the knowledge of the bilinear form defined for all $f, g \in H^{1/2}(\Gamma)$ by

$$ (3.1) \qquad a(f,g) := \sum_{m=1}^{M} R_m^2 \nabla U^f(r_m) \cdot \nabla U^g(r_m) + \sum_{m=1}^{M} R_m^2 \nabla U^g(r_m) \cdot V_m^{\perp}. $$

Let us start by determining the positions which are the hardest parameters to reconstruct. We first note that we can isolate the first sum in the right-hand side of (3.1),

$$ (3.2) \qquad a_0(f,g) := a(f,g) - a(0,g) = \sum_{m=1}^{M} R_m^2 \nabla u^f(r_m) \cdot \nabla u^g(r_m), $$

and hence get a measurement in which the unknown velocities $V_m$ do not appear anymore. Actually, $a_0(f,g)$ coincides to the case of motionless targets. Of course, this does not mean that we have to realize two experiments (one with moving targets and the other one with still ones), but simply that the available data (the full DtN map) can be used to obtain the new data $a_0(\cdot, \cdot)$. The next step is to relate these Laplace type data to a Helmholtz type inverse problem. To do so, following an idea
due to Calderón [8], we use a suitably chosen family of exponential type excitations $f$ and test functions $g$. Next, we make use of the so-called DORT method to recover the unknown positions $r_m$. Based on time-reversal techniques, this method has been successively used for the detection of well separated point like scatterers from far field measurements, for many wave type systems in acoustics [21, 26, 7] and electromagnetics [4]. The underlying idea of this method is that the eigenfunctions of some finite rank integral operator (the so-called time-reversal operator, that can be computed from the measurements) generate waves that selectively focus on each target. For the elliptic Laplace problem considered here, let us emphasize that time reversal is not performed experimentally, but only used numerically in the reconstruction algorithm through the use of suitably chosen excitations and test functions. Of course, one can also use other reconstruction methods which are classically used for the detection of point-like targets in scattering theory, like the MUSIC algorithm [13, 23, 20].

**3.1. Motionless targets.** Here, we exclusively use the quantity $a_0(f,g)$ defined by (3.2), which corresponds to the zero velocity case. Following an idea introduced by Calderón [8], let us choose for every given $\eta \in \mathbb{R}^2$ the excitation data $f$ and the test functions $g$ defined by

$$ (3.3) \qquad f(x) = e^{i(\eta+i\eta^{\perp})\cdot x}, \quad g(x) = e^{i(\eta-i\eta^{\perp})\cdot x}, \quad x \in \Gamma. $$

Noting that these expressions define, in fact, harmonic functions on $\mathbb{R}^2$, we have $U^f(x) = f(x)$ and $U^g(x) = g(x)$ and hence

$$ a_0(f,g) = -\sum_{m=1}^{M} 2|\eta|^2 R_m^2 e^{2i\eta \cdot r_m}. $$

Choosing $k > 0$ and $\eta = \frac{1}{2}k(\alpha - \beta)$, where $\alpha, \beta$ are two given vectors of the unit sphere $\mathbf{S}$ of $\mathbb{R}^2$ (so that $|\eta|^2 = \frac{1}{2}k^2(1 - \alpha \cdot \beta)$), we can compute from the measurements

$$ A(\alpha, \beta) := -\frac{1}{k^2} a_0(f,g) = \sum_{m=1}^{M} R_m^2 (1 - \alpha \cdot \beta) e^{ik(\alpha-\beta)\cdot r_m}. $$

Define the self-adjoint finite rank integral operator $\mathcal{A} \in \mathcal{L}(\mathbf{S})$ associated with the degenerate kernel $A(\alpha, \beta)$,

$$ (\mathcal{A}\varphi)(\beta) := \int_{\mathbf{S}} A(\alpha, \beta)\varphi(\alpha) \, d\alpha. $$

In order to determine its spectrum and its eigenfunctions, we introduce the following functions of $L^2(\mathbf{S})$ for every $m = 1, \dots, M$:

$$ \varphi_{m,0}(\alpha) = \frac{1}{\sqrt{2\pi}} e^{-ik\alpha \cdot r_m}, \quad \varphi_{m,1}(\alpha) = \frac{1}{\sqrt{\pi}} \alpha_1 e^{-ik\alpha \cdot r_m}, \quad \varphi_{m,2}(\alpha) = \frac{1}{\sqrt{\pi}} \alpha_2 e^{-ik\alpha \cdot r_m}. $$

Then, setting $\lambda_{m,0} = 2\pi R_m^2$ and $\lambda_{m,l} = -\pi R_m^2$ for $l = 1, 2$, we can write $A(\alpha, \beta)$ in the more compact form $A(\alpha, \beta) = \sum_{m=1}^{M} \sum_{l=1}^{3} \lambda_{m,l} \varphi_{m,l}(\alpha)\varphi_{m,l}(\beta)$, which yields

$$ \mathcal{A}\varphi = \sum_{m=1}^{M} \sum_{l=1}^{3} \lambda_{m,l} (\varphi, \varphi_{m,l})_{L^2(\mathbf{S})} \varphi_{m,l}. $$
A being at most of rank 3M, it has at most 3M nonzero eigenvalues. More precisely, for every $m = 1, \dots, M$, $\varphi_{m,0}$, $\varphi_{m,1}$, and $\varphi_{m,2}$ are orthonormal in $L^2(S)$ and

$$ (\varphi_{m,l}, \varphi_{n,l'})_{L^2(S)} = \int_S \psi_{l,l'}(\alpha) e^{ik\alpha \cdot (r_n - r_m)} d\alpha, \quad 1 \le m \ne n \le M, $$

where $\psi_{l,l'}(\alpha)$, for $l, l' \in \{0, 1, 2\}$, are given in Table 1 (with $\alpha = (\alpha_1, \alpha_2)$).

<table><caption>TABLE 1
The values of $\psi_{l,l'}(\alpha)$.</caption><thead><tr><th>(l, l')</th><th>l' = 0</th><th>l' = 1</th><th>l' = 2</th></tr></thead><tbody><tr><th>l = 0</th><td>1</td><td>α<sub>1</sub></td><td>α<sub>2</sub></td></tr><tr><th>l = 1</th><td>α<sub>1</sub></td><td>(α<sub>1</sub>)<sup>2</sup></td><td>α<sub>1</sub>α<sub>2</sub></td></tr><tr><th>l = 2</th><td>α<sub>2</sub></td><td>α<sub>1</sub>α<sub>2</sub></td><td>(α<sub>2</sub>)<sup>2</sup></td></tr></tbody></table>

One can easily check that (see, for instance, [1, p. 360])

$$ (\varphi_{m,0}, \varphi_{n,1})_{L^2(S)} = 2 \int_0^\pi \cos \theta e^{ik|r_m - r_n|} \cos \theta \, d\theta = -2i\pi J_1(k|r_m - r_n|), $$

where $J_n$ denotes the Bessel function of the first kind of order $n$, while

$$ (\varphi_{m,0}, \varphi_{n,2})_{L^2(S)} = (\varphi_{m,1}, \varphi_{n,2})_{L^2(S)} = 0. $$

This shows that for $1 \le m \ne n \le M$, we have (see [1, p. 364])

$$ (\varphi_{m,l}, \varphi_{n,l'})_{L^2(S)} = \mathcal{O}((k|r_m - r_n|)^{-1/2}), \quad k \to +\infty. $$

Summing up, this proves that

$$ \mathcal{A}\varphi_{m,l} = \lambda_{m,l} \varphi_{m,l} + \mathcal{O}((kr_{\min})^{-1/2}), \quad k \to +\infty, $$

where $r_{\min} := \min_{m \ne n} |r_m - r_n|$ is the minimal distance between the solids. In other words, for every $m = 1, \dots, M$, $\varphi_{m,0}$, $\varphi_{m,1}$, and $\varphi_{m,2}$ constitute approximate eigenfunctions of A as $k \to +\infty$ (for high frequency and distant disks), the corresponding eigenvalues being, respectively, $2\pi R_m^2$ (simple eigenvalue) and $-\pi R_m^2$ (with multiplicity two). These approximate eigenfunctions can be used to recover the positions of the disks by constructing the corresponding Herglotz wave

$$ (3.4) \qquad u_{m,l}(x) := \int_S \varphi_{m,l}(\alpha) e^{ik\alpha \cdot x} d\alpha, \quad x \in \mathbb{R}^2, $$

we see that

$$ u_{m,l}(x) = \mathcal{O}((k|x - r_m|)^{-1/2}), \quad k \to +\infty, \quad x \neq r_m. $$

Hence, $u_{m,l}$ generates a wave that selectively focuses on the target $m$ as $k \to \infty$ and allows us to recover its location. In particular, for $l = 0$, we have

$$ u_{m,0}(x) = \frac{1}{\sqrt{2\pi}} \int_S e^{ik\alpha \cdot (x - r_m)} d\alpha = 2\pi J_0(k|x - r_m|), $$

and thus, $u_{m,0}$ reaches its maximum exactly at the point $x = r_m$.
**3.2. Moving targets.** Now that we have determined the positions $r_m$ ($m = 1, \dots, M$) of the solids, its remains to compute their rescaled radii $R_m$ and their velocities $V_m$. To this purpose, we introduce the complex polynomials:

$$P_m(z) = \prod_{1 \le l \ne m \le M} (z - r_l)^2$$

and the following harmonic functions (we identify $x$ and $r_j$ with their complex representations):

$$g_m^1(x) := \Re(P_m(x)/P_m'(r_m)), \quad g_m^2(x) := \Im(P_m(x)/P_m'(r_m)).$$

These functions, which are well defined since $P_m'(r_m) \ne 0$, enjoy the properties

$$U g_m^j = g_m^j, \quad \nabla U g_m^1(r_l) = \begin{pmatrix} \delta_{m,l} \\ 0 \end{pmatrix}, \quad \nabla U g_m^2(r_l) = \begin{pmatrix} 0 \\ \delta_{m,l} \end{pmatrix},$$

for every $m, l = 1, \dots, M$ and $j = 1, 2$ ($\delta_{m,l}$ denotes Kronecker's symbol). Substituting these last two relations in (3.1) and (3.2), the expressions of the radii and the velocities follow then easily, as for every $m = 1, \dots, M$:

$$(3.5) \quad R_m = \sqrt{a_0(g_m^1, g_m^1)}, \quad V_m \cdot e_1 = \frac{a(0, g_m^2)}{R_m^2}, \quad V_m \cdot e_2 = -\frac{a(0, g_m^1)}{R_m^2}.$$

**4. Numerical tests.** In this section, we present some numerical results to illustrate our reconstruction method. The solutions of the boundary value problems involved, in particular for the data generation, are computed using a MATLAB boundary integral equation solver.<fnref n="2" /> As a typical configuration, we consider three small disks located in a (smoothened) rectangular domain $\Omega$. The three small disks are of radii $(m+1)\epsilon$, where $\epsilon = 10^{-3}$ and $m = 1, 2, 3$. They are located, respectively, at the points $(0.7, -0.1)$, $(-0.7, 0.3)$, and $(0, -0.3)$ and their velocities are $(0, 1)$, $(1, 1)$, and $(1, 0)$. As a preliminary step, we first generate numerically the data that will be used for the reconstruction. Given a wavenumber $k$ and a uniform discretization of $[0, 2\pi]$ with mesh size $h = 1/(N+1)$, we compute the matrix $A = (A_{i,j})_{1 \le i,j \le N}$ corresponding to the discretization of the kernel $A(\alpha, \beta)$. More precisely, we have $A_{i,j} = -1/k^2(\Lambda^\epsilon f_{i,j}, g_{i,j})$, where $f_{i,j}, g_{i,j}$ are the functions obtained in (3.3) for $\eta = k(\alpha_i - \beta_j)/2$, where $\alpha_i$ and $\beta_j$ belong to the chosen discretization of $[0, 2\pi]$. This requires solving the forward problem (1.4) for $N^2$ different right-hand sides $f_{i,j}$, $1 \le i, j \le N$. The data are generated using $N = 40$ angles and then perturbed them artificially by adding 10% of noise.

Our DORT based reconstruction procedure is applied to these noisy data. We compute the eigenvalues and eigenfunctions of the discretized integral operator $A_h$ associated to the kernel $\mathcal{A}$. As expected from the theoretical analysis of section 3, it turns out that this integral operator has indeed three significant eigenvalues for a wide range of values of $k$.

We show on Figure 2 the dependence of the positive eigenvalues with respect to the wavelength $\lambda = 2\pi/k$ of the oscillating source (without noise). According to this figure, we can recover the number of disks as soon as $k$ is chosen in the shaded region on Figure 2, roughly corresponding to a constant number of significant eigenvalues. The

<footnote n="2">²For more information, see http://iecl.univ-lorraine.fr/~Alexandre.Munnier/IES/.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 2. Dependence of the positive eigenvalues $\lambda_{m,0}$ (in log scale) on the source's wavelength $\lambda = 2\pi/k$. The shaded region corresponds to the region of admissible wavenumbers leading to identifiability of the number of targets.</figcaption></figure>

wavenumber is now chosen to be $k = 10$, which yields accurate location identification. The centers of the disks can be obtained as the points where the Herglotz waves (3.4) generated by the eigenfunctions associated to the largest eigenvalue reach their maximum (in modulus), as shown on Figure 3. We obtain the following estimated positions $(0.0084, -0.2955)$, $(-0.7075, 0.3099)$, and $(0.7123, -0.1017)$, estimated radii $1.9299, 2.9140, 4.0851$ and estimated velocities $(-0.0143, 1.0322)$, $(1.0129, 1.0143)$, and $(0.9569, -0.0594)$. Finally, let us emphasize that although our reconstruction method is theoretically justified only for small disks $(\epsilon \to 0)$, it turns out to still be efficient numerically for “extended” disks as long as the centers are concerned.

## Appendix A. Proof of Lemma 2.3. Define the maps

$$ \mathcal{K}: p \in H^{-1/2}(\Gamma_1) \mapsto \int_{\Gamma_1} K(\cdot, y)p(y) \, d\sigma_y \in H^1(\mathcal{V}_2), $$

$$ \mathcal{L}: u \in H^1(\mathcal{V}_2) \mapsto \int_{\mathcal{V}_2} \nabla_x K(x, \cdot) \cdot \nabla u(x) + K(x, \cdot) u(x) \, dx \in H^1(\mathcal{V}_1), $$

and $\mathcal{L}^\dagger := \gamma_{\Gamma_1} \circ \mathcal{L}$, where $\gamma_{\Gamma_1}$ is the trace operator from $H^1(\mathcal{V}_1)$ into $H^{1/2}(\Gamma_1)$. Let us verify that $\mathcal{L}$ is continuous. For smooth functions $u$, we have

$$ \|\mathcal{L}u\|_{H^1(\mathcal{V}_1)}^2 = \int_{\mathcal{V}_1} \left( \int_{\mathcal{V}_2} \nabla_x K(x, y) \cdot \nabla u(x) + K(x, y)u(x) \, dx \right)^2 dy \\ + \int_{\mathcal{V}_1} \left| \int_{\mathcal{V}_2} D_{xy}^2 K(x, y) \nabla u(x) + \nabla_y K(x, y) u(x) \, dx \right|^2 dy, $$

where $D_{xy}^2 K := (\partial_{x_i y_j} K)_{1 \le i, j \le 2}$. Applying Jensen then the Cauchy–Schwarz inequality, we get

$$ \|\mathcal{L}u\|_{H^1(\mathcal{V}_1)} \le C \|K\|_{H^2(\mathcal{V}_1 \times \mathcal{V}_2)} \|u\|_{H^1(\mathcal{V}_2)}, $$

where $C$ depends only on $\mathcal{V}_1$ and $\mathcal{V}_2$. By density, this estimate remains true for every $u \in H^1(\mathcal{V}_2)$. For smooth $p$ and $u$, one can easily verify that

$$ (Kp, u)_{H^1(\mathcal{V}_2)} = \langle p, \mathcal{L}^\dagger u \rangle_{H^{-1/2}(\Gamma_1) \times H^{1/2}(\Gamma_1)}. $$
<figure><img src="image_4.png" /><figcaption>FIG. 3. Three dimensional (left) and two dimensional (right) visualization of the focusing properties of the eigenfunctions corresponding to the three largest positive eigenvalues for $k = 10.5$. The focusing property leads to an accurate recovery of the centers of the disks.</figcaption></figure>

We deduce from the continuity of the trace operator $\gamma_{\Gamma_1}: H^1(\mathcal{V}_1) \to H^{1/2}(\Gamma_1)$, that

$$|(Kp, u)_{H^1(\mathcal{V}_2)}| \le C \|K\|_{H^2(\mathcal{V}_1 \times \mathcal{V}_2)} \|p\|_{H^{-1/2}(\Gamma_1)} \|u\|_{H^1(\mathcal{V}_2)},$$

where $C > 0$ does not depend on $p, u$, and $K$. This inequality ensures the continuity of $K$. We get the conclusion of the lemma by writing that

$$\int_{\Gamma_2} \int_{\Gamma_1} K(x, y)p(y)q(x) \, d\sigma_y d\sigma_x = \int_{\Gamma_2} (\gamma_{\Gamma_2} \circ Kp)(x)q(x) \, d\sigma_x,$$

where $\gamma_{\Gamma_2}: H^1(\mathcal{V}_2) \to H^{1/2}(\Gamma_2)$ is the continuous trace operator.

## Appendix B. Single layer integral equation for a circle.

LEMMA B.1. Denote by $\gamma$ the circle of radius $R$ centered at the origin and let $c \in \mathbb{R}$ and $V \in \mathbb{R}^2$. Then, the solution of the integral equation

$$(B.1) \qquad S_{\gamma}p(x) = c + V \cdot x \quad \forall x \in \gamma$$
is given by

$$ (B.2) \qquad p(x) = \frac{c}{R} + \frac{2}{R} V \cdot x \quad \forall x \in \gamma. $$

The corresponding simple layer potential is

$$ (B.3) \qquad (S_{\gamma}p)(x) = \begin{cases} c + V \cdot x & \text{for } |x| < R, \\ c \log\left(\frac{eR}{|x|}\right) + \left(\frac{R}{|x|}\right)^2 V \cdot x & \text{for } |x| > R. \end{cases} $$

Furthermore, we have

$$ (B.4) \qquad \int_{\gamma} xp(x) d\sigma_x = 2\pi R^2 V. $$

*Proof.* One can easily check that the function defined by (B.3) satisfies

$$ \begin{align*} -\Delta\varphi &= 0 & \text{in } \mathbb{R}^2 \setminus \gamma, \\ \varphi &= c + V \cdot x & \text{on } \gamma, \\ [\varphi] &= 0 & \text{on } \gamma, \\ \varphi(x) &= O(\log|x|) & \text{as } |x| \to +\infty. \end{align*} $$

Hence $\varphi = S_{\gamma}p$, where $p$ solves (B.1). Equation (B.2) follows then from $p = \left[\frac{\partial\varphi}{\partial n}\right]$ (see Proposition 2.2). Using (B.2), (B.4) follows from a straightforward computation. $\square$

**Acknowledgment.** The authors thank the referees for their careful reading and helpful comments and suggestions.

## REFERENCES

[1] M. ABRAMOWITZ AND I. A. STEGUN, *EDS.*, *Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables*, National Bureau of Standards Applied Mathematics Series 55, U.S. Government Printing Office, Washington, DC, 1964.

[2] H. AMMARI AND H. KANG, *Generalized polarization tensors, inverse conductivity problems, and dilute composite materials: A review*, in Inverse Problems, Multi-scale Analysis and Effective Medium Theory, Contemp. Math. 408, AMS, Providence, RI, 2006, pp. 1–67.

[3] H. AMMARI AND H. KANG, *Polarization and Moment Tensors*, Appl. Math. Sci. 162, Springer, New York, 2007.

[4] X. ANTOINE, B. PINÇON, K. RAMDANI, and B. THIERRY, *Far field modeling of electromagnetic time reversal and application to selective focusing on small scatterers*, SIAM J. Appl. Math., 69 (2008), pp. 830–844.

[5] M. BADRA, F. CAUBET, and M. DAMBRINE, *Detecting an obstacle immersed in a fluid by shape optimization methods*, Math. Models Methods Appl. Sci., 21 (2011), pp. 2069–2101.

[6] L. BOURGEOIS and J. DARDÉ, *The “exterior approach” to solve the inverse obstacle problem for the Stokes system*, Inverse Probl. Imaging, 8 (2014), pp. 23–51.

[7] C. BURKARD, A. MINUT, and K. RAMDANI, *Far field model for time reversal and application to selective focusing on small dielectric inhomogeneities*, Inverse Probl. Imaging, 7 (2013), pp. 445–470.

[8] A.-P. CALDERÓN, *On an inverse boundary value problem*, in Seminar on Numerical Analysis and Its Applications to Continuum Physics (Rio de Janeiro, 1980), Soc. Brasil. Mat., Rio de Janeiro, 1980, pp. 65–73.

[9] F. CAUBET and M. DAMBRINE, *Localization of small obstacles in Stokes flow*, Inverse Problems, 28 (2012), 105007.
[10] F. CAUBET AND M. DAMBRINE, *Stability of critical shapes for the drag minimization problem in Stokes flow*, J. Math. Pures Appl. (9), 100 (2013), pp. 327-346.

[11] F. CAUBET, M. DAMBRINE, AND D. KATEB, *Shape optimization methods for the inverse obstacle problem with generalized impedance boundary conditions*, Inverse Problems, 29 (2013), 115011.

[12] F. CAUBET, M. DAMBRINE, D. KATEB, AND C. Z. TIMIMOUN, *A Kohn-Vogelius formulation to detect an obstacle immersed in a fluid*, Inverse Probl. Imaging, 7 (2013), pp. 123-157.

[13] M. CHENEY, *The linear sampling method and the MUSIC algorithm*, Inverse Problems, 17 (2001), pp. 591-595.

[14] C. CONCA, P. CUMSILLE, J. ORTEGA, AND L. ROSIER, *Corrigendum: On the detection of a moving obstacle in an ideal fluid by a boundary measurement*, Inverse Problems, 24 (2008), 059802.

[15] C. CONCA, P. CUMSILLE, J. ORTEGA, AND L. ROSIER, *On the detection of a moving obstacle in an ideal fluid by a boundary measurement*, Inverse Problems, 24 (2008), 045001.

[16] C. CONCA, M. MALIK, AND A. MUNNIER, *Detection of a moving rigid solid in a perfect fluid*, Inverse Problems, 26 (2010), 095010.

[17] C. CONCA, E. L. SCHWINDT, AND T. TAKAHASHI, *On the identifiability of a rigid body moving in a stationary viscous fluid*, Inverse Problems, 28 (2012), 015005.

[18] M. FINK AND C. PRADA, *Eigenmodes of the time-reversal operator: A solution to selective focusing in multiple-target media*, Wave Motion, 20 (1994), pp. 151-163.

[19] A. FRIEDMAN AND M. VOGELIUS, *Identification of small inhomogeneities of extreme conductivity by boundary measurements: A theorem on continuous dependence*, Arch. Rational Mech. Anal., 105 (1989), pp. 299-326.

[20] R. GRIESMAIER AND M. HANKE, *MUSIC characterization of small scatterers for normal measurement data*, Inverse Problems, 25 (2009), 075012.

[21] C. HAZARD AND K. RAMDANI, *Selective acoustic focusing using time-harmonic reversal mirrors*, SIAM J. Appl. Math., 64 (2004), pp. 1057-1076.

[22] A. M. IL'IN, *Investigation of the asymptotic behavior of the solution of an elliptic boundary-value problem in a domain with a small hole*, J. Soviet Math., 33 (1986), pp. 994-1014.

[23] A. KIRSCH, *The MUSIC algorithm and the factorization method in inverse scattering theory for inhomogeneous media*, Inverse Problems, 18 (2002), pp. 1025-1040.

[24] V. G. MAZ'YA, S. A. NAZAROV, AND B. A. PLAMENEVSKII, *Asymptotic Theory of Elliptic Boundary Value Problems in Singularly Perturbed Domains*. Vol. I, Oper. Theory Adv. Appl. 111, Birkhäuser Verlag, Basel, 2000.

[25] W. MCLEAN, *Strongly Elliptic Systems and Boundary Integral Equations*, Cambridge University Press, Cambridge, 2000.

[26] B. PINÇON AND K. RAMDANI, *Selective focusing on small scatterers in acoustic waveguides using time reversal mirrors*, Inverse Problems, 23 (2007), pp. 1-25.

[27] S. RJASANOW AND O. STEINBACH, *The Fast Solution of Boundary Integral Equations*, Math. Anal. Tech. Appl. Eng., Springer, New York, 2007.

[28] O. STEINBACH, *Numerical Approximation Methods for Elliptic Boundary Value Problems*, Springer, New York, 2008.