# Reconstruction of obstacles in a Stokes flow as a shape-from-moments problem

Alexandre Munnier

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

## To cite this version:

Alexandre Munnier. Reconstruction of obstacles in a Stokes flow as a shape-from-moments problem. Inverse Problems and Imaging, 2024, 18 (6), pp.1343-1365. 10.3934/ipi.2024017. hal-04218495

HAL Id: hal-04218495

https://hal.science/hal-04218495v1

Submitted on 26 Sep 2023

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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Reconstruction of obstacles in a Stokes flow as a shape-from-moments problem

Alexandre Munnier*

Université de Lorraine, CNRS, Inria, IECL, F-54000 Nancy, France

September 26, 2023

## Abstract

We address the classical inverse problem of recovering the position and shape of obstacles immersed in a planar Stokes flow using boundary measurements. We prove that this problem can be transformed into a shape-from-moments problem to which ad hoc reconstruction methods can be applied. The effectiveness of this approach is confirmed by numerical tests that show significant improvements over those available in the literature to date.

**Keywords**— Geometric inverse problem, Stokes equations, non-primitive variables, shape-from-moments problem, bi-harmonic single-layer potential, Prony's system, partial balayage.

## 1 Introduction

A shape-from-moments problem is a geometric inverse problem consisting in recovering the shape of an unknown (possibly multi-connected) domain $O$ from a finite section of its complex moments:

$$ \int_{O} \overline{z^{m}z^{n}} dm(z) \quad (n, m \in \mathbb{N}). $$

In some cases, only the harmonic moments (i.e. for which $m = 0$ in the identity above) are available. The 2D inverse gravimetric problem is inherently a shape-from-moments problem and can therefore be handled with tools developed for this type of problem as it is explained in the recent article [7]. More surprisingly, the geometric Calderón inverse problem has also been shown in [16] to be equivalent to a shape-from-moments problem, leading to an original and efficient reconstruction method.

In this work, we are going to show that the problem of detecting (and reconstructing) obstacles immersed in a Stokes flow from boundary measurements can be seen (and dealt with) as a shape-from-moments problem as well. Let us go into more detail and start by making precise the geometric parameters and the problem under consideration:

Let $\Omega$ and $O_j$ ($j = 1, \dots, N$) be open Jordan domains of class $C^{1,1}$ such that $\bar{O}_j \subset \Omega$ and $\bar{O}_j \cap \bar{O}_k = \emptyset$ for every indices $j, k = 1, \dots, N, j \neq k$. The boundaries of $\Omega$ and $O_j$ are denoted by $\Gamma_0$ and $\Gamma_j$ respectively and we also denote:

$$ O = \cup_{j=1}^{N} O_j \quad \text{and} \quad \Gamma = \cup_{j=1}^{N} \Gamma_j. $$

The unit normal vector field $n$ defined on $\Gamma$ and on $\Gamma_0$ is always assumed to be pointing towards the interior of the domain enclosed by the Jordan curves (see Fig. 1). We assume that the domain $\mathcal{F} = \Omega \setminus \cup_{j=1}^{N} \bar{O}_j$ is filled with a fluid whose flow is governed by the Stokes equations. Thus, the velocity and pressure fields $u$ and $p$ verify:

$$ \left\{ \begin{array}{ll} -\nu\Delta u + \nabla p = 0 & \text{in } \mathcal{F}, \\ \text{div } u = 0 & \text{in } \mathcal{F}, \\ u = 0 & \text{on } \Gamma, \\ u = v & \text{on } \Gamma_0, \end{array} \right. \quad (1a) \quad (1b) \quad (1c) \quad (1d) $$

where $\nu > 0$ stands for the kinematic viscosity of the fluid and $v$ is a prescribed velocity on $\Gamma_0$ satisfying the flux condition:

$$ \int_{\Gamma_0} v \cdot n \, ds = 0. \quad (1e) $$

<footnote>*alexandre.munnier@univ-lorraine.fr</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1: $\mathcal{F}$ and $\mathcal{O}_j$ are respectively the fluid and obstacle domains, and $\Omega$ is the union of $\mathcal{F}$ and $\overline{\mathcal{O}}_j$ ($j = 1, \dots, N$). The boundaries of the domains are $\mathcal{C}^{1,1}$ Jordan curves. At every point of the boundaries, the unit normal vector $n$ is directed towards the interior of the domain enclosed by the curve.</figcaption></figure>

The Cauchy stress tensor is defined in $\mathcal{F}$ by:

$$ T(u, p) = \nu[\nabla u + (\nabla u)^t] - p \operatorname{Id}. \qquad (2) $$

The classical reconstruction problem we focus on in this article consists in determining the shape and position of the obstacles $\mathcal{O}_j$ by observing the Cauchy force $T(u, p)n$ along $\Gamma_0$. More precisely, introducing the Sobolev space:

$$ \tilde{H}^{1/2}(\Gamma_0; \mathbb{R}^2) = \{v \in H^{1/2}(\Gamma_0; \mathbb{R}^2) : \int_{\Gamma_0} v \cdot n \, ds = 0\}, $$

and denoting by $H^{-1/2}(\Gamma_0; \mathbb{R}^2)$ its dual space (using $L^2(\Gamma_0; \mathbb{R}^2)$ as pivot space), we define the operator:

$$ \tilde{\Lambda}_{\Gamma} : \tilde{H}^{1/2}(\Gamma_0; \mathbb{R}^2) \to H^{-1/2}(\Gamma_0; \mathbb{R}^2) \\ v \mapsto T(u, p)n, $$

where $(u, p)$ is the solution of System (1). We can now state:

**Problem 1.** Reconstruct the multi-connected domain $\mathcal{O}$ from $\tilde{\Lambda}_{\Gamma}$.

In [1] it is established that $\tilde{\Lambda}_{\Gamma}$ uniquely determines the multi-connected domain $\mathcal{O}$. Relevant references on this problem can be found in [5]. As far as reconstruction methods are concerned, most are based on the minimization of a cost functional using shape sensitivity analysis techniques. However, as explained in [4], the problem is severely ill-posed and generally these methods, although capable of roughly capturing the position and shape of the obstacles, give mixed results. To achieve a satisfactory reconstruction, the obstacles generally have to satisfy certain constraints: they must be close to the outer boundary (where the measurements are made), sometimes sufficiently small, sometimes their number must be known a priori.

The new approach we propose in this paper includes the following steps: first, we will show that Problem 1 can be equivalently reformulated in terms of vorticity and stream function (the so-called “non-primitive variables”), leading to an inverse problem for the biharmonic operator. This problem will then be transformed into a system of integral equations involving biharmonic single-layer potentials. Following an idea initially introduced in [15], we will explain how the formulation can be modified so that only integrals on the obstacles boundaries appear and this will allow us to compute, for any harmonic functions $h_1, h_2$ in $\Omega$, the quantities:

$$ \int_{\mathcal{O}} h_1 h_2 \, dm. $$

By choosing harmonic monomials for $h_1$ and $h_2$, we will be able to deduce the complex moments of the obstacles:

$$ \int_{\mathcal{O}} \overline{z^m z^n} \, dm(z) \quad \text{for all } n, m \in \mathbb{N}, \qquad (3) $$

thus transforming the initial reconstruction problem into a classical shape-from-moments problem. At this point, two algorithms of reconstruction can be used. The first one, described in [11] requires a finite section of all the complex moments. The boundaries of the obstacles are obtained as the zero level line of a function constructed from the orthogonal polynomials on $\mathcal{O}$. The second requires only a finite section of the harmonic moments. The obstacles are first considered as Dirac masses whose locations and weights are obtained by solving a so-called Prony's system (a non-linear system
associated to a Vandermonde matrix). Then, the shapes of the obstacles are deduced by applying a “partial-balayage” operator (derived from the theory of “balayage” of measures). This method was successfully applied to solve inverse gravimetric problems in [7] and will be seen to provide fairly accurate results regardless of the size, number and location of the obstacles.

It should be noted that one of the original features of the proposed algorithms is that they are not iterative, unlike most of those available to date. They require a finite number of measurements and provide a direct approximation of the target domains.

## 2 The problem in non-primitive variables

For any vector $x = (x_1, x_2)$ in $\mathbb{R}^2$, $x^\perp$ is the vector rotated counter-clockwise, i.e. $x^\perp = (-x_2, x_1)$. The velocity field $u$ in System (1) satisfies $\mathrm{div} \, u = 0$ in $\mathcal{F}$ and $\int_{\mathcal{C}} u \cdot n \, ds = 0$ for any Jordan curve $\mathcal{C}$ included in $\mathcal{F}$. This ensures the existence of a stream function $\psi$ such that $\nabla^\perp \psi = u$ in $\mathcal{F}$. From (1a), we deduce that:

$$ -\nu \nabla^\perp \omega + \nabla p = 0 \quad \text{in } \mathcal{F}, \qquad (4) $$

where $\omega = \Delta\psi$ is the vorticity field. This implies (by taking the rotational of the identity above) that $\omega$ is harmonic in $\mathcal{F}$ and that:

$$ \int_{\Gamma_j} \partial_n \omega \, ds = 0 \quad \text{for all } j = 1, \dots, N. $$

From the prescribed velocity field $v$ in (1d), we define a function $f^d$ on $\Gamma_0$ such that $\partial_\tau f^d = v \cdot n$, where $\tau$ is the unit tangent vector field on $\Gamma_0$ oriented such that $\tau^\perp = n$. The existence of $f^d$ is guaranteed by the flux condition (1e). Let also introduce $f^n = -v \cdot \tau$ on $\Gamma_0$. The boundary value problem (1) can be restated in terms of non-primitive variables $(\psi, \omega)$ as follows:

$$ \Delta\omega = 0 \quad \text{in } \mathcal{F}, \qquad (5a) $$

$$ \Delta\psi = \omega \quad \text{in } \mathcal{F}, \qquad (5b) $$

$$ (\psi, \partial_n \psi) = (c_j, 0) \quad \text{on } \Gamma_j \text{ for } j = 1, \dots, N, \qquad (5c) $$

$$ (\psi, \partial_n \psi) = (f^d, f^n) \quad \text{on } \Gamma_0, \qquad (5d) $$

where the $c_j$ ($j = 1, \dots, N$) are real constants such that:

$$ \int_{\Gamma_j} \partial_n \omega \, ds = 0 \quad \text{for all } j = 1, \dots, N. \qquad (5e) $$

Extending the stream function $\psi$ by $c_j$ in $O_j$, it can be assumed to be defined in the whole domain $\Omega$. In a similar way, the vorticity field $\omega$ is assumed to be defined in $\Omega$, extended by 0 in the obstacles. In the sequel, we shall also require the functions $\xi_j$ ($j = 0, 1, \dots, N$) solving:

$$ \Delta\omega_j = 0 \quad \text{in } \mathcal{F} \qquad (6a) $$

$$ \Delta\xi_j = \omega_j \quad \text{in } \mathcal{F} \qquad (6b) $$

$$ (\xi_j, \partial_n \xi_j) = (1, 0) \quad \text{on } \Gamma_j, \qquad (6c) $$

$$ (\xi_j, \partial_n \xi_j) = (0, 0) \quad \text{on } \Gamma_k, \quad (k = 0, 1, \dots, N, k \neq j). \qquad (6d) $$

The proof of the existence and uniqueness of the functions $\xi_j$ in $H^2(\mathcal{F})$ is classical (see [8, Proposition 1.3]). Once extended by suitable constants in the obstacles they can be considered as functions of the space $H_0^2(\Omega)$. For every function $u \in H^2(\Omega)$, we can define the Dirichlet and Neumann traces on any curve $\Gamma_k$ ($k = 0, \dots, N$) denoted respectively by $\gamma_{\Gamma_k}^d u$ and $\gamma_{\Gamma_k}^n u$ (the Neumann trace is defined taking into account the orientation of $n$). The total trace operator is next given by:

$$ \gamma_{\Gamma_k} : H^2(\Omega) \to H^{3/2}(\Gamma_k) \times H^{1/2}(\Gamma_k) \\ u \mapsto (\gamma_{\Gamma_k}^d u, \gamma_{\Gamma_k}^n u). $$

In the sequel we denote by $H(\Gamma_k)$ the space $H^{3/2}(\Gamma_k) \times H^{1/2}(\Gamma_k)$ and by $H'(\Gamma_k) = H^{-3/2}(\Gamma_k) \times H^{-1/2}(\Gamma_k)$ its dual space, using $L^2(\Gamma_k) \times L^2(\Gamma_k)$ as pivot space.

**Proposition 2.1.** For any $f = (f^d, f^n) \in H(\Gamma_0)$, there exists a unique function $\psi$ in $H^2(\Omega)$ constant in every $O_j$ ($j = 1, \dots, N$) and such that its restriction to $\mathcal{F}$ solves System (5). It is obtained as the unique function achieving:

$$ \min \left\{ \|\Delta u\|_{L^2(\Omega)} : \gamma_{\Gamma_0} u = f, \gamma_{\Gamma_j} u = (c_j, 0), c_j \in \mathbb{R}, (j = 1, \dots, N), u \in H^2(\Omega) \right\}. \qquad (7) $$

*Proof.* Let $u_f$ be a function in $H^2(\Omega)$ such that $\gamma_{\Gamma_0} u_f = f$ and $u_f = 0$ in $O$. The space $H^2(\Omega)$ is provided with the scalar product:

$$ (\theta_1, \theta_2)_\Omega = (\Delta\theta_1, \Delta\theta_2)_{L^2(\Omega)} + (\gamma_{\Gamma_0}\theta_1, \gamma_{\Gamma_0}\theta_2)_{H(\Gamma_0)}, \quad \theta_1, \theta_2 \in H_0^2(\Omega), $$
where $(\cdot, \cdot)_{H(\Gamma_0)}$ stands for the scalar product in $H(\Gamma_0)$. Identifying $H_0^2(\mathcal{F})$ with the subspace of $H_0^2(\Omega)$ consisting in the functions that vanish in $\mathcal{O}$, we define the subspace of $H^2(\Omega)$:

$$B(\Omega) = (H_0^2(\mathcal{F}) \oplus \langle \xi_1, \dots, \xi_N \rangle)^\perp.$$

It is classical to verify that the function achieving (7) is the orthogonal projection of $u_f$ on $B(\Omega)$ in $H^2(\Omega)$ and that it is unique. □

From identities (4) and (5a), we deduce that $\omega$ and $p$ are conjugate harmonic functions, or equivalently that the complex function $\nu\omega + ip$ is holomorphic in $\mathcal{F}$. It means that (up to a constant) $p$ can be deduced from $\omega$ (see [14, §7.1]). We denote by $D^2\psi$ the Hessian matrix of $\psi$ and by $(\perp)$ the rotation matrix of angle $\pi/2$ counter-clockwise. The Cauchy stress tensor (2) reads:

$$T(\psi, \omega) = \nu [ (\perp) D^2\psi - D^2\psi(\perp) ] - p \operatorname{Id}.$$

Introducing $s$, the arc length parameterization on $\Gamma_0$, we recall the relations:

$$\partial_s n(s) = -\kappa(s)\tau(s) \quad \text{and} \quad \partial_s \tau(s) = \kappa(s)n(s),$$

where the function $\kappa$ is the curvature of $\Gamma_0$. Forming the scalar product of $T(\psi, \omega)n$ with $n$ and $\tau$ on $\Gamma_0$ we obtain:

$$T(\psi, \omega)n \cdot n = 2\nu D^2 \psi n \cdot \tau - p \quad \text{and} \quad T(\psi, \omega)n \cdot \tau = \nu [-D^2 \psi \tau \cdot \tau + D^2 \psi n \cdot n].$$

Using the identities $\partial_s(\nabla\psi \cdot n) = D^2\psi\tau \cdot n - \kappa\partial_\tau\psi$, $\omega = D^2\psi\tau \cdot \tau + D^2\psi n \cdot n$ and $\partial_s(\nabla\psi \cdot \tau) = D^2\psi\tau \cdot \tau + \kappa\partial_n\psi$ we end up with:

$$T(\psi, \omega)n \cdot n = 2\nu \kappa \partial_{\tau} \psi + 2\nu \partial_{s}(\partial_{n} \psi) - p \quad \text{and} \quad T(\psi, \omega)n \cdot \tau = -2\nu \kappa \partial_{n} \psi + 2\nu \partial_{s}(\partial_{\tau} \psi) - \nu \omega.$$

Because of (4), we have on $\Gamma_0$, $\partial_s p = -\nu\partial_n\omega$, and hence $p$ is an antiderivative of $-\nu\partial_n\omega$ on $\Gamma_0$ so we denote $p = -\nu \int \partial_n\omega$. Summarizing, taking into account (5d), we obtain eventually:

$$T(\psi, \omega)n \cdot n = \nu \left[ 2\kappa \partial_s f^d + 2\partial_s f^n + \int \partial_n \omega \right] \quad \text{and} \quad T(\psi, \omega)n \cdot \tau = \nu \left[ 2\partial_{ss}^2 f^d - 2\kappa f^n - \omega \right].$$

So, for prescribed $f^d$ and $f^n$, it is equivalent to measure the Cauchy force exerted by the fluid on $\Gamma_0$ and to measure the pair $(\omega, \partial_n\omega)$ on this boundary. The function $\omega$ in System (5) is harmonic and square-integrable on $\Omega$, so it admits a Dirichlet and a Neumann trace on $\Gamma_0$ and we can define:

$$\begin{gathered} \Lambda_{\Gamma}: H(\Gamma_0) \to H'(\Gamma_0) \\ (f^d, f^n) \mapsto (\partial_n\omega|_{\Gamma_0}, -\omega|_{\Gamma_0}). \end{gathered} \qquad (8)$$

From now on, we shall focus on the following problem:

**Problem 2.** Reconstruct the multi-connected domain $\mathcal{O}$ from $\Lambda_{\Gamma}$.

## 3 Identifiability for the biharmonic inverse problem

Identifiability results for Problem 1 are given in [1, Theorem 3.1] and in [2, Theorem 2.1] and can probably be adapted to Problem 2. We prefer to provide an original proof directly for the biharmonic inverse problem.

**Theorem 3.1.** Let $\mathcal{O}$ and $\mathcal{O}'$ be two sets of obstacles as described in Section 1. Assume that there exists $f = (f^d, f^n) \in H(\Gamma_0)$ and an open set $\mathcal{U}$ in $\mathbb{R}^2$ such that $\mathcal{U} \cap \Gamma_0 \neq \emptyset$, $\mathcal{U} \cap (\mathcal{O} \cup \mathcal{O}') = \emptyset$, $f^d$ non-constant in $\mathcal{U} \cap \Gamma_0$ and:

$$\langle \Lambda_{\Gamma} f, \gamma_{\Gamma_0} \theta \rangle_{\Gamma_0} = \langle \Lambda_{\Gamma'} f, \gamma_{\Gamma_0} \theta \rangle_{\Gamma_0} \quad \text{for all } \theta \in \mathcal{D}(\mathcal{U}).$$

Then $\mathcal{O} = \mathcal{O}'$.

In this statement, the bracket $\langle \cdot, \cdot \rangle_{\Gamma_0}$ stands for the duality pairing in $H'(\Gamma_0) \times H(\Gamma_0)$ that extends the scalar product of $L^2(\Gamma_0) \times L^2(\Gamma_0)$.

*Proof.* Denote by $\psi$ and $\psi'$ the solutions to System (5) corresponding respectively to the obstacles $\mathcal{O}$ and $\mathcal{O}'$ and the same right hand side $f = (f^d, f^n)$ in identity (5d). As already explained earlier, these functions extended by suitable constants in the obstacles are in $H^2(\Omega)$. The function $\varphi = \psi - \psi'$ is furthermore extended by zero outside $\Omega$ and since $\gamma_{\Gamma_0}\varphi = 0$, this function is in $H^2(\mathbb{R}^2)$ and is biharmonic in $\mathbb{R}^2 \setminus (\Gamma_0 \cup \Gamma \cup \Gamma')$. For every $\theta \in \mathcal{D}(\mathcal{U})$, an integration by parts yields:

$$(\Delta\varphi, \Delta\theta)_{L^2(\mathbb{R}^2)} = \langle \Lambda_{\Gamma} f - \Lambda_{\Gamma'} f, \gamma_{\Gamma_0} \theta \rangle_{\Gamma_0} = 0,$$

which means that $\Delta\varphi$ is harmonic in $\mathcal{U}$. But since $\varphi = 0$ in the open set $\mathcal{U} \setminus \bar{\Omega}$, the unique continuation principle entails that the function $\varphi$ vanishes in the whole set $\mathbb{R}^2 \setminus (\bar{\mathcal{O}} \cup \mathcal{O}')$.
For every $j = 1, \dots, N$, we can decompose the boundary $\Gamma_j$ into two parts:

$$ \Gamma_j = (\Gamma_j \cap \overline{\mathcal{O}'} \cup (\Gamma_j \cap (\mathbb{R}^2 \setminus \overline{\mathcal{O}'}))). $$

On $\Gamma_j \cap \overline{\mathcal{O}'}$, $\psi'$ is piecewise constant by definition and this is also true on $\Gamma_j \cap (\mathbb{R}^2 \setminus \overline{\mathcal{O}'}$) because $\psi' = \psi$ on this set. Since $\psi'$ is continuous, $\psi'$ is constant on every curve $\Gamma_j$. According to the Proposition 2.1, this implies that $\|\Delta\psi'\|_{L^2(\Omega)} \le \|\Delta\psi\|_{L^2(\Omega)}$. The role played by $\psi$ and $\psi'$ being symmetric, it turns out that $\|\Delta\psi'\|_{L^2(\Omega)} = \|\Delta\psi\|_{L^2(\Omega)}$ and then $\psi = \psi'$ by uniqueness of the minimum in Proposition 2.1. If there were an index $j \in \{1, \dots, N\}$ such that the open set $\mathcal{O}'_j \setminus \overline{\mathcal{O}}$ were non-empty, the function $\psi$ would be constant on this set (because $\psi'$ is) and therefore constant on the whole set $\Omega \setminus \overline{\mathcal{O}}$. This would imply that $f^d$ is constant, which is excluded by hypothesis. We deduce that $\mathcal{O} \subset \mathcal{O}'$ and, because $\mathcal{O}$ and $\mathcal{O}'$ play symmetric roles, that $\mathcal{O} = \mathcal{O}'$. $\square$

# 4 Biharmonic single-layer potential

This section is devoted to establishing (and recalling) some results on the biharmonic single-layer potential. So we leave aside for a moment the inverse problem we are dealing with and we consider a general framework in which $\Gamma$ represents a finite, disjoint union of $C^{1,1}$ Jordan curves. Sticking to our convention, the unit normal vector $n$ on $\Gamma$ is directed towards the interior of the domain enclosed by the curves. If $u$ is a function defined on both sides of $\Gamma$ and admitting one-sided Dirichlet and Neumann traces on $\Gamma$, we can define the jump of these traces across the curve:

$$ [\partial_n u]_{\Gamma} = \gamma_{\Gamma}^{n} u^{+} - \gamma_{\Gamma}^{n} u^{-} \quad \text{and} \quad [u]_{\Gamma} = \gamma_{\Gamma}^{d} u^{+} - \gamma_{\Gamma}^{d} u^{-}. $$

These quantities are defined piecewise on each connected component of $\Gamma$ and each component is by definition a Jordan curve. The notation $u^{-}$ represents the restriction of $u$ to the domain bounded by the curve, while $u^{+}$ designates the part of the function $u$ outside the curve.

The fundamental solution of the Bilaplacian is defined by:

$$ G(x) = \frac{1}{8\pi} \left[ |x|^2 \ln \frac{|x|}{\kappa_0} + \kappa_1 \right] \quad \text{for all } x \in \mathbb{R}^2, \qquad (9) $$

where $\kappa_0$ and $\kappa_1$ are real constants that will be fixed later on. Using the usual abuse of notation to identify $G(x-y)$ with a two-variables function $G(x, y)$, the biharmonic single-layer potential is defined for every $q = (q_n, q_d) \in H'(\Gamma)$ by:

$$ \mathcal{S}_{\Gamma} q(x) = \int_{\Gamma} G(x, y) q_n(y) + \partial_{n(y)} G(x, y) q_d(y) \, ds(y) \quad \text{for all } x \in \mathbb{R}^2. \qquad (10) $$

The operator $\mathcal{S}_{\Gamma}: H'(\Gamma) \to H_{\text{loc}}^2(\mathbb{R}^2)$ is bounded so the same conclusion applies to the operator:

$$ S_{\Gamma}: H'(\Gamma) \to H(\Gamma) \\ q \mapsto \gamma_{\Gamma} \circ \mathcal{S}_{\Gamma} q. \qquad (11) $$

The results presented without proof in this section are borrowed from [17]:

**Theorem 4.1.** Let $R > 0$ be the radius of a circle $C_R$ that enclosed $\Gamma$. If we choose $\kappa_0 > eR$ and $\kappa_1 > R^2$ in the definition (9), then the operator $S_{\Gamma}$ is strongly elliptic on $H'(\Gamma)$ and therefore invertible.

The so-called jump relations allow recovering $q = (q_n, q_d)$ from the function $\mathcal{S}_{\Gamma}q$:

$$ q_n = -[\partial_n \mathcal{S}_{\Gamma} q]_{\Gamma} = -(\partial_n \mathcal{S}_{\Gamma}^{+} q - \partial_n \mathcal{S}_{\Gamma}^{-} q) \qquad (12a) $$

$$ q_d = [\Delta \mathcal{S}_{\Gamma} q]_{\Gamma} = \Delta \mathcal{S}_{\Gamma}^{+} q - \Delta \mathcal{S}_{\Gamma}^{-} q. \qquad (12b) $$

Under the hypotheses of Theorem 4.1, the spaces $H'(\Gamma)$ and $H(\Gamma)$ can be provided with the scalar products:

$$ (q, q')_{H'(\Gamma)} = \langle q, S_{\Gamma} q' \rangle_{\Gamma}, \quad \text{for all } q, q' \in H'(\Gamma), $$

$$ (p, p')_{H(\Gamma)} = \langle S_{\Gamma}^{-1} p, p' \rangle_{\Gamma}, \quad \text{for all } p, p' \in H(\Gamma), $$

where $\langle \cdot, \cdot \rangle_{\Gamma}$ stands for the duality pairing on $H'(\Gamma) \times H(\Gamma)$ that extends the scalar product on $L^2(\Gamma) \times L^2(\Gamma)$. All the inclusions below are continuous and dense:

$$ H(\Gamma) \subset L^2(\Gamma) \times L^2(\Gamma) \subset H'(\Gamma), $$

and $S_{\Gamma}$ is an isometric operator from $H'(\Gamma)$ onto $H(\Gamma)$. In the following, for every $p \in H(\Gamma)$, we denote $\hat{p} = S_{\Gamma}^{-1} p$ so that $\gamma_{\Gamma} \circ \mathcal{S}_{\Gamma} \hat{p} = p$ and $\|p\|_{H(\Gamma)} = \|\hat{p}\|_{H'(\Gamma)}$.

**Lemma 4.1.** Let $\mathcal{A}$ be the three-dimensional subspace of $H(\Gamma)$ spanned by the total traces of the affine functions in $\mathbb{R}^2$. Then, for every $p, p' \in \mathcal{A}^{\perp}$, $\mathcal{S}_{\Gamma} \hat{p}$ and $\mathcal{S}_{\Gamma} \hat{p}'$ are in $L^2(\mathbb{R}^2)$ and:

$$ (p, p')_{H(\Gamma)} = \int_{\mathbb{R}^2} \Delta(\mathcal{S}_{\Gamma} \hat{p}) \Delta(\mathcal{S}_{\Gamma} \hat{p}') dm. \qquad (13) $$
**Lemma 4.2.** Let $\Gamma$ and $\Gamma'$ be two sets of non-intersecting $C^{1,1}$ Jordan curves. Then for every $p_\Gamma \in H(\Gamma)$ and $p_{\Gamma'} \in H(\Gamma')$, we have:

$$ (\gamma_{\Gamma} \circ \mathcal{S}_{\Gamma'} \hat{p}_{\Gamma'}, S_{\Gamma} \hat{p}_{\Gamma})_{H(\Gamma)} = (S_{\Gamma'} \hat{p}_{\Gamma'}, \gamma_{\Gamma'} \circ \mathcal{S}_{\Gamma} \hat{p}_{\Gamma})_{H(\Gamma')}. $$

*Proof.* Let $\hat{p}_{\Gamma} = (\hat{p}_{\Gamma}^n, \hat{p}_{\Gamma}^d) \in L^2(\Gamma) \times L^2(\Gamma)$ and $\hat{p}_{\Gamma'} = (\hat{p}_{\Gamma'}^n, \hat{p}_{\Gamma'}^d) \in L^2(\Gamma') \times L^2(\Gamma')$. Then:

$$ (\gamma_{\Gamma} \circ \mathcal{S}_{\Gamma'} \hat{p}_{\Gamma'}, S_{\Gamma} \hat{p}_{\Gamma})_{H(\Gamma)} = \int_{\Gamma} \left( \int_{\Gamma'} G(x, y) \hat{p}_{\Gamma'}^n(y) + \partial_{n(y)} G(x, y) \hat{p}_{\Gamma'}^d(y) \right) \mathrm{d}s(y) \hat{p}_{\Gamma}^n(x) \\ + \left( \int_{\Gamma'} \partial_{n(x)} G(x, y) \hat{p}_{\Gamma'}^n(y) + \partial_{n(x)} \partial_{n(y)} G(x, y) \hat{p}_{\Gamma'}^d(y) \right) \mathrm{d}y \hat{p}_{\Gamma}^d(x) \mathrm{d}s(x). $$

Since $\Gamma$ and $\Gamma'$ do not intersect, neither kernel is singular. So, we can reverse the order of integration and conclude with a density argument. $\square$

# 5 The reconstruction problem as a shape-from-moments problem

We return now to Problem 2 with the notation introduced in Section 2. In Definition 9, we choose the constants $\kappa_0$ and $\kappa_1$ in such a way that the operators $S_{\Gamma_0}$ and $S_{\Gamma}$ be strongly elliptic as explained in Theorem 4.1. The solution $\psi$ to System 5 can therefore be represented as a sum of biharmonic single-layer potentials:

$$ \psi = \mathcal{S}_{\Gamma_0} \hat{p}_{\Gamma_0} + \mathcal{S}_{\Gamma} \hat{p}_{\Gamma}. $$

The densities $\hat{p}_{\Gamma_0} = (\hat{p}_{\Gamma_0}^n, \hat{p}_{\Gamma_0}^d) \in H'(\Gamma_0)$ and $\hat{p}_{\Gamma} = (\hat{p}_{\Gamma}^n, \hat{p}_{\Gamma}^d) \in H'(\Gamma)$ satisfy the system:

$$ S_{\Gamma_0} \hat{p}_{\Gamma_0} + \gamma_{\Gamma_0} \circ \mathcal{S}_{\Gamma} \hat{p}_{\Gamma} = f \quad \text{in } H(\Gamma_0) \quad (14a) $$

$$ \gamma_{\Gamma} \circ \mathcal{S}_{\Gamma_0} \hat{p}_{\Gamma_0} + S_{\Gamma} \hat{p}_{\Gamma} = (c, 0) \quad \text{in } H(\Gamma). \quad (14b) $$

In the second identity, $c$ is a function equal to a real constant $c_j$ on each connected component $\Gamma_j$ ($j = 1, \dots, N$) of $\Gamma$. These constants are uniquely determined by the conditions:

$$ \int_{\Gamma_j} \hat{p}_{\Gamma}^n \mathrm{d}s = 0, \quad (j = 1, \dots, N). \quad (14c) $$

For every $j = 1, \dots, N$, let $\mathbf{1}_{\Gamma_j}$ be the piecewise constant function in $H(\Gamma)$, equal to $(1,0)$ on $\Gamma_j$ and $(0,0)$ on $\Gamma_k$ for $k \neq j$ and define the subspaces of $H'(\Gamma)$ and $H(\Gamma)$ of codimension $N$:

$$ \mathcal{H}'(\Gamma) = \{q \in H'(\Gamma) : \langle q, \mathbf{1}_{\Gamma_j} \rangle_{\Gamma} = 0, \quad \forall j = 1, \dots, N\} \quad \text{and} \quad \mathcal{H}(\Gamma) = \{p \in H(\Gamma) : \langle \hat{\mathbf{1}}_{\Gamma_j}, p \rangle = 0, \quad \forall j = 1, \dots, N\}. $$

Notice that identities (14c) mean that $\hat{p}_{\Gamma}$ belongs to $\mathcal{H}'(\Gamma)$ and that the operator $S_{\Gamma}$ isometrically maps $\mathcal{H}'(\Gamma)$ onto $\mathcal{H}(\Gamma)$. We introduce the operators:

$$ K_{\Gamma_0}^{\Gamma} : H(\Gamma_0) \to H(\Gamma) \quad \text{and} \quad K_{\Gamma}^{\Gamma_0} : H(\Gamma) \to H(\Gamma_0) \\ p \mapsto \gamma_{\Gamma} \circ \mathcal{S}_{\Gamma_0} \hat{p} \qquad p \mapsto \gamma_{\Gamma_0} \circ \mathcal{S}_{\Gamma} \hat{p}, \quad (15) $$

and we deduce from Lemma 4.2 that:

$$ (K_{\Gamma_0}^{\Gamma} p_{\Gamma_0}, p_{\Gamma})_{H(\Gamma)} = (p_{\Gamma_0}, K_{\Gamma}^{\Gamma_0} p_{\Gamma})_{H(\Gamma_0)} \quad \text{for all } p_{\Gamma} \in H(\Gamma), p_{\Gamma_0} \in H(\Gamma_0). \quad (16) $$

We denote by $\Pi_{\Gamma}$ the orthogonal projection from $H(\Gamma)$ onto $\mathcal{H}(\Gamma)$ and we apply this operator to equation (14b). Rewriting System (14) in terms of traces instead of densities we obtain:

$$ p_{\Gamma_0} + K_{\Gamma}^{\Gamma_0} p_{\Gamma} = f \quad \text{in } H(\Gamma_0) \quad (17a) $$

$$ \Pi_{\Gamma} \circ K_{\Gamma_0}^{\Gamma} p_{\Gamma_0} + p_{\Gamma} = 0 \quad \text{in } \mathcal{H}(\Gamma). \quad (17b) $$

From the identity (16), specifying that $p_{\Gamma_0}$ is equal to $\mathbf{1}_{\Gamma_0}$, we deduce that if $p_{\Gamma}$ is in $\mathcal{H}(\Gamma)$, then $K_{\Gamma}^{\Gamma_0} p_{\Gamma}$ is in $\mathcal{H}(\Gamma_0)$. It follows from equation (17a) that if $f$ belongs to $\mathcal{H}(\Gamma_0)$, $p_{\Gamma_0}$ belongs to $\mathcal{H}(\Gamma_0)$ as well. We conclude that for every $f \in \mathcal{H}(\Gamma_0)$, there exists $(p_{\Gamma_0}, p_{\Gamma}) \in \mathcal{H}(\Gamma_0) \times \mathcal{H}(\Gamma)$ such that:

$$ p_{\Gamma_0} + K_{\Gamma}^{\Gamma_0} p_{\Gamma} = f \quad \text{in } \mathcal{H}(\Gamma_0) \quad (18a) $$

$$ \Pi_{\Gamma} \circ K_{\Gamma_0}^{\Gamma} p_{\Gamma_0} + p_{\Gamma} = 0 \quad \text{in } \mathcal{H}(\Gamma). \quad (18b) $$

We now rewrite the measurement function (8) in terms of biharmonic single-layer potentials:

$$ \Lambda_{\Gamma} : \mathcal{H}(\Gamma_0) \to \mathcal{H}'(\Gamma_0) \\ f \mapsto (\gamma_{\Gamma_0}^n (\Delta \mathcal{S}_{\Gamma}^{+} \hat{p}_{\Gamma} + \Delta \mathcal{S}_{\Gamma_0}^{-} \hat{p}_{\Gamma_0}), -\gamma_{\Gamma_0}^d (\Delta \mathcal{S}_{\Gamma}^{+} \hat{p}_{\Gamma} + \Delta \mathcal{S}_{\Gamma_0}^{-} \hat{p}_{\Gamma_0})), \quad (19a) $$
where $(p_{\Gamma_0}, p_{\Gamma})$ is the solution of System (18). Defining also:

$$ \Lambda_0 : \mathcal{H}(\Gamma_0) \to \mathcal{H}'(\Gamma_0) \\ f \mapsto (\gamma_{\Gamma_0}^n (\Delta \mathcal{S}_{\Gamma_0}^- \hat{f}), -\gamma_{\Gamma_0}^d (\Delta \mathcal{S}_{\Gamma_0}^- \hat{f})), \qquad (19b) $$

and using the identity $\mathcal{S}_{\Gamma_0}^+ \hat{f} = \mathcal{S}_{\Gamma}^+ \hat{p}_{\Gamma} + \mathcal{S}_{\Gamma_0}^+ \hat{p}_{\Gamma_0}$, we first establish that, for every $f \in \mathcal{H}(\Gamma_0)$:

$$ (\Lambda_{\Gamma} - \Lambda_{0}) f = \left( -[\partial_{n} \Delta \mathcal{S}_{\Gamma_{0}} \hat{p}_{\Gamma_{0}}]_{\Gamma_{0}} + [\partial_{n} \Delta \mathcal{S}_{\Gamma_{0}} \hat{f}]_{\Gamma_{0}}, [\Delta \mathcal{S}_{\Gamma_{0}} \hat{p}_{\Gamma_{0}}]_{\Gamma_{0}} - [\Delta \mathcal{S}_{\Gamma_{0}} \hat{f}]_{\Gamma_{0}} \right). $$

Then, from the relations (12), we deduce that $S_{\Gamma_0} \circ (\Lambda_\Gamma - \Lambda_0) f = p_{\Gamma_0} - f$ and therefore, denoting by $R_\Gamma$ the operator $S_{\Gamma_0} \circ (\Lambda_\Gamma - \Lambda_0)$, we finally obtain that:

$$ (\mathrm{Id} + R_{\Gamma}) f = p_{\Gamma_0}. \qquad (20a) $$

On the other hand, according to (18a), we have $p_{\Gamma_0} = f - K_{\Gamma_0}^{\Gamma_0} p_{\Gamma}$ and applying the operator $K_{\Gamma_0}^{\Gamma_0}$ to (18b) yields $K_{\Gamma_0}^{\Gamma_0} p_{\Gamma} = -K_{\Gamma_0}^{\Gamma_0} \circ \Pi_{\Gamma} \circ K_{\Gamma_0}^{\Gamma_0} p_{\Gamma_0}$. Therefore, denoting $K_{\Gamma} = K_{\Gamma_0}^{\Gamma_0} \circ \Pi_{\Gamma} \circ K_{\Gamma_0}^{\Gamma_0}$ we get:

$$ (\mathrm{Id} - K_{\Gamma}) p_{\Gamma_0} = f. \qquad (20b) $$

The theorem below follows directly from the identities (20a) and (20b):

**Theorem 5.1.** The operators $(\mathrm{Id} + R_\Gamma) : \mathcal{H}(\Gamma_0) \to \mathcal{H}(\Gamma_0)$ and $(\mathrm{Id} - K_\Gamma) : \mathcal{H}(\Gamma_0) \to \mathcal{H}(\Gamma_0)$ are invertible and they are each other's inverse. Moreover:

$$ K_{\Gamma} = R_{\Gamma} \circ (\mathrm{Id} + R_{\Gamma})^{-1}. \qquad (21) $$

From the point of view of solving the inverse problem we are dealing with, it is important to note that the right-hand member in the identity (21) can be computed from $\Lambda_\Gamma$ and $\Lambda_0$, two operators we assume to be known. We now turn our attention to the operator $K_\Gamma$. From equality (16), we deduce straightforwardly:

**Theorem 5.2.** For every $f, g \in \mathcal{H}(\Gamma_0)$:

$$ (K_{\Gamma} f, g)_{H(\Gamma_0)} = (\Pi_{\Gamma} \circ K_{\Gamma_0}^{\Gamma_0} f, K_{\Gamma_0}^{\Gamma_0} g)_{H(\Gamma)}. \qquad (22) $$

This theorem is the keystone of the reconstruction method. Indeed, let $F$ et $G$ be biharmonic functions in $H^2(\Omega)$, denote $f_{\Gamma_0} = \gamma_{\Gamma_0} F$ and $g_{\Gamma_0} = \gamma_{\Gamma_0} G$ and assume that $f_{\Gamma_0}$ and $g_{\Gamma_0}$ are in $\mathcal{H}(\Gamma_0)$ (simply add an appropriate constant to $F$ and $G$ to satisfy this assumption). Since $\mathcal{S}_{\Gamma_0} \hat{f}_{\Gamma_0} = F$ and $\mathcal{S}_{\Gamma_0} \hat{g}_{\Gamma_0} = G$ in $\Omega$, we deduce from (22) that:

$$ (K_{\Gamma} f_{\Gamma_0}, g_{\Gamma_0})_{H(\Gamma_0)} = (\Pi_{\Gamma} f_{\Gamma}, g_{\Gamma})_{H(\Gamma)}, \qquad (23) $$

where $f_\Gamma = \gamma_\Gamma F$ and $g_\Gamma = \gamma_\Gamma G$. In other words, we have access to the scalar product in $\mathcal{H}(\Gamma)$ of the traces (up to an additive constant because of the projector $\Pi_\Gamma$) of any two functions $F$ and $G$ that are biharmonic in $\Omega$. This will enable us to calculate the complex moments of $\mathcal{O}$.

Let us denote by $\mathcal{H}(\Omega)$ the space of the harmonic functions in $H^2(\Omega)$.

**Lemma 5.1.** Let $p$ be in $\mathcal{H}(\Gamma)$ such that $(K_{\Gamma_0}^{\Gamma_0} p, \gamma_{\Gamma_0} h)_{H(\Gamma_0)} = 0$ for every function $h$ in $\mathcal{H}(\Omega)$. Then $\mathcal{S}_{\Gamma} \hat{p}$ is harmonic outside $\mathcal{O}$.

*Proof.* Let $D_\Omega$ stands for the Dirichlet-to-Neumann operator on $I_0^-$ (see Section A in the Appendix), let $p$ be in $\mathcal{H}(\Gamma)$ as in the statement of the Lemma and denote $g_{\Gamma_0} = K_{\Gamma_0}^{\Gamma_0} p$. The hypothesis means that:

$$ -\langle [\partial_n \Delta \mathcal{S}_{\Gamma_0} \hat{g}_{\Gamma_0}]_{\Gamma_0}, q \rangle_{-\frac{3}{2}, \frac{3}{2}} + \langle [\Delta \mathcal{S}_{\Gamma_0} \hat{g}_{\Gamma_0}]_{\Gamma_0}, D_{\Omega} q \rangle_{-\frac{1}{2}, \frac{1}{2}} = 0 \quad \text{for all } q \in H^{3/2}(\Gamma_0). $$

The operator $D_\Omega$ being self-adjoint (see Proposition A.1), it follows that:

$$ \langle \gamma_{\Gamma_0}^n \Delta \mathcal{S}_{\Gamma_0}^- \hat{g}_{\Gamma_0} - D_{\Omega} \Delta \mathcal{S}_{\Gamma_0}^- \hat{g}_{\Gamma_0}, q \rangle_{-\frac{3}{2}, -\frac{3}{2}} - \langle \gamma_{\Gamma_0}^n \Delta \mathcal{S}_{\Gamma_0}^+ \hat{g}_{\Gamma_0} - D_{\Omega} \Delta \mathcal{S}_{\Gamma_0}^+ \hat{g}_{\Gamma_0}, q \rangle_{-\frac{1}{2}, \frac{1}{2}} = 0 \quad \text{for all } q \in H^{3/2}(\Gamma_0). $$

The first term vanishes by definition of $D_\Omega$, which allows us to deduce that $\gamma_{\Gamma_0}^n \Delta \mathcal{S}_{\Gamma_0}^+ \hat{g}_{\Gamma_0} = D_\Omega \Delta \mathcal{S}_{\Gamma_0}^+ \hat{g}_{\Gamma_0}$. Let $u_\Omega$ be the function in $L^2(\Omega)$, harmonic and such that $\gamma_{\Gamma_0}^d u_\Omega = \gamma_{\Gamma_0}^d \Delta \mathcal{S}_{\Gamma_0}^+ \hat{g}_{\Gamma_0}$ (existence and uniqueness of such a function is asserted in Proposition B.1). The function $u$ defined by $u = u_\Omega$ in $\Omega$ and $u = \Delta \mathcal{S}_{\Gamma_0}^+ \hat{g}_{\Gamma_0}$ in $\mathbb{R}^2 \setminus \bar{\Omega}$ is therefore harmonic in $\mathbb{R}^2$. According to [17, Lemma 4.1], $u(x) = \mathcal{O}(\ln|x|)$ as $|x|$ goes to $+\infty$, which together with [3, Theorem 9.10] (generalized Liouville Theorem) implies that $u$ is constant in $\mathbb{R}^2$. Referring again to [17, Lemma 4.1], the only possible constant is zero. Since $\mathcal{S}_{\Gamma}^+ p = \mathcal{S}_{\Gamma_0}^+ \hat{g}_{\Gamma_0}$ in $\mathbb{R}^2 \setminus \bar{\Omega}$, we conclude that $\Delta \mathcal{S}_{\Gamma}^+ p = 0$ in $\mathbb{R}^2 \setminus \bar{\mathcal{O}}$. $\square$

**Theorem 5.3.** Let $h$ be a harmonic function in $L^2(\Omega)$. For every $\varepsilon > 0$ there exists $f_{\Gamma_0}^\varepsilon$ in $\mathcal{H}(\Gamma_0)$ such that

$$ \left| (K_{\Gamma} f_{\Gamma_0}^{\varepsilon}, f_{\Gamma_0}^{\varepsilon})_{H(\Gamma_0)} - \int_{\mathcal{O}} h^2 dm \right| < \varepsilon. $$
Proof. Let $g_{\Gamma_0}$ be in $H(\Gamma_0)$ such that $\Delta \mathcal{S}_{\Gamma_0}^- \hat{g}_{\Gamma_0} = h$ (for instance choose $g_{\Gamma_0} = \gamma_{\Gamma_0} u$ with $u$ such that $\Delta u = h$ in $H^2(\Omega) \cap H_0^1(\Omega)$). Let $g_{\Gamma_0}^h$ be the orthogonal projection in $H(\Gamma)$ of $K_{\Gamma_0}^\Gamma g_{\Gamma_0}$ onto $\overline{K_{\Gamma_0}^\Gamma \mathcal{H}(\Omega)}$. Thus, for every $\varepsilon > 0$, there exists $h^\varepsilon \in \mathcal{H}(\Omega)$ such that

$$ \|g_{\Gamma}^{h} - K_{\Gamma_{0}}^{\Gamma} \circ \gamma_{\Gamma_{0}} h^{\varepsilon}\|_{H(\Gamma_{0})} < \varepsilon [1 + 2 \|K_{\Gamma_{0}}^{\Gamma} g_{\Gamma_{0}} - g_{\Gamma}^{h}\|_{H(\Gamma)}]^{-1}. \quad (24) $$

Define $f_{\Gamma_0}^\varepsilon = \Pi_{\Gamma_0}(g_{\Gamma_0} - \gamma_{\Gamma_0}h^\varepsilon)$ where $\Pi_{\Gamma_0}$ is the orthogonal projection onto $\mathcal{H}(\Gamma_0)$ in $H(\Gamma_0)$. We have:

$$ (K_{\Gamma} f_{\Gamma_0}^{\varepsilon}, f_{\Gamma_0}^{\varepsilon})_{H(\Gamma_0)} = (\Pi_{\Gamma} \circ K_{\Gamma_0}^{\Gamma} f_{\Gamma_0}^{\varepsilon}, K_{\Gamma_0}^{\Gamma} f_{\Gamma_0}^{\varepsilon})_{H(\Gamma)} = (\Pi_{\Gamma} \circ K_{\Gamma_0}^{\Gamma} f_{\Gamma_0}^{\varepsilon}, \Pi_{\Gamma} \circ K_{\Gamma_0}^{\Gamma} f_{\Gamma_0}^{\varepsilon})_{H(\Gamma)}. $$

Notice that, for every $p \in H(\Gamma_0)$, $\Pi_{\Gamma} \circ K_{\Gamma_0}^{\Gamma} \circ \Pi_{\Gamma_0} p = \Pi_{\Gamma} \circ K_{\Gamma_0}^{\Gamma} p$ because the functions $\mathcal{S}_{\Gamma_0}p$ and $\mathcal{S}_{\Gamma_0} \circ \Pi_{\Gamma_0} p$ differ only up to a constant in $\Omega$. This implies that:

$$ (K_{\Gamma} f_{\Gamma_0}^{\varepsilon}, f_{\Gamma_0}^{\varepsilon})_{H(\Gamma_0)} = \| \Pi_{\Gamma} \circ K_{\Gamma_0}^{\Gamma} (g_{\Gamma_0} - \gamma_{\Gamma_0} h^{\varepsilon}) \|_{H(\Gamma)}^2. \quad (25a) $$

On the other hand, observing that $\Pi_{\Gamma}(K_{\Gamma_0}^{\Gamma} g_{\Gamma_0} - g_{\Gamma}^{h}) = K_{\Gamma_0}^{\Gamma} g_{\Gamma_0} - g_{\Gamma}^{h}$, we get:

$$ \begin{aligned} \left\| \Pi_{\Gamma} \circ K_{\Gamma_0}^{\Gamma} (g_{\Gamma_0} - \gamma_{\Gamma_0} h^{\varepsilon}) \right\|_{H(\Gamma)}^2 - \left\| K_{\Gamma_0}^{\Gamma} g_{\Gamma_0} - g_{\Gamma}^{h} \right\|_{H(\Gamma)}^2 &\le \| \Pi_{\Gamma} (g_{\Gamma}^{h} - K_{\Gamma_0}^{\Gamma} \circ \gamma_{\Gamma_0} h^{\varepsilon}) \|_{H(\Gamma)} \left[ \| \Pi_{\Gamma} (g_{\Gamma}^{h} - K_{\Gamma_0}^{\Gamma} \circ \gamma_{\Gamma_0} h^{\varepsilon}) \|_{H(\Gamma)} + 2 \| K_{\Gamma_0}^{\Gamma} g_{\Gamma_0} - g_{\Gamma}^{h} \|_{H(\Gamma)} \right], \end{aligned} $$

which, with the estimate (24) gives for every $\varepsilon < 1$:

$$ \left\| \Pi_{\Gamma} \circ K_{\Gamma_0}^{\Gamma} (g_{\Gamma_0} - \gamma_{\Gamma_0} h^{\varepsilon}) \right\|_{H(\Gamma)}^2 - \left\| K_{\Gamma_0}^{\Gamma} g_{\Gamma_0} - g_{\Gamma}^{h} \right\|_{H(\Gamma)}^2 < \varepsilon. \quad (25b) $$

By construction:

$$ (K_{\Gamma_0}^{\Gamma} g_{\Gamma_0} - g_{\Gamma}^{h}, K_{\Gamma_0}^{\Gamma} \circ \gamma_{\Gamma_0} h)_{H(\Gamma)} = 0, $$

but also, according to (16), for every $h \in \mathcal{H}(\Omega)$:

$$ (K_{\Gamma_0}^{\Gamma} g_{\Gamma_0} - g_{\Gamma}^{h}, K_{\Gamma_0}^{\Gamma} \circ \gamma_{\Gamma_0} h)_{H(\Gamma)} = (K_{\Gamma}^{\Gamma_0} (K_{\Gamma_0}^{\Gamma} g_{\Gamma_0} - g_{\Gamma}^{h}), \gamma_{\Gamma_0} h)_{H(\Gamma_0)}. $$

This implies, with Lemma 5.1 that $\mathcal{S}_{\Gamma}(K_{\Gamma_0}^{\Gamma} g_{\Gamma_0} - g_{\Gamma}^{h})$ is harmonic outside $\mathcal{O}$. On the other hand, inside $\mathcal{O}$, $\mathcal{S}_{\Gamma} \circ K_{\Gamma_0}^{\Gamma} g_{\Gamma_0} = \mathcal{S}_{\Gamma_0} g_{\Gamma_0}$ and $\mathcal{S}_{\Gamma} g_{\Gamma_0}^h$ is harmonic (because $g_{\Gamma}^h$ is in the space $\overline{K_{\Gamma_0}^\Gamma \mathcal{H}(\Omega)}$) and therefore $\Delta \mathcal{S}_{\Gamma}(K_{\Gamma_0}^{\Gamma} g_{\Gamma_0} - g_{\Gamma}) = h$. According to (13), we have:

$$ \| K_{\Gamma_0}^{\Gamma} g_{\Gamma_0} - g_{\Gamma}^{h} \|_{H(\Gamma)}^2 = \int_{\mathcal{O}} h^2 dm. \quad (25c) $$

Combining the equations (25) leads to the conclusion. $\square$

Using the polarization identity, we prove:

**Corollary 5.1.** For every $\varepsilon > 0$ and every $k, m \in \mathbb{N}$, there exist $f_k^\varepsilon$ and $f_m^\varepsilon$ in the complex space $\mathcal{H}(\Gamma_0) + i\mathcal{H}(\Gamma_0)$ such that:

$$ \left| (K_{\Gamma} f_{k}^{\varepsilon}, f_{m}^{\varepsilon})_{H(\Gamma_{0})} - \int_{\mathcal{O}} z^{k} \bar{z}^{m} dm(z) \right| < \varepsilon. $$

We have now reached our goal: turning the reconstruction problem 2 into a shape-from-moments problem. Indeed, Theorem 5.1 explains how to compute the operator $K_{\Gamma}$ from the measurement operator $\Lambda_{\Gamma}$ and Corollary 5.1 shows that a suitable choice of inputs can be used to calculate with an arbitrary precision the complex moments of the obstacles.

# 6 Using a finite section of all the complex moments

For a positive integer $n$ we assume known the complex moments of $\mathcal{O}$:

$$ \int_{\mathcal{O}} \overline{z^k} z^l dm(z) \quad \text{for all } k \le n, l \le n. $$

Applying a Gram-Schmidt process to the family of monomials $\{1, z, z^2, \dots, z^n\}$ we can compute from these moments the so-called Bergman polynomials $P_0, P_1, \dots, P_n$, i.e. the polynomials orthonormalized for the scalar product of $L^2(\mathcal{O})$ and such that $P_k$ is of degree $k$ for every $k = 0, \dots, n$. Following the ideas developed in [11], we introduce the function $\Theta_n$ defined by:

$$ \Theta_n(z) = \frac{1}{\sqrt{\pi \sum_{j=0}^{n} |P_j(z)|^2}} \quad \text{for all } z \in \mathbb{C}. \quad (26) $$

When all the curves $\Gamma_j$ ($j = 1, \dots, n$) are analytic, this function is shown to approximate the distance to $\Gamma$ in $\mathcal{O}$ while it decays to zero at certain rates, as $n$ goes to $+\infty$ on $\Gamma$ and in $\mathbb{C} \setminus \overline{\mathcal{O}}$. More precisely, $\text{dist}(z, \Gamma) \le \Theta_n(z)$ for $z \in \mathcal{O}$ and there exists two positive constants $C_1$ and $C_2$ such that $C_1/n \le \Theta_n(z) \le C_2/n$ for $z \in \Gamma$. We will see in Section 8 that plotting the level sets $\{z \in \mathbb{C} : \Theta_n(z) = \lambda/n\}$ for certain values of $\lambda$ gives a pretty good approximation of the boundaries of the obstacles.
# 7 Using a finite section of the harmonic moments only

For a positive integer $n$ we assume known the harmonic moments:

$$ \tau_{\ell} = \int_{\mathcal{O}} z^{\ell} dm(z) \quad \text{for all } \ell = 0, \dots, 2n-1. $$

An algorithm of reconstruction of $\mathcal{O}$ from the $\tau_{\ell}$ is detailed in [7]. The key idea is to approximate $\mathcal{O}$ by so-called quadrature domains. The method consists of two stages:

1. Determine complex weights $c_j$ and complex nodes $z_j$ ($j = 1, \dots, n$) such that

$$ \sum_{j=0}^{n} c_j z_j^{\ell} = \tau_j \quad \text{for all } \ell = 0, \dots, 2n-1. \qquad (27a) $$

Such a system of equations is called a Prony's system.

2. Determine a domain $\mathcal{O}_n$ satisfying the quadrature identity:

$$ \int_{\mathcal{O}_n} z^{\ell} dm(z) = \sum_{j=0}^{n} c_j z_j^{\ell} \quad \text{for all } \ell \in \mathbb{N}. \qquad (27b) $$

By construction, $\mathcal{O}$ satisfies the identities above for $\ell = 0, \dots, 2n-1$ but we emphasize that $\mathcal{O}_n$ is required to satisfy these identities for all the integers $\ell$. Such a domain is called a quadrature domain.

Regarding the first stage, we introduce the polynomial:

$$ P_n(z) = \begin{vmatrix} \tau_0 & \tau_1 & \cdots & \tau_{n-1} & \tau_n \\ \tau_1 & \tau_2 & \cdots & \tau_n & \tau_{n+1} \\ \vdots & \vdots & & \vdots & \vdots \\ \tau_{n-1} & \tau_n & \cdots & \tau_{2n-2} & \tau_{2n-1} \\ 1 & z & \cdots & z^{n-1} & z^n \end{vmatrix}, $$

which enters the statement of the following result proved in [7]:

**Theorem 7.1.** Equations (27a) admit a solution if and only if the polynomial $P_n$ admits $n$ simple roots. In this case, this solution is unique and the nodes $z_1, \dots, z_n$ are the roots of $P_n$.

The explicit determination of the nodes $z_j$ and weights $c_j$ ($j = 1, \dots, n$) is carried out by means of the matrix pencil method described in [9]. Introducing the Hankel matrices:

$$ \mathbb{H}_0^{(n)} = \begin{pmatrix} \tau_0 & \tau_1 & \cdots & \tau_{n-1} \\ \tau_1 & \tau_2 & \cdots & \tau_n \\ \vdots & \vdots & \ddots & \vdots \\ \tau_{n-1} & \tau_n & \cdots & \tau_{2n-2} \end{pmatrix} \quad \text{and} \quad \mathbb{H}_1^{(n)} = \begin{pmatrix} \tau_1 & \tau_2 & \cdots & \tau_n \\ \tau_2 & \tau_3 & \cdots & \tau_{n+1} \\ \vdots & \vdots & \ddots & \vdots \\ \tau_n & \tau_{n+1} & \cdots & \tau_{2n-1} \end{pmatrix}, \qquad (28a) $$

the nodes $z_j$ are the solutions of the following generalized eigenvalue problem:

$$ \exists \xi \in \mathbb{C}^n \setminus \{0\}, \quad \mathbb{H}_0^{(n)} \xi = z \mathbb{H}_1^{(n)} \xi. \qquad (28b) $$

If the $z_j$ are pairwise distinct, the weights $c_j$ are obtained by solving the Vandermonde linear system:

$$ \begin{pmatrix} 1 & 1 & \cdots & 1 \\ z_1 & z_2 & \cdots & z_n \\ z_1^2 & z_2^2 & \cdots & z_n^2 \\ \vdots & \vdots & & \vdots \\ z_1^{n-1} & z_2^{n-1} & \cdots & z_n^{n-1} \end{pmatrix} \begin{pmatrix} c_1 \\ c_2 \\ c_3 \\ \vdots \\ c_n \end{pmatrix} = \begin{pmatrix} \tau_0 \\ \tau_1 \\ \tau_2 \\ \vdots \\ \tau_{n-1} \end{pmatrix}. \qquad (28c) $$

One can verify that $z_1, \dots, z_n$ and $c_1, \dots, c_n$ solve (27a) (see [9, Section 3]). Since the inverse problem we are dealing with is ill-posed, it is not surprising that it leads to a numerical method requiring the solution of two ill-conditioned problems: a generalized eigenvalue problem with Hankel matrices and a Vandermonde system. It is worth mentioning also that numerical instabilities increase with $n$.

We now turn our attention to the second step of the reconstruction method. Disks are the simplest examples of quadrature domains, and can be shown to be the only ones for which $n = 1$. Rather counter-intuitively, there are many quadrature domains. Actually, every domain whose boundary consists in non-intersecting $C^{\infty}$ Jordan curves is arbitrarily close to a quadrature domain.
Given a set of nodes and weights, existence of a quadrature domain satisfying (27b) is asserted in [10, Theorem 2.4, (vi)], providing that the weights are real and positive (uniqueness does not hold in general). The construction of this domain is strongly related to partial balayage of measures (see [12]) and free boundary problems (see [13]). In summary, if the weights $c_j$ ($j = 1, \dots, n$) are positive in (27b), a quadrature domain $O_n$ satisfying (27b) can be obtained by means of the following identity:

$$1_{O_n} = -\Delta V^n, \qquad (29a)$$

where $V^n$ is the unique function achieving:

$$\min_{v \in K^n} \frac{1}{2} \int_{\Omega} |\nabla v|^2 dm - \int_{\Omega} v dm, \qquad (29b)$$

with:

$$U^n = -\frac{1}{2\pi} \sum_{j=1}^{n} c_j \ln |\cdot - z_j| \quad \text{and} \quad K^n = \{v \in H^1(\Omega) : \gamma_{\Gamma_0}^d v = \gamma_{\Gamma_0}^d U^n \text{ and } v \le U^n \text{ in } \Omega\}. \qquad (29c)$$

The function that associates to any set of nodes and weights the corresponding quadrature domain (by means of the steps (29)) is called the partial balayage operator.

# 8 Algorithm and numerical tests

In this section we describe an algorithm derived from the results of Sections 6 and 7 and we provide some numerical tests to illustrate the efficiency of the reconstruction methods.

## Generation of the measurements

Fix a positive integer $m \ge 3$, denote $m' = 2m - 1$ and for $k = 1, \dots, m'$ define the functions:

$$F_k(z) = \begin{cases} z^k & \text{if } k = 1, \dots, m \\ \bar{z}^{k-m}/(4(k-m)) & \text{if } k = m+1, \dots, m' \end{cases} \quad \text{for all } z \in \mathbb{C}.$$

Notice that $\Delta \bar{z}^j/(4j) = z^{j-1}$ for every $j \ge 1$ and hence $\Delta^2 F_k = 0$ for every $k = 1, \dots, m'$. Our first objective is to generate measurements corresponding to the boundary data $f_k = \tilde{\gamma}_{\Gamma_0} F_k$ where $\tilde{\gamma}_{\Gamma_0} = \Pi_{\Gamma_0} \circ \gamma_{\Gamma_0}$.

Computations are based on a BEM involving biharmonic single-layer potentials. For every $k = 1, \dots, m'$, we compute the complex densities $q_j^k$ ($j = 0, \dots, n$) that solve the system of integral equations:

$$\gamma_{\Gamma_0} [ \mathcal{S}_{\Gamma_0} q_0^k + \mathcal{S}_{\Gamma_1} q_1^k + \dots + \mathcal{S}_{\Gamma_N} q_N^k ] = \gamma_{\Gamma_0} F_k \quad \text{on } \Gamma_0 \qquad (30a)$$

$$\gamma_{\Gamma_j} [ \mathcal{S}_{\Gamma_0} q_0^k + \mathcal{S}_{\Gamma_1} q_1^k + \dots + \mathcal{S}_{\Gamma_N} q_N^k ] = 0 \quad \text{on } \Gamma_j \quad (j = 1, \dots, N). \qquad (30b)$$

The function $\psi_k = \sum_{j=0}^N \mathcal{S}_{\Gamma_j} q_j^k$ is not the solution we are looking for as it satisfies neither the conditions $\gamma_{\Gamma_0} \psi = f_k$ nor the boundary condition (5c) (such that (5e) holds). So, we compute for every $j, k = 0, \dots, N$, the densities $p_j^k$ such that:

$$\gamma_{\Gamma_j} [ \mathcal{S}_{\Gamma_0} p_0^k + \mathcal{S}_{\Gamma_1} p_1^k + \dots + \mathcal{S}_{\Gamma_N} p_N^k ] = (\delta_{jk}, 0) \quad \text{on } \Gamma_j \quad (j = 0, 1, \dots, N).$$

It follows that $\xi_k = \sum_{j=0}^N \mathcal{S}_{\Gamma_j} p_j^k$ is the function introduced in Section 2 and that solves System (6). Then, we determine the constants $\alpha_{j,k}$ ($j = 0, \dots, N, k = 1, \dots, m'$) by inverting the linear system:

$$[ \langle \hat{1}_{\Gamma_j}, q_j^k \rangle ]_{\substack{0 \le j \le N \\ 1 \le k \le m'}} = [ \langle \hat{1}_{\Gamma_j}, p_j^k \rangle ]_{\substack{0 \le j \le N \\ 0 \le k \le N}} [ \alpha_{j,k} ]_{\substack{0 \le j \le N \\ 1 \le k \le m'}} ,$$

and we define $\tilde{q}_j^k = q_j^k - \sum_{i=0}^N \alpha_{i,k} p_i^k$. This time, the function $\tilde{\psi}_k = \sum_{j=0}^N \mathcal{S}_{\Gamma_j} \tilde{q}_j^k$ satisfies all the conditions mentioned above. Recall that we assume we have the operator $\Lambda_{\Gamma}$ (defined in (8)) at our disposal. Since the definition of the operator $\Lambda_0$ depends only on the boundary $\Gamma_0$, we have access, as already mentioned, to the operator $R_{\Gamma} = S_{\Gamma_0} \circ (\Lambda_{\Gamma} - \Lambda_0)$ (see Theorem 5.1) and then also to the operator $V_{\Gamma} = \text{Id} + R_{\Gamma}$. We easily verify that:

$$V_{\Gamma} f_k = \tilde{q}_0^k \quad \text{for all } k = 1, \dots, m',$$

and this will be the measurements we shall use for the reconstruction.
## Computing the scalar product on Γ

The objective of this subsection is to explain how to compute the scalar products:

$$ (\tilde{\gamma}_{\Gamma} \overline{F_j}, \tilde{\gamma}_{\Gamma} F_k)_{H(\Gamma)} \quad \text{for all } j, k = 1, \dots, m'. $$

According to Theorem 5.1, $K_Γ = \text{Id} - V_Γ^{-1}$ and it follows that:

$$ (K_Γ \overline{f}_j, f_k)_{H(Γ_0)} = (\overline{f}_j, f_k)_{H(Γ_0)} - (V_Γ^{-1} \overline{f}_j, f_k)_{H(Γ_0)} \quad \text{for all } j, k = 1, \dots, m'. \quad (31) $$

Notice that, by definition

$$ (\overline{f}_j, f_k)_{H(Γ_0)} = (\tilde{\gamma}_{Γ_0} \overline{F}_j, \tilde{\gamma}_{Γ_0} F_k)_{H(Γ_0)}, $$

and according to Theorem 5.2, for every $j, k = 1, \dots, m'$ we have:

$$ (K_Γ \overline{f}_j, f_k)_{H(Γ_0)} = (\Pi_Γ \circ K_{Γ_0}^Γ \overline{f}_j, K_{Γ_0}^Γ f_k)_{H(Γ)} = (\tilde{\gamma}_Γ \overline{F}_j, \tilde{\gamma}_Γ F_k)_{H(Γ)}. $$

We define the positive-definite symmetric matrices:

$$ Q_m^{\Gamma_0} = [(\overline{f}_j, f_k)_{H(Γ_0)}]_{\substack{1 \le j \le m' \\ 1 \le k \le m'}} \quad \text{and} \quad Q_m^{\Gamma} = [(K_Γ \overline{f}_j, f_k)_{H(Γ)}]_{\substack{1 \le j \le m' \\ 1 \le k \le m'}} $$

and also:

$$ V_m = [(V_Γ \overline{f}_j, f_k)_{H(Γ_0)}]_{\substack{1 \le j \le m' \\ 1 \le k \le m'}} \quad \text{and} \quad W_m = [(V_Γ^{-1} \overline{f}_j, f_k)_{H(Γ_0)}]_{\substack{1 \le j \le m' \\ 1 \le k \le m'}} $$

The matrix $W_m$ is approximated by the matrix $Q_m^{\Gamma_0} V_m^{-1} Q_m^{\Gamma_0}$ so that (31) can be rewritten in matrix form:

$$ Q_m^{\Gamma} = Q_m^{\Gamma_0} - Q_m^{\Gamma_0} V_m^{-1} Q_m^{\Gamma_0}. $$

## Computing the complex moments of the obstacles

According to Lemma 5.1, if $p$ is in $\mathcal{H}(\Gamma)$ and satisfies $(p, K_{\Gamma_0}^\Gamma \circ \gamma_{\Gamma_0} h)_{H(\Gamma)} = 0$ for every $h \in \mathcal{H}(\Omega)$, then $\Delta \mathcal{S}_\Gamma^+ \hat{p} = 0$ and in this case:

$$ \| \mathcal{S}_\Gamma \hat{p} \|_{H(\Gamma)}^2 = \| \Delta \mathcal{S}_\Gamma \hat{p} \|_{L^2(\mathcal{O})}^2. $$

To put these principles into practice, let us denote $G_k = F_{m+k}$ and recall that $\Delta G_k = z^{k-1}$ for $k = 1, \dots, m-1$. We decompose the matrix $Q_m^\Gamma$ into 4 sub-matrices:

$$ Q_m^{\Gamma} = \begin{bmatrix} X_m & Y_m \\ Y_m^* & Z_m \end{bmatrix}, $$

where, by construction:

$$ X_m = [(\tilde{\gamma}_\Gamma \overline{z^j}, \tilde{\gamma}_\Gamma z^k)_{H(\Gamma)}]_{\substack{1 \le j \le m \\ 1 \le k \le m}} \quad Y_m = [(\tilde{\gamma}_\Gamma \overline{z^j}, \tilde{\gamma}_\Gamma G_k)_{H(\Gamma)}]_{\substack{1 \le j \le m \\ 1 \le k \le m-1}} \quad \text{and} \quad Z_m = [(\tilde{\gamma}_\Gamma \overline{G_j}, \tilde{\gamma}_\Gamma G_k)_{H(\Gamma)}]_{\substack{1 \le j \le m-1 \\ 1 \le k \le m-1}}. $$

It follows that the entries of the matrix $M_m^\Gamma = Z_m - Y_m^* X_m^{-1} Y_m$ are

$$ (\Pi_m \circ \gamma_\Gamma \overline{G}_j, \Pi_m \circ \gamma_\Gamma G_k)_{H(\Gamma)} \quad \text{for all } j, k = 1, \dots, m-1, $$

where $\Pi_m$ is the orthogonal projection in $H(\Gamma)$ onto the orthogonal of the subspace spanned by $\{1, z, z^2, \dots, z^m\}$. As a conclusion, for $m$ large enough, the entries of the matrix $M_m^\Gamma$ are such that:

$$ (\Pi_m \circ \gamma_\Gamma \overline{G}_j, \Pi_m \circ \gamma_\Gamma G_k)_{H(\Gamma)} \approx \int_{\mathcal{O}} \overline{z^{j-1} z^{k-1}} \, dm(z) \quad \text{for all } j, k = 1, \dots, m-1. $$

## Reconstruction with a finite section of all the complex moments

Let $L_m$ be the upper triangular matrix entering the Cholevsky factorization of $M_m^\Gamma$ i.e. $M_m^\Gamma = L_m^* L_m$. Then, it is easy to verify that the columns of $L_m$ are the coefficients of the Bergman polynomials. It is then possible to evaluate the function $\Theta_{m-2}$ (defined in (26)) at any point of $\Omega$ and draw the set corresponding to $\Theta_{m-2} = \lambda/(m-2)$ for some values of $\lambda$.

## Reconstruction with a finite section of the harmonic moments only

With this approach, we use only the first row of the matrix $M_m^\Gamma$. We fix an integer $n \le (m-1)/2$ and we apply the method described in Section 7. As explained there, to solve the Prony's system (27a) we construct the Hankel matrices (28a) and solve the generalized eigenvalue problem (28b) which provides the values of the nodes $z_j$. The weights $c_j$ are obtained by solving the linear system (28c). We replace $c_j$ by $\text{Re}(c_j)$ when $\text{Im}(c_j)$ is non-zero and we solve the convex minimization problem (29) (partial balayage process) to obtain (hopefully) approximations of the obstacles.
## Numerical tests

For all the tests, the domain $\Omega$ is a disk centered at the origin and of radius 1.

**Example 1.** We consider a single cross-shaped obstacle whose boundary is parameterized as:

$$ \begin{pmatrix} x_1(\theta) \\ x_2(\theta) \end{pmatrix} = \left[ 0.25(1 + 0.4 \cos(4\theta)) \begin{pmatrix} \cos(\theta) \\ \sin(\theta) \end{pmatrix} + 0.2 \begin{pmatrix} 1 \\ 1 \end{pmatrix} \right], \quad (\theta \in [0, 2\pi]). $$

Note that this is a non-trivial example (neither small, nor too close to the boundary $\Gamma_0$, nor a ball, nor even convex). The reconstruction is performed with $m = 13$ (corresponding to $2m - 1 = 25$ measurements). It means that we have at our disposal the moments:

$$ \int_{\mathcal{O}} \bar{z}^k z^n dm(z) \quad \text{for all } k, n = 0, \dots, 11. $$

Numerical instabilities make it difficult to increase $m$ much further in this case (for larger $m$ the matrix $M_m^\Gamma$ exhibits spurious negative eigenvalues). First we apply the method described in Section 6. Several level lines corresponding to $\Theta_{11} = \lambda/(m-2)$ for $\lambda = 0, 2, 0, 4, 0.6, 0.8$ and $1$ are plotted on Fig. 2. The position of the obstacle is accurately found and the level lines outline clearly its convex hull. However, the details of the branches are not captured by the reconstruction. This would probably require $m$ to be larger.

<figure><img src="image_4.png" /><figcaption>Figure 2: Reconstruction of a single cross-shaped obstacle using the method described in Section 6 (borrowed from [11]) with 25 measurements (i.e. $m = 13$). The image on the right is a close-up of the obstacle and its reconstruction.</figcaption></figure>

We apply now the method explained in Section 7 and for which only the harmonic moments:

$$ \int_{\mathcal{O}} z^k dm(z) \quad \text{for all } k = 0, \dots, 10, $$

are involved. We solve the related Prony's system (27a) for $n = 1, 2, \dots, 5$ and for each case, we represent in Fig. 3 the disks centered at the nodes $z_j$ and of radii $\sqrt{\text{Re}(c_j)/\pi}$. Indeed, when these disk are pairwise disjoint, the mean value property for harmonic functions ensures that their union is a quadrature domain that satisfies equality (27b). For $n = 1, 2, 3$, only one weight is non zero and hence only one disk is non-degenerated (obviously no partial balayage step is needed). For $n = 4$, we obtain four disjoint disks, so no partial balayage step is necessary either in this case. On the other hand, the disks overlap for $n = 5$, so we solve the convex minimization problem (29). The resulting quadrature domain is represented on Fig. 4. The reconstruction is much better than on Fig. 2, as the outline of the cross can be clearly seen. Such reconstruction accuracy is remarkable for this type of problem.

**Example 2.** We consider now two obstacles. The cross-shaped obstacle is retained but slightly translated and is now set as follows:

$$ \begin{pmatrix} x_1(\theta) \\ x_2(\theta) \end{pmatrix} = \left[ 0.25(1 + 0.4 \cos(4\theta)) \begin{pmatrix} \cos(\theta) \\ \sin(\theta) \end{pmatrix} + 0.35 \begin{pmatrix} 1 \\ 1 \end{pmatrix} \right], \quad (\theta \in [0, 2\pi]). $$

We add an ellipse-shaped obstacle parameterized as:

$$ \begin{pmatrix} x_1(\theta) \\ x_2(\theta) \end{pmatrix} = \left[ \begin{pmatrix} 0.25 \cos(\theta) \\ 0.1 \sin(\theta) \end{pmatrix} - 0.45 \begin{pmatrix} 1 \\ 1 \end{pmatrix} \right], \quad (\theta \in [0, 2\pi]). $$

In this case, we can increase the parameter $m$ to $m = 19$ (corresponding to 37 measurements). Some level lines of the function $\Theta_{17}$ (corresponding to values of order $1/17$) are displayed in Fig. 5. Obstacle positions and sizes are recovered
<figure><img src="image_5.png" /><figcaption>Figure 3: Reconstruction of a single obstacle with disks obtained by solving the Prony's system corresponding to the harmonic moments for $n = 1, 2, 3$ (first picture), $n = 4$ (second picture) and $n = 5$ (third picture).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: Reconstruction of the cross-shaped obstacle after applying the partial balayage operator to the case $n = 5$.</figcaption></figure>

satisfactorily. Reconstruction quality is better for the ellipse than for the crosse (convexity seems to play a role for this method).

Using the first 18 harmonic moments only, we apply now the method of Section 7 and plot the corresponding disks (as explained for the first example) for $n = 2, 5$ and $8$ in Fig 6. For $n = 2$ and $n = 5$, no partial balayage step is needed as the disks do not overlap. For $n = 8$, we apply the partial balayage operator by solving the convex minimization problem (29), which provides the quadrature domains represented in Fig. 7 (the two small spurious disks are neglected). Two branches of the cross can be identified and the shape of the ellipse is correctly recovered. Once again, this is a very satisfactory reconstruction for this demanding example.

*Example 3.* The last case we consider is borrowed from [5]. The obstacles are two squares centered at $(-0.6, 0.3)$ and $(0.6, 0.3)$ with a distance between the center and the vertices equal to 0.2. The interest of studying this case is twofold: the obstacle boundaries are not smooth (unlike examples 1 and 2), and results can be compared with those provided in [5] where the obstacles are a priori assumed to be star-shaped. The parameter $m$ is set to $m = 18$ (35 measurements). The results obtained with the method of Section 6 are displayed in Fig. 8 and those obtained with the method of Section 7 in Fig. 9 (no partial balayage step is applied). Both methods deliver highly satisfactory reconstructions.

# 9 Conclusion

By transforming the classical inverse problem of reconstructing immersed obstacles using boundary measurements into a shape-from-moment problem, we have been able to develop algorithms leading to remarkably accurate reconstructions in this context. The methods we propose are not iterative and require no a priori information about the obstacles. A number of legitimate points remain to be addressed in a future paper:

1. Stability issues. The use of noisy data has not been rigorously tested in this work, but some simple principles can nevertheless be stated. The reconstruction methods are based on a certain number of measurements. *The higher the number, the more accurate the reconstruction, but also the more sensitive to noise.* For example, the reconstruction of Fig. 6 with only two disks can be achieved with 9 measurements, and the reconstruction is quite stable. However, this needs to be precisely quantified.
<figure><img src="image_7.png" /><figcaption>Figure 5: Reconstruction of two obstacles using the method described in Section 6. Some level lines of the function $\Theta_{17}$ are represented.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: Reconstruction with 2, 5 and 8 disks. Two small spurious disks are observed in the case $n = 8$. They are neglected for the partial balayage step.</figcaption></figure>

2. Incomplete data. One may want to take measurements only along part of the outer boundary $\Gamma_0$ and Theorem 3.1 states that this is theoretically sufficient to identify the obstacles. This additional difficulty has not been taken into account in this work and deserves to be studied further.

## A The Dirichlet-to-Neumann operator

Let $\Omega$ be a $C^{0,1}$ bounded domain and denote by $\Gamma$ its boundary. For every $p \in H^{1/2}(\Gamma)$, let $u_p$ be the unique function achieving:

$$ \min \{ \|\nabla \theta\|_{L^2(\Omega)} : \theta \in H^1(\Omega), \gamma_{\Gamma}^d \theta = p \}. $$

The Dirichlet-to-Neumann operator is the operator $D_{\Omega}: H^{1/2}(\Gamma) \to H^{-1/2}(\Gamma)$, $D_{\Omega}p = \gamma_{\Gamma}^n u_p$. Since, for every $p, q \in H^{1/2}(\Gamma)$:

$$ \langle D_{\Omega}p, q \rangle = -(\nabla u_p, \nabla u_q)_{L^2(\Omega)}, $$

the operator $D_{\Omega}$ is self-adjoint.

If $\Omega$ is $C^{1,1}$, by elliptic regularity, $D_{\Omega}$ maps continuously $H^{3/2}(\Omega)$ into $H^{1/2}(\Omega)$. We deduce:

**Proposition A.1.** The operator $D_{\Omega}$ extends by density as a self-adjoint operator from $H^{-1/2}(\Gamma)$ to $H^{-3/2}(\Gamma)$, i.e. for every $p \in H^{-1/2}(\Gamma)$ and every $q \in H^{3/2}(\Gamma)$ (with obvious notations):

$$ \langle D_{\Omega}p, q \rangle_{-\frac{3}{2}, \frac{3}{2}} = \langle p, D_{\Omega}q \rangle_{-\frac{1}{2}, \frac{1}{2}}. $$
<figure><img src="image_9.png" /><figcaption>Figure 7: Quadrature domains obtained by applying the partial balayage operator to the case $n = 8$ of Fig. 6.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8: Reconstruction of two squares using the method described in Section 6. Some level lines of the function $\Theta_{16}$ are represented. The picture on the right is a close-up of the right square and its reconstruction.</figcaption></figure>

## B The harmonic Dirichlet problem in $L^2$

Let $\Omega$ be a $C^{1,1}$ bounded domain, denote by $\Gamma$ its boundary and define $V = H^2(\Omega) \cap H_0^1(\Omega)$. For any $q \in H^{1/2}(\Gamma)$ define $u_q$ as the unique function in $V$ achieving:

$$ \min \{ \|\Delta\theta\|_{L^2(\Omega)} : \gamma_{\Gamma}^{n}\theta = q \}. $$

Then define the operator $T: H^{1/2}(\Gamma) \to H^{-1/2}(\Gamma)$, $Tq = \gamma_{\Gamma}^{d}\Delta u_q$. For every $q, p \in H^{1/2}(\Gamma)$, we have:

$$ \langle Tq, p \rangle_{-\frac{1}{2}, \frac{1}{2}} = -(\Delta u_q, \Delta u_p)_{L^2(\Omega)}, $$

and therefore $T$ is an isomorphism from $H^{1/2}(\Gamma)$ onto $H^{-1/2}(\Gamma)$.

**Proposition B.1.** For every $p \in H^{-1/2}(\Gamma)$, there exists a unique function $v$ in $L^2(\Omega)$, harmonic in $\Omega$ and such that $\gamma_{\Gamma}^{d}v = p$.

*Proof.* The function $v = \Delta u_{T^{-1}p}$ satisfies the conditions. To prove uniqueness, assume that $v$ is in $L^2(\Omega)$, harmonic in $\Omega$ and such that $\gamma_{\Gamma}^{d}v = 0$. Let $u$ be the unique solution in $H_0^1(\Omega)$ to the problem $\Delta u = v$. By elliptic regularity, $u$ is in $V$ and, integrating by parts, $\|\Delta u\|_{L^2(\Omega)}^2 = 0$. The conclusion follows. $\square$

Note that in a $C^{0,1}$ domain, uniqueness is not guaranteed (an example of square integrable non-zero harmonic function with vanishing Dirichlet boundary conditions in a $C^{0,1}$ bounded domain is provided in [6]).

## References

[1] C. Alvarez, C. Conca, L. Friz, O. Kavian, and J. H. Ortega. Identification of immersed obstacles via boundary measurements. *Inverse Problems*, 21(5):1531–1552, 2005.
<figure><img src="image_11.png" /><figcaption>Figure 9: Reconstruction of two squares with 2 and 8 disks. Small spurious disks are observed in the case $n = 8$.</figcaption></figure>

[2] C. J. S. Alves, R. Kress, and A. L. Silvestre. Integral equations for an inverse boundary value problem for the two-dimensional Stokes equations. *J. Inverse Ill-Posed Probl.*, 15(5):461–481, 2007.

[3] S. Axler, P. Bourdon, and W. Ramey. *Harmonic function theory*, volume 137 of *Grad. Texts Math*. New York, NY: Springer, 2nd ed. edition, 2001.

[4] M. Badra, F. Caubet, and M. Dambrine. Detecting an obstacle immersed in a fluid by shape optimization methods. *Math. Models Methods Appl. Sci.*, 21(10):2069–2101, 2011.

[5] F. Caubet, C. Conca, and M. Godoy. On the detection of several obstacles in 2d Stokes flow: topological sensitivity and combination with shape derivatives. *Inverse Probl. Imaging*, 10(2):327–367, 2016.

[6] M. Dauge, M.-S. Lubuma, and S. Nicaise. Coefficients des singularités pour le problème de Dirichlet sur un polygone. (Coefficients of singularities for the Dirichlet problem in a polygon). *C. R. Acad. Sci., Paris, Sér. I*, 304:483–486, 1987.

[7] A. Gerber-Roth, A. Munnier, and K. Ramdani. A reconstruction method for the inverse gravimetric problem. To appear in *SMAI J. of Comp. Math*. https://inria.hal.science/hal-04157036/

[8] V. Girault and P.-A. Raviart. *Finite element methods for Navier-Stokes equations*, volume 5 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, 1986. Theory and algorithms.

[9] G. H. Golub, P. Milanfar, and J. Varah. A stable numerical method for inverting shape from moments. *SIAM J. Sci. Comput.*, 21(4):1222–1243, 2000.

[10] B. Gustafsson. On quadrature domains and an inverse problem in potential theory. *J. Anal. Math.*, 55:172–216, 1990.

[11] B. Gustafsson, M. Putinar, E. B. Saff, and N. Stylianopoulos. Bergman polynomials on an archipelago: estimates, zeros and shape reconstruction. *Adv. Math.*, 222(4):1405–1460, 2009.

[12] B. Gustafsson and M. Sakai. Properties of some balayage operators, with applications to quadrature domains and moving boundary problems. *Nonlinear Anal., Theory Methods Appl.*, 22(10):1221–1245, 1994.

[13] B. Gustafsson and H. Shahgholian. Existence and geometric properties of solutions of a free boundary problem in potential theory. *J. Reine Angew. Math.*, 473:137–179, 1996.

[14] W. E. Langlois and M. O. Deville. *Slow viscous flow*. Cham: Springer, 2nd ed. edition, 2014.

[15] A. Munnier and K. Ramdani. Conformal mapping for cavity inverse problem: an explicit reconstruction formula. *Applicable Analysis*, 96(1):108–129, 2017.

[16] A. Munnier and K. Ramdani. Calderón cavities inverse problem as a shape-from-moments problem. *Quarterly of Applied Mathematics*, 2018.

[17] A. Munnier. Invertibility criteria for the biharmonic single-layer potential. Preprint 2023. https://arxiv.org/abs/2307.03160