# A reconstruction method for the inverse gravimetric problem

Anthony Gerber-Roth, Alexandre Munnier, Karim Ramdani

► To cite this version:

Anthony Gerber-Roth, Alexandre Munnier, Karim Ramdani. A reconstruction method for the inverse gravimetric problem. SMAI Journal of Computational Mathematics, 2023, 9, pp.197-225. 10.5802/smai-jcm.99. hal-04157036v2

HAL Id: hal-04157036

https://inria.hal.science/hal-04157036v2

Submitted on 9 Oct 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A reconstruction method for the
inverse gravimetric problem

ANTHONY GERBER-ROTH, ALEXANDRE MUNNIER & KARIM RAMDANI

Volume 9 (2023), p. 197-225.

https://doi.org/10.5802/smai-jcm.99

© The authors, 2023.
# A reconstruction method for the inverse gravimetric problem

ANTHONY GERBER-ROTH<fnref n="1" />

ALEXANDRE MUNNIER<fnref n="2" />

KARIM RAMDANI<fnref n="3" />

**Abstract.** We propose a reconstruction method to solve the inverse gravimetric problem with constant mass density. The method is based on the computation of the harmonic moments of the unknown domain. Convergence results are proved and numerical experiments are provided to illustrate the method and show its efficiency.

**2020 Mathematics Subject Classification.** 31A25, 86A22, 35R30.

**Keywords.** inverse gravimetric problem, shape from moments problem, Prony's system, quadrature domains, partial balayage of measure.

## 1. Introduction

### 1.1. Problem setting

Let $\omega$ be an open and bounded set in $\mathbb{R}^2$ with Lipschitz boundary. We assume that $\bar{\omega}$ is included in the open unit disk $B$ and we denote by $\Gamma$ the boundary of $B$ (see Figure 1.1).

<figure><img src="image_2.png" /><figcaption>FIGURE 1.1. Geometry of the problem</figcaption></figure>

The gravitational potential generated by a uniform mass distribution in $\omega$ is given by the Newtonian potential

$$U_{\omega}(x) = (G * 1_{\omega})(x) = \int_{\omega} G(x - y) \, dm(y), \quad \text{for all } x \in \mathbb{R}^2,$$

<footnote n="1">¹ Université de Lorraine, CNRS, Inria, IECL, F-54000 Nancy, France
*E-mail address:* anthony.gerber-roth@univ-lorraine.fr</footnote>

<footnote n="2">² Université de Lorraine, CNRS, Inria, IECL, F-54000 Nancy, France
*E-mail address:* alexandre.munnier@univ-lorraine.fr</footnote>

<footnote n="3">³ Université de Lorraine, CNRS, Inria, IECL, F-54000 Nancy, France
*E-mail address:* karim.ramdani@inria.fr.</footnote>

<footnote>https://doi.org/10.5802/smai-jcm.99</footnote>

<footnote>© The authors, 2023</footnote>
where $G(x) := -1/(2\pi) \ln|x|$ is the fundamental solution of the Laplacian in $\mathbb{R}^2$ and $m$ is the Lebesgue measure. We can easily verify that the function $U_\omega$ satisfies the Laplace equation

$$ -\Delta U_{\omega} = \mathbf{1}_{\omega} \quad \text{in } \mathbb{R}^{2}, \qquad (1.1) $$

and admits the following behavior at infinity ($|\omega|$ stands for the area of $\omega$):

$$ U_{\omega}(x) = -\frac{1}{2\pi}|\omega|\ln|x| + O(1/|x|). \qquad (1.2) $$

The inverse problem we are interested in is often referred to as the inverse gravimetric problem and reads as follows: *how to reconstruct $\omega$ from the knowledge of $\nabla U_\omega$ on $\Gamma$?*

**Remark 1.1.** Two comments are in order. First, it is natural to assume that the available measurement is $\nabla U_\omega$ (and not $U_\omega$) on $\Gamma$ since this quantity represents the gravitational field, the actual quantity measured in real experiments, typically in geophysics. Second, it is possible to recover $U_\omega$ on $\Gamma$ from the knowledge of $\nabla U_\omega$. Indeed, integrating the identity (1.1) over $B$ and applying Green's formula leads to the expression of the area of $\omega$:

$$ |\omega| = - \int_{\Gamma} \partial_{n} U_{\omega} d\sigma. $$

Then, the expression of $U_\omega$ in $B^+ = \mathbb{R}^2 \setminus \bar{B}$ follows by remarking that $U_\omega = u + |\omega|G$ in $B^+$ where the function $u$ is the unique solution of the following exterior Neumann problem (see [29, Theorem 8.18] for the well-posedness):

$$ \left\{ \begin{array}{ll} -\Delta u = 0 & \text{in } B^+ \\ \partial_n u = \partial_n U_\omega - |\omega| \partial_n G & \text{on } \Gamma \\ u = O(1/|x|) & \text{as } |x| \to +\infty. \end{array} \right. $$

## 1.2. Background on uniqueness and stability

The inverse problem stated above is known to be severely ill-posed as proved in [20, Section 3.4], since even uniqueness is not guaranteed in general. Nonetheless, uniqueness results can be established under either one of the following assumptions:

(S) $\omega$ is star-shaped with respect to its center of gravity.

(C) $\omega$ is convex in one direction (say $x_1$), i.e the intersection of any straight line parallel to the $x_1$-axis with $\omega$ is an interval.

More precisely, it holds [20, Theorem 2.2.1]:

**Theorem 1.2.** Assume that $\omega_1$ and $\omega_2$ are two domains both satisfying either (S) or (C). If $\nabla U_{\omega_1} = \nabla U_{\omega_2}$ on $\Gamma$, then $\omega_1 = \omega_2$.

It is worth mentioning that these uniqueness results have been improved by Gardiner and Sjödin [9, 10] that proved that for solid domains, uniqueness is ensured if we know that one of the two domains is convex.

A stability result can also be found in the literature, at least in the specific class of star-shaped domains having the same center of gravity. This common center of gravity is assumed to be 0 for the rest of this subsection for the sake of simplicity. In polar coordinates, any such a domain $\omega$ can be described by means of a function $\rho_\omega: [0, 2\pi[ \to \mathbb{R}^+$ defined by:

$$ \rho_{\omega}(\theta) = \max\{\rho \ge 0 : (\rho \cos \theta, \rho \sin \theta) \in \omega\} \quad \text{for all } \theta \in [0, 2\pi[. $$

Using this notation we can state the following logarithmic-type stability estimate [20, Theorem 3.6.1]:
**Theorem 1.3.** Let $M, \rho_-, \rho_+$ and $\alpha$ be four positive constants and let:

$$ \mathcal{R} = \{ \rho \in C^{0,2+\alpha}([0, 2\pi[) : \rho_- < \rho < \rho_+, \quad \|\rho\|_{C^{0,2+\alpha}} \le M \} $$

Then, there exists $\kappa > 0$ such that for all star-shaped domains $\omega_1$ and $\omega_2$ satisfying $\rho_{\omega_1}, \rho_{\omega_2} \in \mathcal{R}$ we have

$$ \|\rho_{\omega_1} - \rho_{\omega_2}\|_{L^{\infty}} \le \kappa |\ln \varepsilon|^{-1/\kappa}, \quad (1.3) $$

provided $\|\nabla U_{\omega_1} - \nabla U_{\omega_2}\|_{L^{\infty}(\Gamma)} \le \varepsilon$.

### 1.3. Overview of the method

The aim of this paper is to develop an algorithm that generates a sequence of domains $(\omega_N)_N$ converging (in some sense) to the unknown domain $\omega$, or at least whose gravitational fields fit the targeted one (i.e, is graviequivalent to $\omega$). To achieve this goal, we only make use of the harmonic moments of $\omega$, quantities that can easily be deduced from the measurements. Indeed, for every function $v$ harmonic in $B$, Green's formula yields:

$$ \int_{\omega} v \, dm = \int_{\Gamma} (\partial_n v U_{\omega} - \partial_n U_{\omega} v) \, d\sigma, \quad (1.4) $$

where we recall that $\nabla U_{\omega}$ and $U_{\omega}$ are both supposed to be known on $\Gamma$ (see Remark 1.1). Thus, we are led to study a general shape-from-moments problem which can be stated as follows: *how to reconstruct a domain $\omega$ from the knowledge of its harmonic moments?*

Throughout the paper, any complex number $z = x_1 + ix_2$ will be identified with the point $(x_1, x_2)$ of the plane. Let us denote by $z^\ell$ ($\ell \ge 0$) the harmonic monomials in the plane. Given an integer $N \ge 1$, the main idea of this work is to construct a sequence of domains $(\omega_N)_N$ such that:

$$ \int_{\omega} z^{\ell} \, dm(z) = \int_{\omega_N} z^{\ell} \, dm(z), \quad \text{for all } \ell = 0, \dots, 2N - 1. $$

To achieve this goal we proceed in two steps. For any positive integer $N$:

Step 1 Construction of a quadrature formula associated with the data, i.e. computation of weights $c_1, \dots, c_N \in \mathbb{C}$ and nodes $z_1, \dots, z_N \in \mathbb{C}$ such that:

$$ \int_{\omega} z^{\ell} \, dm(z) = \sum_{k=1}^{N} c_k z_k^{\ell}, \quad \text{for all } \ell = 0, \dots, 2N - 1. \quad (1.5) $$

The values of the $c_j$ and $z_j$ ($j = 1, \dots, N$) will be obtained by means of the so-called Prony's method.

Step 2 Construction of a quadrature domain associated with the previous quadrature formula, i.e. determination of a domain $\omega_N$ (if any) satisfying identity (1.5) but for any harmonic polynomial this time:

$$ \int_{\omega_N} z^{\ell} \, dm(z) = \sum_{k=1}^{N} c_k z_k^{\ell}, \quad \text{for all } \ell \ge 0. $$

Such a domain is called a “quadrature domain” and for given $c_k$ and $z_k$ ($k = 1, \dots, N$), its existence is not guaranteed in general. We shall test two approaches to reconstruct this domain: the first one rests on purely complex analysis tools while the second one, called “partial balayage”, is formulated as a convex minimization problem.
Notice that our method may be applied to any shape-from-moments problem, though we are able to obtain convergence results only for the gravimetric problem.

Our method is based on an original combination of classical notions. The relation between gravitational potential and quadrature domains is natural (see for example [12]). On the other hand, although Prony's method (and related approaches) is a classical tool in signal processing (particularly in exponential fitting), its use for the reconstruction of quadrature domains is not very common.

## 1.4. Related works

Inverse gravimetric problems may play an important role in future practical applications (primarily in geophysics), more especially since recent advances in gravitational field measurements [39]. As far as theoretical questions are concerned, the book by Isakov [20] remains an essential reference. Let us also mention the monograph of Cherednichenko [6]. Following Potthast in his survey paper [32], it is relevant to classify the methods for solving geometric inverse problems in two main categories: direct and iterative ones.

Our strategy belongs to the first category, and has similarities with those presented in [11, 13]. In these articles, the authors address the problem of reconstructing an unknown domain from its moments. Only polygonal shapes are considered in the former while all the moments (and not only the harmonic ones) are required for the reconstruction in the latter.

Moments methods were also used in [33] to recover the Fourier expansion of the boundaries' parametrizations of domains having the property (S). Finally, let us also mention [5], in which a ball fitting approach, evoking our ball representation of quadrature domains, is used.

In the class of iterative methods, a first glance shows the predominance of those based on level-set optimization [21, 22, 26, 27, 28, 40]. Before the spreading of this powerful general tool, Hettlich and Rundell [18] tackled the problem using shape derivatives coupled to a Newton-type method, while Zidarov implemented a least squares approximation, as described in his classic book [41]. Finally, Kress and Rundell [25] proposed an original rewriting of this problem using boundary integral equations.

## 1.5. Outline of the paper

The paper is organized as follows: Step 1 and Step 2 are elaborated in Section 2 and 3 respectively. Next, we establish convergence results in Section 4 and finally we present some numerical experiments and give some details about their implementation in Section 5.

# 2. Construction of the quadrature formula

Recall that Step 1 of our method consists in seeking, for any positive integer $N$, weights $c_1, \dots, c_N \in \mathbb{C}$ and nodes $z_1, \dots, z_N \in \mathbb{C}$ such that:

$$ \sum_{k=1}^{N} c_k z_k^{\ell} = \tau_{\ell}, \quad \text{for all } \ell = 0, \dots, 2N-1, \qquad (2.1) $$

where the quantities in the right hand side are the harmonic moments of the domain $\omega$ defined by:

$$ \tau_{\ell} = \int_{\omega} z^{\ell} dm(z), \quad \text{for all } \ell = 0, \dots, 2N-1. \qquad (2.2) $$

The non-linear system (2.1) turns out to be a so-called Prony's system. Prony's systems have been extensively studied, as they appear in several fields of theoretical and applied mathematics, such as Padé's approximants, signal processing or error correction codes.
**Remark 2.1.** Observe that (2.1) can be rewritten as:

$$ \int z^l d\mu_N = \int z^l \mathbb{1}_\omega(z) dm(z), \quad \text{for all } l = 0, \dots, 2N-1, \qquad (2.3) $$

where the measure $\mu_N$ is equal to the linear combination of Dirac masses $\sum_{k=1}^N c_k \delta_{z_k}$. Intuitively, identity (2.3) suggests that we are trying to approximate the measure $\mathbb{1}_\omega$ with point masses. It is thus natural to expect the weights $c_k$ to be positive for all $k$. However we do not add this physical constraint, for the resulting constrained system may have no solution.

We formulate the problem of solving (2.1) as follows:

$$ \left\{ \begin{array}{l} \text{Let } N \text{ be a positive integer and } \tau_0, \dots, \tau_{2N-1} \in \mathbb{C} \text{ be given. Determine for some} \\ \text{positive integer } p, z_1, \dots, z_p \in \mathbb{C} \text{ pairwise distincts and } c_1, \dots, c_p \in \mathbb{C} \setminus \{0\} \text{ such} \\ \text{that} \\ \displaystyle \sum_{k=1}^p c_k z_k^l = \tau_l, \quad \text{for all } l = 0, \dots, 2N-1. \end{array} \right. \qquad (2.4) $$

Let us emphasize that in this statement, the $\tau_l$ are arbitrary complex numbers and do not need to be harmonic moments as in (2.2).

## 2.1. Existence and uniqueness

When $p$ is allowed to be as large as $2N$, Problem (2.4) admits always a solution, as proved for instance in the first theorem in [31, Section 2]. However, this does not meet our requirement in (2.1), which is $p \le N$. Moreover, the nodes $z_k$ of this solution lie on the unit circle, which is not satisfactory for our purpose. Indeed, we want to reconstruct the unknown domain $\omega$ (compactly contained in the unit disk by assumption), so we expect the nodes to be located in $\omega$ (this intuition will become obvious in Section 3).

The next results stated in this subsection might have been already proved, but we have not been able to find them in the literature. The proofs are given in Appendix A.

We begin with introducing the matrix:

$$ \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}. \qquad (2.5a) $$

This Hankel structure plays a central role in the analysis, as illustrated by the following uniqueness result:

**Proposition 2.2.** Assume that Problem (2.4) admits a solution with $p \le N$. Then, if $\mathbb{H}_0^{(N)}$ is invertible, this solution is the unique solution for which $p \le N$ and necessarily $p = N$.

From this proposition, we deduce:

**Corollary 2.3.** Problem (2.4) admits at most one solution such that $p = N$.

The statement of the existence result requires introducing the following 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 & \ddots & \vdots & \vdots \\ \tau_{N-1} & \tau_N & \cdots & \tau_{2N-2} & \tau_{2N-1} \\ 1 & z & \cdots & z^{N-1} & z^N \end{vmatrix}. $$
**Theorem 2.4.** *Problem (2.4) admits a solution for which $p = N$ 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$.*

## 2.2. Numerical solutions

A review of available methods to solve (2.1) can be found in [19]. Roughly speaking (see, e.g., [38, Section 4]), one can distinguish direct nonlinear minimization methods (typically nonlinear least squares methods), recurrence-based methods (such as the original Prony's method), and subspace methods (such as Pisarenko's method, MUSIC, ESPRIT and matrix pencil methods). In this work, we apply the matrix pencil method, as described in [11].

Let us introduce another Hankel matrix of size $N$ (obtained by applying a shift to the indices of $\mathbb{H}_0^{(N)}$):

$$ \mathbb{H}_1^{(N)} = \begin{pmatrix} \tau_1 & \tau_2 & \dots & \tau_N \\ \tau_2 & \tau_3 & \dots & \tau_{N+1} \\ \vdots & \vdots & \ddots & \vdots \\ \tau_N & \tau_{N+1} & \dots & \tau_{2N-1} \end{pmatrix}. \qquad (2.5b) $$

First we denote by $z_1, \dots, z_N$ the (complex) eigenvalues that verify 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. $$

Then, if the linear system:

$$ \begin{pmatrix} 1 & 1 & \dots & 1 \\ z_1 & z_2 & \dots & z_N \\ z_1^2 & z_2^2 & \dots & z_N^2 \\ \vdots & \vdots & \ddots & \vdots \\ z_1^{N-1} & z_2^{N-1} & \dots & 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}, $$

admits solutions $c_1, \dots, c_N$, one can verify that $z_1, \dots, z_N$ and $c_1, \dots, c_N$ solve (2.1) (see [11, Section 3]).

Since the inverse gravimetric problem is ill-posed, it is not surprising that it leads to Prony's systems which are generally ill-conditioned as well (see [3]). Actually, the above numerical method involves solving two ill-conditioned problems: a generalized eigenvalue problem with Hankel matrices [35] and a Vandermonde system. To improve the overall condition number, we use in our implementation the regularization techniques detailed in [11, Section 4]. It is worth mentioning also that numerical instabilities increase along with $N$.

## 2.3. An example of Prony's system without solution

In this section, we provide an example of domain whose moments lead to a Prony's system that does not have a solution (and hence for which the first step of our reconstruction algorithm fails). Consider the domain $\omega$ whose boundary is described in polar coordinates by $\rho(\theta) = 1 + 2a \cos(\theta)$ with $0 < a < 1/2$ (see Figure 2.1).

For every $N \ge 1$, let the Hankel matrix $\mathbb{H}_0^{(N)}$ and the polynomial $P_N$ be defined as in Section 2.1, where the measurements $(\tau_\ell)_{\ell \ge 0}$ are given by (2.2).

**Proposition 2.5.** *The following assertions hold true.*

(1) $P_2$ has a double root.

(2) $\text{rank} \mathbb{H}_0^{(N)} = 2$, for all $N \ge 2$.
<figure><img src="image_3.png" /><figcaption>FIGURE 2.1. A “pathological” domain with $a = 0.35$.</figcaption></figure>

**Proof.** For all $\ell \ge 0$, direct computations lead to:

$$ \tau_{\ell} = \frac{1}{\ell + 2} \int_{0}^{2\pi} (1 + a(e^{i\theta} + e^{-i\theta}))^{\ell+2} e^{i\ell\theta} d\theta, \qquad (2.6) $$

and since:

$$ (1 + a(e^{i\theta} + e^{-i\theta}))^{\ell+2} = \sum_{j=0}^{\ell+2} \sum_{k=0}^{j} \binom{\ell+2}{j} \binom{j}{k} a^j e^{-i(2k-j)\theta}, $$

we end up with:

$$ \tau_{\ell} = \pi a^{\ell} (\ell + 1 + 2a^2). $$

Consequently, $\det \mathbb{H}_0^{(2)} = -\pi^2 a^2$ and

$$ P_2(z) = \pi^2 \begin{vmatrix} 1 + 2a^2 & a(2+2a^2) & a^2(3+2a^2) \\ a(2+2a^2) & a^2(3+2a^2) & a^3(4+2a^2) \\ 1 & z & z^2 \end{vmatrix} = -\pi^2 a^2 (z-a)^2, $$

which proves the first point. Concerning the second point, we have:

$$ \mathbb{H}_0^{(N)} = \pi \begin{pmatrix} 1 + 2a^2 & a(2+2a^2) & \cdots & a^{N-1}(N+2a^2) \\ a(2+2a^2) & a^2(3+2a^2) & \cdots & a^N(N+1+2a^2) \\ \vdots & \vdots & \ddots & \vdots \\ a^{N-1}(N+2a^2) & a^N(N+1+2a^2) & \cdots & a^{2N-2}(2N-1+2a^2) \end{pmatrix}. $$

Multiplying the k-th row by $a^{-(k-1)}$ and the j-th column by $a^{-(j-1)}$, we obtain that this matrix has the same rank as

$$ \begin{pmatrix} 1 + 2a^2 & (2 + 2a^2) & \cdots & (N + 2a^2) \\ (2 + 2a^2) & (3 + 2a^2) & \cdots & (N + 1 + 2a^2) \\ \vdots & \vdots & \ddots & \vdots \\ (N + 2a^2) & (N + 1 + 2a^2) & \cdots & (2N - 1 + 2a^2) \end{pmatrix}, $$

which has the same rank as

$$ \begin{pmatrix} 1 + 2a^2 & 1 & \cdots & 1 \\ (2 + 2a^2) & 1 & \cdots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ (N + 2a^2) & 1 & \cdots & 1 \end{pmatrix}, $$

as it can be seen by replacing column $C_j$ by $C_j - C_{j-1}$ for all $j \ge 2$.

<figure><img src="image_4.png" /></figure>
We are now in position to prove:

**Proposition 2.6.** For every $N \ge 2$, Problem (2.4) has no solution for which $p \le N$.

**Proof.** The result is proved by induction on the index $N$. For $N = 2$, $\mathbb{H}_0^{(2)}$ is invertible, so it follows from Proposition 2.2 that every solution for which $p \le 2$ would actually be such that $p = 2$. However, since $P_2$ has a double root, Theorem 2.4 asserts that Problem (2.4) has no such solution. Assume now that the claim of the Proposition holds true for some $N' \ge 2$. We first notice that the coefficient in front of the monomial $z^{N'+1}$ in $P_{N'+1}$ is

$$(-1)^{N'+1} \det \mathbb{H}_0^{(N'+1)} = 0,$$

since $\mathrm{rank} \mathbb{H}_0^{(N'+1)} = 2$. Consequently, $P_{N'+1}$ cannot have $N' + 1$ simple roots and it follows that Problem (2.4) for $N = N' + 1$ has no solution for which $p = N' + 1$. If a solution with $p < N' + 1$ existed, it would be also a solution to Problem (2.4) for $N = N'$, yielding a contradiction with the assumption. ■

Pathological domains, like the example provided in this subsection, seem to be interesting mainly from a theoretical point of view. In practice, we have never encountered this type of situation. However, the cancellation of $\det \mathbb{H}_0^{(N)}$ for some values of $N$ (at least in the floating-point precision sense) is common and arises typically when $N$ is too large (see also Remark 4.8 for a related discussion).

### 3. Construction of the quadrature domain

We address now the second step of our method summarized in Section 1.3. This step consists in constructing a domain associated with a given quadrature formula. More precisely, given $N$ nodes $z_1, \dots, z_N \in \mathbb{C}$ and $N$ weights $c_1, \dots, c_N \in \mathbb{C}$ (those computed in the first step of the reconstruction algorithm detailed in Section 2), we want to construct a domain $\omega_N$ such that for all harmonic functions $v$ on $\omega_N$:

$$\int_{\omega_N} v \, dm = \sum_{k=1}^{N} c_k v(z_k). \qquad (3.1)$$

Obviously, such a domain will satisfy in particular:

$$\int_{\omega_N} z^\ell dm(z) = \int_{\omega} z^\ell dm(z) = \tau_\ell, \quad \text{for all } \ell = 0, \dots, 2N-1. \qquad (3.2)$$

The domains satisfying quadrature identities in the genus of (3.1) are called quadrature domains. The simplest quadrature domain is provided by the disk $D(a, R)$ of center $a \in \mathbb{C}$ and radius $R > 0$ for which, according to the mean value property:

$$\int_{D(a,R)} v \, dm = \pi R^2 v(a),$$

for every function $v$ harmonic in $D(a, R)$. Let us now collect some of the main properties of quadrature domains that will be needed in the sequel (for a complete introduction to quadrature domains, we refer the interested reader to [16] and references therein).
## 3.1. Definitions and properties

**Definition 3.1.** A planar domain $\Omega$ is called a quadrature domain associated to a complex Radon measure $\mu$ when:

$$ \int_{\Omega} v \, dm = \mu(v), \qquad (3.3) $$

for all complex analytic integrable function $v$ in $\Omega$. The set of all domains satisfying this property is denoted by $Q(\mu, AL^1)$. Replacing “complex analytic” by “harmonic” in this definition leads to the definition of harmonic quadrature domains which form the set $Q(\mu, HL^1)$.

Eventually, a domain $\Omega$ satisfying:

$$ \int_{\Omega} v \, dm \geq \mu(v), \qquad (3.4) $$

for all subharmonic integrable function $v$, is called a subharmonic quadrature domain. The set of all the subharmonic quadrature domains is denoted by $Q(\mu, SL^1)$.

**Remark 3.2.** The following inclusions hold:

$$ Q(\mu, SL^1) \subset Q(\mu, HL^1) \subset Q(\mu, AL^1), \qquad (3.5) $$

(the latter is obvious and for the former, notice that when $v$ is a harmonic function, $v$ and $-v$ are both subharmonic).

Returning to our main concern, that is to construct domains satisfying (3.1), we shall focus in the sequel on the particular case where $\mu$ is an atomic measure. More precisely, we assume from now on that:

$$ \mu = \mu_N := \sum_{k=1}^{N} c_k \delta_{z_k}, \qquad (3.6) $$

where $N$ is a positive integer, $z_1, \dots, z_N$ are pairwise distinct complex numbers (the nodes), and $c_1, \dots, c_N$ are complex numbers (the weights).

**Remark 3.3.** Depending on the type of quadrature domain considered, the weights $c_k$ must satisfy additional constraints:

(1) If $\Omega \in Q(\mu_N, HL^1)$, then all the weights $c_k$ are real. Indeed, taking the imaginary part in (3.3) leads to:

$$ \sum_{k=1}^{N} \mathrm{Im}(c_k) v(z_k) = 0, $$

for all harmonic function $v$, and thus $\mathrm{Im}(c_k) = 0$.

(2) If $\Omega \in Q(\mu_N, SL^1)$, then $c_k \geq 0$ for all $k$. Indeed, let $k$ be given in $\{1, \dots, N\}$ and for every $\varepsilon > 0$ define the subharmonic function (borrowed from [36, p. 37]):

$$ u_{\varepsilon}(z) = \begin{cases} \ln|z - z_k| & \text{if } |z - z_k| \geq \varepsilon \\ \ln \varepsilon & \text{if } |z - z_k| \leq \varepsilon. \end{cases} $$

Using this function in (3.4) (with $\mu = \mu_N$) and letting $\varepsilon$ go to 0, we get the result.

As already mentioned, the disks are quadrature domains and it can be shown that they are the only ones for which $N = 1$ [34, Example 1.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 (see [4, Section 2]).
Given a set of nodes and weights, existence of a subharmonic quadrature domain is asserted by Gustafsson in [12, Theorem 2.4(vi)], providing that the weights are positive (recall that every sub-harmonic quadrature domain is also a harmonic quadrature domain by (3.5)). The construction of this domain from $\mu_N$ is strongly related to partial balayage of measures [14] and free boundary problems [15]. Things are more complicated when it comes to the question of uniqueness. Focusing again on subharmonic quadrature domains, uniqueness is proved in the already mentioned paper [12, Theorem 2.2]. In view of (3.4), uniqueness means up to a set of null Lebesgue measure and in [12, Theorem 2.2], the author claims that this domain can be chosen open. This uniqueness property is lost when we move to harmonic quadrature domains. Indeed, even in the simplest case of an atomic measure like $\mu_N$ (with positive weights), uniqueness does not hold in general. An example is provided in the paper [1], to which we refer for general questions related to uniqueness.

We shall now present two numerical methods allowing the construction of a quadrature domain from the knowledge of a given set of nodes and weights: a conformal mapping method and a convex minimization method. It is worth noticing that these methods won't always succeed if we do not assume the $c_k$ to be positive. In particular, characterizing nodes and weights for which there exists a quadrature domain is an open problem for general complex weights $c_k$.

## 3.2. Conformal mapping

Let $N$ be a positive integer, and $z_1, \dots, z_N$ be complex numbers pairwise distincts and located in the open unit disk $B$. Consider in addition complex weights $c_1, \dots, c_N$. The method used in this subsection to recover a quadrature domain from these numbers is taken from the recent article [1]. It allows for complex weights but is able to yield only quadrature domains $\omega_N$ that are connected and “solid” (i.e. an open set such that $\omega_N^e = \mathbb{R}^2 \setminus \overline{\omega_N}$ is connected and $\partial\omega_N^e = \partial\omega_N$).

More precisely, let us consider $\mu_N$ as introduced in (3.6) and assume that at least one solid domain $\omega_N$ containing the origin belongs to $Q(\mu_N, AL^1)$. Then it can be obtained as the image of the unit disk $B$ by a conformal map $\varphi_N$ taking the form:

$$ \varphi_N(z) = \sum_{k=1}^{N} \frac{z}{1 - \bar{\lambda}_k z} \quad \text{for all } z \in B, \qquad (3.7) $$

where $x_1, \dots, x_N$ and $\lambda_1, \dots, \lambda_N$ are complex parameters (with $|\lambda_k| < 1$ for every $k$). These parameters are related to the nodes $z_1, \dots, z_N$ and the weights $c_1, \dots, c_N$ through the identities:

$$ z_k = \sum_{j=1}^{N} \frac{\overline{x_j} \lambda_k}{1 - \lambda_k \overline{\lambda_j}} \quad \text{and} \quad c_k = \sum_{j=1}^{N} \frac{\overline{x_j} x_k}{(1 - \lambda_k \overline{\lambda_j})^2} \quad \text{for all } k = 1, \dots, N. \qquad (3.8) $$

To solve this system we use the pre-implemented Matlab function `fsolve`, which is based on an iterative method. Since the nodes and weights are supposed to be the harmonic moments of the domain $\omega_N$, it follows that $\sum_{k=1}^N c_k = |\omega_N|$ is real and positive. We initialize the method with the parameters:

$$ \lambda_k^0 = 0, \quad x_k^0 = \frac{r_k \sqrt{|\omega_N|}}{\sum_{j=1}^{N} |c_j|}, \quad \text{for all } k = 1, \dots, N, \qquad (3.9) $$

where $r_k := \sqrt{|c_k|/\pi}$. With these values for the parameters, the image of the unit disk by the conformal mapping $\varphi_N$ is a disk of area $|\omega_N|$ i.e. with the same first harmonic moment as $\omega_N$.

## 3.3. Convex minimization

Consider now in addition that the weights $c_1, \dots, c_N$ are real numbers and define the measure $\mu_N$ as in (3.6). We aim to construct a solid element of $Q(\mu_N, HL^1)$. For this purpose we introduce a balayage
operator that will provide, under some conditions, the characteristic function of the quadrature domain $\omega_N$. This balayage operator can be evaluated by solving a convex minimization problem. This will be done using Finite Element Method (FEM).

In the following, we denote by $M_c$ the set of real Radon measures with compact support and by $L_c^\infty$ the set of compactly supported functions in $L^\infty(\mathbb{R}^N)$. As already mentioned, we will identify two domains when they differ only by a set of measure zero.

It will be more convenient to deal with the following regularized version of the measure $\mu_N$:

$$ \tilde{\mu}_N = \sum_{k=1}^{N} \frac{c_k}{\alpha_k^2 \pi} \mathbb{1}_{B(z_k, \alpha_k)}, \qquad (3.10) $$

for given $\alpha_1, \dots, \alpha_N > 0$. The next result shows that $\tilde{\mu}_N$ and $\mu_N$ yield the same quadrature domains.

**Proposition 3.4.** We have:

$$ Q(\mu_N, HL^1) = Q(\tilde{\mu}_N, HL^1). $$

**Proof.** By the mean value property, for every harmonic function $v$ we can write

$$ \begin{aligned} \int v \, d\tilde{\mu}_N &= \sum_{k=1}^{N} \frac{c_k}{\alpha_k^2 \pi} \int_{B(z_k, \alpha_k)} v \, dm \\ &= \sum_{k=1}^{N} c_k v(z_k) \\ &= \int v \, d\mu_N, \end{aligned} $$

and the result follows.

■

We define now the balayage operator $F$:

**Definition 3.5.** For every $\mu \in M_c$, we set

$$ \mathcal{F}_{\mu} = \{u \in \mathcal{D}'(\mathbb{R}^N) : u \le U^{\mu} \text{ and } -\Delta u \le 1\}, $$

and we define:

(1) $V^\mu = \sup \mathcal{F}_\mu$,

(2) $F(\mu) = -\Delta V^\mu$,

(3) $\Omega(\mu) = \mathbb{R}^N \setminus \text{supp}(1 - F(\mu))$.

Some properties of the balayage operator $F$ are collected in the next proposition (a more complete analysis can be found in [12, 14]).

**Proposition 3.6.** Let $\mu \in M_c$, then

(1) $F(\mu) \le 1$

(2) if $\mu \in L_c^\infty$, $F(\mu) = \mathbb{1}_{\Omega(\mu)} + \mu \mathbb{1}_{\Omega(\mu)^c}$

(3) if $\mu \in L_c^\infty$ and $F(\mu) = \mathbb{1}_{\Omega}$, then $\Omega = \Omega(\mu)$ and $\Omega(\mu) \in Q(\mu, SL^1)$.
## Proof.

(1). This point follows from the definition.

(2). This is [14, Formula (2.40)].

(3). From, [14, Example 2.2], we obtain $\Omega = \Omega(\mu)$ up to a null set, and $U_{\Omega} \le U^{\mu}$ in $\mathbb{R}^2$ and $U_{\Omega} = U^{\mu}$ outside $\Omega$, which implies that $\Omega \in Q(\mu, SL^1)$ from [12, Proposition 1.2]. ■

The cornerstone of our method consists in noticing that when $F(\tilde{\mu}_N)$ is a characteristic function of a domain $\omega_N$:

$$F(\tilde{\mu}_N) = 1_{\omega_N}, \qquad (3.11)$$

then we can deduce that:

$$\omega_N \in Q(\tilde{\mu}_N, SL^1) \subset Q(\tilde{\mu}_N, HL^1) = Q(\mu_N, HL^1).$$

If in addition $\omega_N$ is solid, it is a solution to our problem. Let us emphasize that it is crucial to work with the regularized measure $\tilde{\mu}_N$ because results from [12, 14] do not apply to the singular measure $\mu_N$.

From now, in the definition (3.10) of $\tilde{\mu}_N$, we set $\alpha_k = \alpha\sqrt{|c_k|/\pi} = \alpha r_k$, with $\alpha \in (0, 1]$, for $k = 1, \dots, N$. To emphasize the dependence on the parameter $\alpha$, we will use the notation $\tilde{\mu}_N^\alpha$ for the measure $\tilde{\mu}_N$. When $\alpha = 1$, $\tilde{\mu}_N^1$ can be thought of as an overlap of weighted disks with a uniform mass distribution equal to 1.

**Proposition 3.7.** Let $\alpha \in (0, 1)$.

$$ (1) \quad \text{If } \mu_N \ge 0, \text{ then} \qquad F(\tilde{\mu}_N^\alpha) = F(\mu_N) = 1_{\Omega(\mu_N)}, \qquad (3.12) \\ \text{and } \Omega(\mu_N) \in Q(\mu_N, SL^1). $$

$$ (2) \quad \text{Conversely, if (3.12) holds, then necessarily } \mu_N \ge 0. $$

## Proof.

(1). Following [14, Example 3.2] we obtain $F(\tilde{\mu}_N^\alpha) = F(\mu_N)$, which according to Proposition 3.6, leads to:

$$F(\tilde{\mu}_N^\alpha) = 1_{\Omega(\tilde{\mu}_N^\alpha)} + \tilde{\mu}_N^\alpha 1_{\Omega(\tilde{\mu}_N^\alpha)^c} \le 1.$$

But $\tilde{\mu}_N^\alpha > 1$ on its support, and thus $\tilde{\mu}_N^\alpha 1_{\Omega(\tilde{\mu}_N^\alpha)^c} = 0$. Since $\tilde{\mu}_N^\alpha \ge 1$ on its support, it follows from [12, Corollary 2.3] that $\Omega(\tilde{\mu}_N^\alpha) \in Q(\tilde{\mu}_N^\alpha, SL^1)$. Finally, taking $v$ subharmonic, we have

$$ \begin{aligned} \int_{\Omega(\tilde{\mu}_N^\alpha)} v \, dm &\ge \int v \, d\tilde{\mu}_N^\alpha \\ &= \frac{1}{\alpha^2} \sum_{k=1}^N \int_{B(z_k, \alpha r_k)} v \, dm \\ &\ge \sum_{k=1}^N c_k v(z_k), \end{aligned} $$

by the mean value property for subharmonic functions. Hence, $\Omega(\tilde{\mu}_N^\alpha) \in Q(\mu_N, SL^1)$ which implies $\Omega(\tilde{\mu}_N^\alpha) = \Omega(\mu_N)$ up to a null set by [12, Theorem 2.2] and completes the first point of the Proposition.
(2). If (3.12) holds, then

$$U_{\Omega(\mu_N)} = U^{F(\mu_N)} = V^{\mu_N} \le U^{\mu_N},$$

which implies that:

$$\int_{\Omega(\mu_N)} G(x-y) \, dx \le \sum_{k=1}^{N} c_k G(x_k - y), \quad \text{for all } y \in \mathbb{R}^N \setminus \{x_1, \dots, x_N\}.$$

Now, take $1 \le k \le N$ and choose $y$ such that $|x_k - y| < 1$. Remarking that $G(x_k - y) \ge 0$ we can write

$$\frac{\int_{\Omega(\mu_N)} G(x-y) \, dx - \sum_{j \ne k} c_j G(x_j - y)}{G(x_k - y)} \le c_k,$$

and letting $y \to x_k$ yields $c_k \ge 0$.

Thus, when $\mu_N \ge 0$ (i.e. $c_k \ge 0$ for all $k$), computing $F(\tilde{\mu}_N^\alpha)$ directly provides the characteristic function of an element of $Q(\mu_N, HL^1)$ (and even of $Q(\mu_N, SL^1)$).

Now we focus on the practical evaluation of $F$. Keeping in mind that $F(\mu) = -\Delta V^\mu$, this can be achieved by solving a convex minimization problem to determine $V^\mu$.

**Proposition 3.8.** Let $\mu \in M_c$, and assume that $\text{supp}(\mu)$ and $\Omega(\mu)$ are compactly contained in the unit ball $B$. Then $V^\mu|_B$ is the unique minimizer of the convex functional $J$ defined by

$$J(v) = \frac{1}{2} \int_B |\nabla v|^2 \, dm - \int_B v \, dm, \quad \text{for all } v \in H^1(B), \qquad (3.13a)$$

over the following closed convex subset of $H^1(B)$:

$$K^\mu = \{v \in H^1(B) : v|_{\Gamma} = U^\mu|_{\Gamma} \text{ and } v \le U^\mu \text{ on } B\}. \qquad (3.13b)$$

**Proof.** We follow the lines of the proof of [12, Theorem 2.3]. Let us denote by $u$ the unique minimizer of (3.13), then according to (for instance) [2, Example 11.1.1], $u$ admits the following equivalent characterization: find $u \in K^\mu$ such that

$$\int_B \nabla u \cdot \nabla (v-u) \, dm \ge \int_B (v-u) \, dm, \quad \text{for all } v \in K^\mu. \qquad (3.14)$$

Taking $v = u - \varphi$ in (3.14), with $\varphi \in \mathcal{D}(B)$, $\varphi \ge 0$ leads to $-\Delta u \le 1$. Thus, recalling the definition of $K^\mu$, we see that $u \in \mathcal{F}_\mu$ and hence $u \le V^\mu$. Since $\Omega(\mu)$ is compactly contained in $B$, by [14, Section 2] we can write that

$$V^\mu = U^\mu \quad \text{on } \Gamma. \qquad (3.15)$$

Thus, $V^\mu - u$ is nonnegative and belongs to $H_0^1(B)$. As a consequence, since $-\Delta V^\mu \le 1$, we have:

$$\int_B \nabla V^\mu \cdot \nabla (V^\mu - u) \, dm = \int_B -\Delta V^\mu (V^\mu - u) \, dm \\ \le \int_B (V^\mu - u) \, dm.$$

Now, choosing $v = V^\mu$ in (3.14) shows that

$$\int_B \nabla u \cdot \nabla (V^\mu - u) \, dm \ge \int_B (V^\mu - u) \, dm.$$

Taking the difference between the two last inequalities leads to:

$$\int_B |\nabla (V^\mu - u)|^2 \, dm \le 0,$$

which implies $u = V^\mu|_B$ and concludes the proof.
We recognize in (3.13) the classical obstacle problem which is known to be well-posed [8]. Summarizing, if $\mu_N \ge 0$, we have:

$$1_{\Omega(\mu_N)} = F(\mu_N) = -\Delta V^{\tilde{\mu}_N^\alpha}, \qquad (3.16a)$$

where $V^{\tilde{\mu}_N^\alpha}|_B$ is the unique solution of the variational inequality:

$$\int_B \nabla u \cdot \nabla (v-u) dm \ge \int_B (v-u) dm, \quad \text{for all } v \in K^{\tilde{\mu}_N^\alpha}, \qquad (3.16b)$$

with:

$$K^{\tilde{\mu}_N^\alpha} = \{v \in H^1(B) : v|_{\Gamma} = U^{\tilde{\mu}_N^\alpha}|_{\Gamma} \text{ and } v \le U^{\tilde{\mu}_N^\alpha} \text{ on } B\}. \qquad (3.16c)$$

We are now in position to numerically evaluate $F$ using a FEM on $B$. We refer to [23] for a Newton-like method designed to find an approximate solution for this kind of problem. One can notice that $U^{\tilde{\mu}_N^\alpha}$ belongs to $H^1(B)$ (and hence can be discretized on a finite elements space) whereas $U^{\mu_N}$ does not, this is another argument for considering a regularization of $\mu_N$.

We implemented this method using FreeFEM++ [17]. We initialize the method with the gravitational potential $U^{\tilde{\mu}_N^1}$. This choice is motivated by the fact that, in the special case where the disks $B(z_k, r_k)$ are pairwise disjoint, $U^{\tilde{\mu}_N^1}$ turns out to be the minimizer we are looking for (indeed, in this case $U^{\tilde{\mu}_N^1} \in \mathcal{F}_{\tilde{\mu}_N^1}$).

When $\mu_N$ is not positive, $F(\tilde{\mu}_N^\alpha)$ may not be an indicator function. However, it is sometimes possible to tune $\alpha$ in such way that $\text{supp}(\tilde{\mu}_N^\alpha) \subset \Omega(\tilde{\mu}_N^\alpha)$ (see Proposition 3.6). Furthermore, adjusting this parameter can result in a solid domain, whereas in general, a domain resulting from an evaluation of $F$ may not be simply connected. Figure 3.1 illustrates the influence of the parameter $\alpha$ on $F(\tilde{\mu}_N^\alpha)$. We choose for $\mu_N$ a sum of two Dirac masses, one with positive weight and the other with negative weight. The measure $\tilde{\mu}_N^\alpha$ is therefore the weighted sum of two indicator functions of disks (with the same centers as the Dirac masses). The boundaries of these disks (and hence $\text{supp}(\tilde{\mu}_N^\alpha)$) is represented with a green dashed line. On Figure 3.1a and Figure 3.1b, $F(\tilde{\mu}_N^\alpha)$ is the indicator function of a non simply connected domain, while on Figure 3.1c, the balayage of $\tilde{\mu}_N^\alpha$ is the indicator function of a solid domain.

<figure><img src="image_5.png" /><figcaption>FIGURE 3.1. Representation of $F(\tilde{\mu}_N^\alpha)$ for different values of $\alpha$.</figcaption></figure>
## 4. Convergence results

### 4.1. Basic elements from potential theory

For reader's convenience, we recall in this subsection some basic properties from potential theory used in the sequel. We shall focus on the single layer potential and we refer to the monograph [29] for a thorough presentation of potential theory and boundary integral operators.

Let $\mathbf{D}$ be a bounded Lipschitz domain whose boundary $\mathcal{C}$ is a Jordan curve. Let $n$ be the unit normal directed towards the exterior of $\mathbf{D}$. We denote by $\gamma^-$ and $\gamma^+$ the one-sided trace operators respectively from the interior and the exterior of $\mathbf{D}$. We use a similar convention for the normal derivatives $\partial_n^\pm$.

**Definition 4.1.** The single layer potential is the weakly singular integral operator defined for any $q \in L^2(\mathcal{C})$ by:

$$ \tilde{S}_{\mathcal{C}}q(x) = \int_{\mathcal{C}} G(x-y)q(y) \, dy \quad \text{for all } x \in \mathbb{R}^2 \setminus \mathcal{C}, $$

and extended by density as a bounded operator $\tilde{S}_{\mathcal{C}}: H^{-\frac{1}{2}}(\mathcal{C}) \to H_{\text{loc}}^{1}(\mathbb{R}^2)$. For any $q \in H^{-1/2}(\mathcal{C})$, the function $\tilde{S}_{\mathcal{C}}q$ is harmonic in $\mathbb{R}^2 \setminus \mathcal{C}$.

According to [29, p. 261], the single layer potential admits the following asymptotic expansion as $|x|$ goes to infinity:

$$ \tilde{S}_{\mathcal{C}}q(x) = -\frac{1}{2\pi} \langle q, 1 \rangle_{\mathcal{C}} \ln|x| + O(1/|x|), \quad (4.1) $$

where $\langle \cdot, \cdot \rangle_{\mathcal{C}}$ stands for the duality pairing on $H^{-\frac{1}{2}}(\mathcal{C}) \times H^{\frac{1}{2}}(\mathcal{C})$ extending the $L^2$ scalar product. For every $q \in H^{-1/2}(\mathcal{C})$, the outer trace $\gamma^+ \tilde{S}_{\mathcal{C}}q$ and inner trace $\gamma^- \tilde{S}_{\mathcal{C}}q$ coincide on $\mathcal{C}$. We denote by $S_{\mathcal{C}}q$ this common value and we recall that the operator $S_{\mathcal{C}}: H^{-1/2}(\mathcal{C}) \to H^{1/2}(\mathcal{C})$ is bounded. Regarding the inner and outer normal traces, the single layer potential satisfies the so-called “jump relation”, namely:

$$ \partial_n^- \tilde{S}_{\mathcal{C}}q - \partial_n^+ \tilde{S}_{\mathcal{C}}q = q \quad \text{for all } q \in H^{-1/2}(\mathcal{C}). \quad (4.2) $$

Let us recall the following result involving the operator $S_{\mathcal{C}}$, referred to as Theorem 8.16 in [29].

**Theorem 4.2.** If $\text{diam}(\mathbf{D}) < 1$, then $S_{\mathcal{C}}$ is an isomorphism from $H^{-\frac{1}{2}}(\mathcal{C})$ onto $H^{\frac{1}{2}}(\mathcal{C})$.

Let us emphasize that the assumption on the diameter of $\mathbf{D}$ is not restrictive (otherwise, it suffices to rescale the problem) and simply ensures that the capacity of $\mathcal{C}$ is less than 1 (see, for instance, [37, p. 143] and references therein). We shall also require the following uniqueness result (see [29, Lemma 8.14]):

**Theorem 4.3.** Given $f \in H^{\frac{1}{2}}(\mathcal{C})$ and $b \in \mathbb{R}$, there exists a unique $(q, c) \in H^{-\frac{1}{2}}(\mathcal{C}) \times \mathbb{R}$ such that

$$ S_{\mathcal{C}}q + c = f \quad \text{and} \quad \langle q, 1 \rangle_{\mathcal{C}} = b. $$

### 4.2. Main results

The main steps of our algorithm are summarized in Algorithm 4.1.

At the end of this process we then obtain a domain $\omega_N$ such that

$$ \int_{\omega} z^{\ell} dm(z) = \int_{\omega_N} z^{\ell} dm(z), \quad \text{for all } \ell = 0, \dots, 2N-1. \quad (4.3) $$

We denote by $U_{\omega_N}$ the Newtonian potential $G * \mathbb{1}_{\omega_N}$ and we claim:
**Algorithm 4.1 Summary of the method**

**Require:** $N \in \mathbb{N}$ and the Cauchy data $U_\omega$ and $\partial_n U_\omega$ on $\Gamma$, where $U_\omega = G * 1_\omega$.

(1) Using (1.4), compute the harmonic moments of $\omega$:

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

(2) Find the nodes $z_{1,N}, \dots, z_{N,N} \in \mathbb{C}$ and the weights $c_{1,N}, \dots, c_{N,N} \in \mathbb{C}$ such that

$$ \sum_{k=1}^{N} c_{k,N} z_{k,N}^{\ell} = \tau_{\ell} \quad \text{for all } \ell = 0, \dots, 2N - 1. $$

(3) Find a quadrature domain $\omega_N$ associated to $z_{1,N}, \dots, z_{N,N}$ and $c_{1,N}, \dots, c_{N,N}$ using one of the methods described Subsection 3.2 (conformal mapping) and in Subsection 3.3 (partial balayage).

**Theorem 4.4.** Assume that there exists a compact set $K$ included in the open unit disk $B$ and a constant $M > 0$ such that, for all $N \ge 1$

(1) $\omega_N \subset K$,

(2) $\omega_N \in Q(\mu_N, HL^1)$,

(3) $\sum_{k=1}^{N} |c_{k,N}| \le M$.

Then $\|\nabla U_{\omega_N} - \nabla U_{\omega}\|_{L^{\infty}(\Gamma)} \to 0$ as $N \to +\infty$.

**Proof.** Let $0 < r < 1$ be such that $\omega$ and $K$ are contained in $B_r = B(0, r)$ and denote by $\Gamma_r$ the boundary of $\partial B_r$. Let $O_r$ and $O$ be two bounded open sets such that $\Gamma_r \subset O_r$, $\Gamma \subset O$ and the sets $\bar{O}_r$, $\bar{O}$ and $K$ are pairwise disjoint (see Figure 4.1).

<figure><img src="image_6.png" /><figcaption>FIGURE 4.1. The compact $K$, the open sets $O$ and $O_r$ and the boundaries $\Gamma$ and $\Gamma_r$.</figcaption></figure>

The function $G(\cdot, \cdot)$ (we use here the notation $G(x, y)$ for $G(x-y)$) and all its derivatives are bounded on the compact set $K \times \bar{O}_r$. It follows that there exists a constant $C_0 > 0$ such that:

$$ \|G(x, \cdot)\|_{H^2(\bar{O}_r)} \le C_0, \quad \text{for all } x \in K, $$
A RECONSTRUCTION METHOD FOR THE INVERSE GRAVIMETRIC PROBLEM

and next, by continuity of the trace operator, there exists a constant $C_1 > 0$ such that:

$$ \|G(x, \cdot)\|_{H^{\frac{3}{2}}(\Gamma_r)} \le C_1, \quad \text{for all } x \in K. $$

Now, since $x \to G(x, y)$ is harmonic when $y \in \Gamma_r$, we can write (recall that $c_{k,N} \in \mathbb{R}$ from Remark 3.3)

$$ \|U_{\omega_N}\|_{H^{\frac{3}{2}}(\Gamma_r)} = \left\| \sum_{k=1}^{N} c_{k,N} G(x_{k,N}, \cdot) \right\|_{H^{\frac{3}{2}}(\Gamma_r)} \le MC_1, \quad \text{for all } N \ge 1. \quad (4.4) $$

On the other hand, according to Theorem 4.3, there exists a unique $(q_\omega, c_\omega) \in H^{-\frac{1}{2}}(\Gamma_r) \times \mathbb{R}$ such that
$$ S_{\Gamma_r} q_\omega + c_\omega = U_\omega \quad \text{on } \Gamma_r \quad \text{and} \quad \langle q_\omega, 1 \rangle_{\Gamma_r} = |\omega|. $$

Comparing (4.1) and (1.2) we deduce that:

$$ \tilde{S}_{\Gamma_r} q_\omega - U_\omega = O(1/|x|), \quad \text{as } |x| \to +\infty, \quad (4.5) $$

and thanks to [29, Theorem 8.10] (uniqueness result for the exterior Dirichlet problem), we conclude that $c_\omega = 0$ and

$$ \tilde{S}_{\Gamma_r} q_\omega(x) = U_\omega(x) \quad \text{for all } |x| > r. $$

It follows that for every harmonic function $v$ in $H^1(B)$:

$$ \int_{\Gamma_r} U_\omega \partial_n v - \partial_n U_\omega v \, d\sigma = \int_{\Gamma_r} (\tilde{S}_{\Gamma_r} q_\omega) \partial_n^- v - (\partial_n^+ \tilde{S}_{\Gamma_r} q_\omega) v \, d\sigma. $$

Integrating by parts over $B_r$ the first term in the integral, and taking into account the jump relation (4.2), we end up with:

$$ \int_{\Gamma_r} U_\omega \partial_n v - \partial_n U_\omega v \, d\sigma = \langle q_\omega, v \rangle_{\Gamma_r}. \quad (4.6a) $$

The very same computations can be carried out replacing $\omega$ with $\omega_N$ and we get as well:

$$ \int_{\Gamma_r} U_{\omega_N} \partial_n v - \partial_n U_{\omega_N} v \, d\sigma = \langle q_{\omega_N}, v \rangle_{\Gamma_r}, \quad (4.6b) $$

where $q_{\omega_N} \in H^{-1/2}(\Gamma_r)$ is such that $\tilde{S}_{\Gamma_r} q_{\omega_N}(x) = U_{\omega_N}(x)$ for every $|x| \ge r$. Setting $v_\ell(x_1, x_2) = (x_1 + ix_2)^\ell$ for every positive integer $\ell$ and every $(x_1, x_2) \in \mathbb{R}^2$, we obtain, combining (1.4) and (4.3):

$$ \langle q_\omega, v_\ell \rangle_{\Gamma_r} = \langle q_{\omega_N}, v_\ell \rangle_{\Gamma_r}, \quad \text{for all } \ell = 0, \dots, 2N-1. $$

Returning to (4.4) and invoking the compact embedding of $H^{3/2}(\Gamma_r)$ into $H^{1/2}(\Gamma_r)$, we can assume that, up to a subsequence extraction, the traces of the functions $U_{\omega_N}$ strongly converge in $H^{\frac{1}{2}}(\Gamma_r)$ as $N$ goes to $+\infty$. Since $S_{\Gamma_r}$ has a bounded inverse (as asserted by Theorem 4.2 above), the sequence $(q_{\omega_N})_N$ converges in $H^{-\frac{1}{2}}(\Gamma_r)$ to some density $q$ and for every $\ell \ge 0$ we have:

$$ \langle q_\omega, v_\ell \rangle_{\Gamma_r} = \langle q_{\omega_N}, v_\ell \rangle_{\Gamma_r}, \quad \text{for all } N > \frac{\ell+1}{2}. $$

Letting $N$ go to $+\infty$, we get:

$$ \langle q_\omega - q, v_\ell \rangle_{\Gamma_r} = 0 \quad \text{for all } \ell \ge 0. \quad (4.7) $$

The subspace of trigonometric polynomials being a dense subset in $H^{\frac{1}{2}}(\Gamma_r)$ (see [24, Theorem 8.2]), the identity (4.7) leads to $q = q_\omega$.

Notice now that the function $G(\cdot, \cdot)$ and all its derivatives are bounded on the compact set $\bar{\mathcal{O}} \times \bar{\mathcal{O}}_r$, which entails the existence of a constant $C_2 > 0$ such that:

$$ \|\nabla_x G(x, \cdot)\|_{H^1(\mathcal{O}_r)} \le C_2 \quad \text{for all } x \in \bar{\mathcal{O}}. $$
By continuity of the trace operator, it follows that there exists a constant $C_3 > 0$ such that:

$$|(\nabla U_{\omega_N} - \nabla U_{\omega})(x)| = \left|\int_{\Gamma_r} \nabla_x G(x, y)(q_{\omega_N} - q_{\omega})(y) d\sigma_y\right| \\ \le \|\nabla_x G(x, \cdot)\|_{H^{1/2}(\Gamma_r)} \|q_{\omega_N} - q_{\omega}\|_{H^{-1/2}(\Gamma_r)} \\ \le C_3 \|q_{\omega_N} - q_{\omega}\|_{H^{-1/2}(\Gamma_r)}.$$

The conclusion follows by letting $N \to +\infty$ and noticing, by a classical argument, that the uniqueness of the limit ensures that the convergence holds for the whole sequence $(U_{\omega_N})_N$ and not only for a subsequence. ■

**Remark 4.5.** Condition 3 in Theorem 4.4 holds, for example, when all the $c_{k,N}$ are positive, since

$$\sum_{k=1}^{N} c_{k,N} = |\omega|, \quad \text{for all } N \ge 1.$$

As expected, the previous theorem gives a convergence result only for the gravitational fields. However, adding the hypotheses of Theorem 1.3 we can prove the following result regarding the domains.

**Corollary 4.6.** Under the hypotheses of Theorem 4.4 assume in addition that $\omega$ and $\omega_N$ both satisfy (S) from Section 1.2 and that their polar representations $\rho_\omega$ and $\rho_{\omega_N}$ belong to the set $\mathcal{R}$ for all $N \ge 1$, where $\mathcal{R}$ is defined in Theorem 1.3. Then, $(\omega_N)_N$ converges to $\omega$ in the sense of

$$\rho_{\omega_N} \xrightarrow[N \to +\infty]{} \rho_\omega, \quad \text{in } L^\infty.$$

**Proof.** Note that, thanks to (4.3), $\omega$ and all the $\omega_N$ have the same center of gravity and then apply Theorems 1.3 and 4.4. ■

In the special case where the unknown domain is a subharmonic quadrature domain the above result can be strongly improved, as an exact reconstruction of $\omega$ can be achieved in a finite number of steps.

**Theorem 4.7.** Assume that $\omega$ is a subharmonic quadrature domain associated to a finite sum of point masses. Let $(\omega_N)_N$ be the sequence of subharmonic quadrature domains constructed with the algorithm 4.1, using partial balayage. Then there exists a positive integer $M$ such that $\omega_M = \omega$.

**Proof.** By assumption, there exists a positive integer $M$, nodes $z_1, \dots, z_M$ and weights $c_1, \dots, c_M$ such that $\omega \in Q(\mu, SL^1)$ with $\mu = \sum_{k=1}^M c_k \delta_{z_k}$. Next, by Corollary 2.3, applying step (2) of Algorithm 4.1, with $N = M$ yields exactly the nodes $z_1, \dots, z_M$ and weights $c_1, \dots, c_M$. Then, step (3) generates $\omega_M \in Q(\mu, SL^1)$, by Proposition 3.7, and the claimed result follows by uniqueness in $Q(\mu, SL^1)$. ■

According to [7, Lemma 2], the value of $M$ can be determined by noticing that:

$$\mathrm{rank}(\mathbb{H}_0^{(n)}) = M, \quad \forall n \ge M, \qquad (4.8)$$

where $\mathbb{H}_0^{(n)}$ is defined in (2.5a).

**Remark 4.8.** Theoretically, increasing the number $N$ of nodes in our algorithm should increase the accuracy of the reconstruction. In practice, it also increases the numerical errors (mainly in the resolution of Prony's system). Hence taking $N$ as large as possible does not always result in a successful strategy. Inspired by (4.8), we could take for $N$ the number of “significant” singular values of the matrix $\mathbb{H}_0^{(n)}$, with $n$ large enough.
## 5. Numerical experiments

We present in this section some numerical experiments. For every simulation, we choose a domain $\omega$ to be reconstructed and compute its gravitational potential $U_\omega$ and gravitational field $\nabla U_\omega$ to which we apply a white gaussian noise. These data are then taken as inputs in Algorithm 4.1 for different values of $N$.

The harmonic moments are computed from these boundary data using the Fast Fourier Transform (FFT). In order to reduce the impact of noise on the measurements, we use a large number of quadrature points to evaluate the moments.

From the harmonic moments, we construct the Hankel matrices (2.5) and solve the corresponding Prony's system which provides the (complex) nodes $z_{1,N}, \dots, z_{N,N}$ and (complex) weights $c_{1,N}, \dots, c_{N,N}$.

Finally we compute a quadrature domain using one of the two methods previously introduced (in Section 3.2 and Section 3.3). The results are presented in the following two subsections.

### 5.1. Conformal mapping

As this method presents convergence issues, we restrict our experiments to the cases where $\omega$ satisfies Property (S) and we use noiseless data. The boundary of the target domain is parameterized in polar coordinates by a function $\rho : ] -\pi, \pi] \to \mathbb{R}_+$ taking the form of a Fourier series with a finite number of modes:

$$ \rho(\theta) = 1 + \sum_{k=1}^{n} a_k \cos(k\theta) + b_k \sin(k\theta) \quad \text{for all } \theta \in ] -\pi, \pi]. $$

We test the conformal mapping reconstruction method with the domains depicted on Figure 5.1. We represent on the first row of Figure 5.2 and Figure 5.3, the union of the disks centered at the nodes $z_{k,N}$ and with radii $\sqrt{|\text{Re}(c_{k,N})|}/\pi$ (if $\text{Re}(c_{k,N}) \le 0$, we color the corresponding disk in red). On the second row, we plot the image of the unit disk by the computed conformal mapping $\varphi_N$, as introduced in (3.7), for several values of $N$. We can notice that the quality of the reconstruction improves with $N$ and that we obtain a fairly accurate approximation of the target shape.

However, let us emphasize that the nonlinear solver is quite sensitive to the shape of the target domain and to the initial guess. For example, on Figure 5.2d and Figure 5.3e, the algorithm failed to converge. It is not clear whether this failure is due to the lack of existence of solutions to System (3.8) or to a poor choice of initial guess. Nevertheless, one may notice that a partial balayage applied to the measure obtained by replacing the complex weights by their real parts, would provide a fairly good approximation of the target shape. This can be seen by considering the disks representation (see next subsection).

<figure><img src="image_7.png" /><figcaption>FIGURE 5.1. Two star-shaped domains to be reconstructed.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIGURE 5.2. Reconstruction of the domain presented in Figure 5.1a</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIGURE 5.3. Reconstruction of the domain presented in Figure 5.1b</figcaption></figure>

Besides its sensitivity to the initial guess, this algorithm is limited to the cases of solid domains. Its main advantage is that it allows for non-real (and non-positive) weights. However, in our tests, the imaginary parts of the weights were systematically close to zero. For all these reasons, in most of the cases, the partial balayage method should be preferred.

## 5.2. Convex minimization

We focus hereafter on the partial balayage method (described in Section 3.3). Since we need the computed weights $c_{k,N}$ to be real, we replace them by their real part if necessary. In our tests, though, this choice turned out to be of little importance because we observed that the imaginary part was always very small with respect to the real part. We recall that the problem to be solved is summarized in (3.16). We follow the approach described in [23] to solve the variational inequality involved in this
A RECONSTRUCTION METHOD FOR THE INVERSE GRAVIMETRIC PROBLEM

convex minimization problem and we implemented this algorithm using FreeFEM++ [17]. We use $P_2$ finite elements in order to compute the laplacian of the minimizer, and hence evaluate an approximation of the operator $F$. In the tests of this subsection, we use the regularization parameter $\alpha = 0.99$.

We begin with a domain having the property (S) (see Figure 5.4). On each picture is displayed the targeted domain (in red), the disks $(B(z_k, r_k))_k$ associated with the solution of Prony's system (the centers are materialized by small markers), and a finite element approximation of $F(\tilde{u}_N^\alpha)$, for different values of $N$. These tests are performed without noise and allow a fairly accurate reconstruction of the domain. One may notice that some centers lay outside the domain (for $N = 15$), but they correspond to nodes for which the weight is zero.

<figure><img src="image_10.png" /><figcaption>FIGURE 5.4. Reconstruction of a star-shaped domain (without noise).</figcaption></figure>

We wish now to investigate the influence of noise through some examples: simple ones (the previous star-shaped domain in Figure 5.5, a domain convex in one direction in Figure 5.6 and a set of disjoint disks in Figure 5.7) and more sophisticated ones (a rabbit-shaped domain in Figure 5.8 and Figure 5.9, and a plane-shaped one in Figure 5.10 and Figure 5.11). For these two last examples, we represent only the tests with a reasonable ratio between the truncation parameter $N$ and the noise. For all these cases, even in the presence of noise, we observe that the algorithm generates very few positive mass outside the targeted domain. Although the noise has a deteriorating effect on the details, the raw reconstruction of the domain contour seems to be relatively robust. In Figure 5.7, note the good quality of the reconstruction on the last image, which is unusual at this noise level.

One limitation of our method is the condition number of the first step, as invoked in Section 2.3. Indeed, for some values of $N$ the computations may conduct to totally wrong reconstructions. We illustrate this fact in Figure 5.12, for which we deal with the same unknown domain as in Figure 5.6. We remark that in Figure 5.12a, all the weights are almost equal to zero. Of course, searching for an associated quadrature domain do not yield a satisfactory approximation of $\omega$. However, such situations can be avoided by choosing another value for $N$, as visible in Figure 5.12b.
<figure><img src="image_11.png" /><figcaption>FIGURE 5.5. Reconstruction of a star-shaped domain with noisy data (N = 15)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 5.6. Reconstruction of a domain convex in one direction with noisy data (N = 12)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 5.7. Reconstruction of disks with noisy data (N = 10)</figcaption></figure>

## 6. Conclusion

In this work, we propose an algorithm to solve the planar gravitational inverse problem. The approach is based on an original combination of Prony's method and generation of quadrature domains. Under some technical hypotheses, we are able to prove that the sequence of reconstructed domains converges
<figure><img src="image_14.png" /><figcaption>FIGURE 5.8. A rabbit-shaped domain.</figcaption></figure>

<table><caption>FIGURE 5.9. Reconstruction of the rabbit-shaped domain shown in Figure 5.8</caption><thead><tr><th></th><th>N = 10</th><th>N = 15</th><th>N = 20</th><th>N = 25</th></tr></thead><tbody><tr><th>0 %</th><td></td><td></td><td></td><td></td></tr><tr><th>1 %</th><td></td><td></td><td></td><td></td></tr><tr><th>3 %</th><td></td><td></td><td></td><td></td></tr><tr><th>5 %</th><td></td><td></td><td></td><td></td></tr></tbody></table>
<figure><img src="image_15.png" /><figcaption>FIGURE 5.10. A plane-shaped domain.</figcaption></figure>

<table><caption>FIGURE 5.11. Reconstruction of the plane-shaped domain shown in Figure 5.10</caption><thead><tr><th></th><th>N = 10</th><th>N = 15</th><th>N = 20</th><th>N = 25</th></tr></thead><tbody><tr><th>0 %</th><td></td><td></td><td></td><td></td></tr><tr><th>1 %</th><td></td><td></td><td></td><td></td></tr><tr><th>3 %</th><td></td><td></td><td></td><td></td></tr><tr><th>5 %</th><td></td><td></td><td></td><td></td></tr></tbody></table>
<figure><img src="image_16.png" /><figcaption>FIGURE 5.12. False reconstruction with $N = 31$ (all the weights are almost zero)</figcaption></figure>

to a domain which is graviequivalent to the targeted one (we recall that the gravitational inverse problem is severely ill-posed because even the uniqueness is not guaranteed in general). The numerical tests suggest that the hypotheses required for the algorithm to converge are generally satisfied. The reconstructed domains are quite convincing even with noisy data.

We conjecture that the assumptions of Theorem 4.4 can be lifted. Nonetheless, the first one (asserting that the nodes deriving from Prony's system have to remain in a compact domain) is at the crossroads of several research areas in complex analysis (such as Padé approximants) and it seems therefore rather difficult to get around it.

Extension to the three dimensional case is uncertain. Although the partial balayage method works the same in 3D, we can no longer rely on Prony's method which is a purely complex analysis tool. Another generalization is to consider non uniform and possibly surface mass distribution. This seems to be within reach and require only minor adaptation of the method.

## Appendix A. Appendix

**Proof of Proposition 2.2.** Assume that $(z_1, \dots, z_p)$ and $(c_1, \dots, c_p)$ is a solution of Problem (2.4) with $p \le N$. Then, one easily checks that:

$$ \mathbb{H}_0^{(N)} = VDV^T, \qquad (A.1) $$

where $V$ and $D$ are matrices of size $N \times N$ with $D = \text{diag}(c_1, \dots, c_p, 0, \dots, 0)$ and:

$$ V = \begin{pmatrix} 1 & \dots & 1 & 1 & \dots & 1 \\ z_1 & \dots & z_p & z_p & \dots & z_p \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ z_1^{N-1} & \dots & z_p^{N-1} & z_p^{N-1} & \dots & z_p^{N-1} \end{pmatrix}. $$

Since $\mathbb{H}_0^{(N)}$ is invertible, $V$ and $D$ are necessarily invertible as well, and thus $p = N$.

Concerning uniqueness, we follow the lines of the proof given in [30]. Let $(z_1, \dots, z_N)$ be a solution of Problem (2.4) with $p = N$ and define the polynomial

$$ Q(z) = \prod_{j=1}^{N} (z - z_j) = z^N + \sum_{j=1}^{N} q_j z^{N-j}. $$

One can verify that

$$ \sum_{k=0}^{N-1} \tau_{k+i} q_{N-k} = -\tau_{N+i} \quad \text{for all } i = 0, \dots, N-1. $$
These equations can be reformulated as the following linear system:

$$ \mathbb{H}_0^{(N)} q = - \begin{pmatrix} \tau_N \\ \vdots \\ \tau_{2N-1} \end{pmatrix} \quad \text{with} \quad q = \begin{pmatrix} q_1 \\ \vdots \\ q_N \end{pmatrix}. $$

Since $\mathbb{H}_0^{(N)}$ is invertible, the vector $q$ is uniquely determined and hence this is also true for $Q$ and for its roots ($z_1, \dots, z_N$). Uniqueness for the weights follows because they are solution of a linear system with an invertible Vandermonde matrix. ■

**Proof of Corollary 2.3.** Assume that there exists a solution with $p = N$ to Problem (2.4). In that case, we deduce from (A.1) that $\mathbb{H}_0^{(N)}$ is invertible. The conclusion follows from Proposition 2.2. ■

We briefly recall Prony's algorithm in Algorithm 1.1.

<caption>**Algorithm 1.1 Prony's algorithm**</caption>

**Require:** $N \in \mathbb{N}$ and $\tau_0, \dots, \tau_{2N-1}$

(1) Construct the Hankel matrix $\mathbb{H}_0^{(N)}$ (defined in (2.5a)) and solve the linear system:

$$ \mathbb{H}_0^{(N)} \begin{pmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_{N-1} \end{pmatrix} = - \begin{pmatrix} \tau_N \\ \tau_{N+1} \\ \vdots \\ \tau_{2N-1} \end{pmatrix}. $$

(2) Find the $N$ roots $z_1, \dots, z_N$ of the polynomial:

$$ z^N + \sum_{j=0}^{N-1} \beta_j z^j. $$

(3) Find $c_1, \dots, c_N$ by solving the linear system:

$$ \begin{pmatrix} 1 & \cdots & 1 \\ z_1 & \cdots & z_N \\ \vdots & & \vdots \\ z_1^{N-1} & \cdots & z_N^{N-1} \end{pmatrix} \begin{pmatrix} c_1 \\ c_2 \\ \vdots \\ c_N \end{pmatrix} = \begin{pmatrix} \tau_0 \\ \tau_1 \\ \vdots \\ \tau_{N-1} \end{pmatrix}. $$

**Proof of Theorem 2.4.** Assume that Problem (2.4) admits a solution with $p = N$ and denote by $z_1, \dots, z_N$ the nodes. We introduce $\beta_0, \beta_1, \dots, \beta_{N-1}$ such that:

$$ R(z) = \prod_{j=1}^{N} (z - z_j) = \sum_{k=0}^{N-1} \beta_k z^k + z^N. $$

In the other hand, let

$$ M(z) = \begin{pmatrix} \tau_0 & \tau_1 & \cdots & \tau_{N-1} & \tau_N \\ \tau_1 & \tau_2 & \cdots & \tau_N & \tau_{N+1} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ \tau_{N-1} & \tau_N & \cdots & \tau_{2N-2} & \tau_{2N-1} \\ 1 & z & \cdots & z^{N-1} & z^N \end{pmatrix} \quad \text{for all } z \in \mathbb{C}. $$
A RECONSTRUCTION METHOD FOR THE INVERSE GRAVIMETRIC PROBLEM

One can verify that

$$ M(z_k) \begin{pmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_{N-1} \\ 1 \end{pmatrix} = 0, \quad \text{for all } k = 1, \dots, N $$

and thus,

$$ \det(M(z_k)) = 0 = P_N(z_k), $$

which proves that $z_k$ are the roots of $P_N$ (remark that the degree of $P_N$ is $N$).

Conversely suppose that $P_N$ admits $N$ simple roots $z_1, \dots, z_N$. Keeping the same notation, we have:

$$ M(z_k) = \left( \begin{array}{c|c} \mathbb{H}_0^{(N)} & \begin{array}{c} \tau_N \\ \tau_{N+1} \\ \vdots \\ \tau_{2N-1} \end{array} \\ \hline 1 & z_k \cdots z_k^{N-1} \end{array} \right). $$

Observe that $\deg(P_N) = N$ and since its leading coefficient is $(-1)^N \det \mathbb{H}_0^{(N)}$, we deduce that $\det \mathbb{H}_0^{(N)} \neq 0$. In addition,

$$ \ker(M(z_k)) \neq \{0\}, \quad \text{for all } k = 1, \dots, N, $$

and we can verify that every nonzero element of $\ker(M(z_k))$ is proportional to $(\beta_0, \dots, \beta_{N-1}, 1)$ with:

$$ \begin{pmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_{N-1} \end{pmatrix} = -[\mathbb{H}_0^{(N)}]^{-1} \begin{pmatrix} \tau_N \\ \tau_{N+1} \\ \vdots \\ \tau_{2N-1} \end{pmatrix}. \quad (A.2) $$

Now, writing

$$ M(z_k) \begin{pmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_{N-1} \\ 1 \end{pmatrix} = 0, $$

and using the last line of this system, we obtain that $z_1, \dots, z_N$ are the roots of the polynomial

$$ \sum_{k=0}^{N-1} \beta_k z^k + z^N. $$

Since the coefficients of this polynomial satisfy (A.2), we can apply Algorithm 1.1 to construct a solution for which $p = N$, whose nodes are exactly $z_1, \dots, z_N$. Finally, since $\det(\mathbb{H}_0^{(N)}) \neq 0$, the uniqueness follows from Proposition 2.2. ■

References

[1] Yacin Ameur, Martin Helmer, and Felix Tellander. On the uniqueness problem for quadrature domains. *Comput. Methods Funct. Theory*, 21(3):473–504, 2021.

[2] Kendall Atkinson and Weimin Han. *Theoretical numerical analysis. A functional analysis framework*, volume 39 of *Texts in Applied Mathematics*. Springer, second edition, 2005.
[3] Dmitry Batenkov and Yosef Yomdin. On the accuracy of solving confluent Prony systems. *SIAM J. Appl. Math.*, 73(1):134–154, 2013.

[4] Steven R. Bell. Density of quadrature domains in one and several complex variables. *Complex Var. Elliptic Equ.*, 54(3-4):165–171, 2009.

[5] Alfredo Canelas, Antoine Laurain, and Antonio A. Novotny. A new reconstruction method for the inverse potential problem. *J. Comput. Phys.*, 268:417–431, 2014.

[6] Victor G. Cherednichenko. *Inverse logarithmic potential problem*. Inverse and Ill-Posed Problems Series. Walter de Gruyter, 1996.

[7] Abdellatif El Badia and Tuong Ha-Duong. An inverse source problem in potential analysis. *Inverse Probl.*, 16(3):651–663, 2000.

[8] Avner Friedman. Variational principles and free-boundary problems, 1982.

[9] Stephen J. Gardiner and Tomas Sjödin. Convexity and the exterior inverse problem of potential theory. *Proc. Am. Math. Soc.*, 136(5):1699–1703, 2008.

[10] Stephen J. Gardiner and Tomas Sjödin. Partial balayage and the exterior inverse problem of potential theory. In *Potential theory and stochastics in Albac. Aurel Cornea memorial volume. Conference proceedings, Albac, Romania, September 4–8, 2007*, volume 11 of *Theta Series in Advanced Mathematics*, pages 111–123. The Theta Foundation, 2008.

[11] Gene H. Golub, Peyman Milanfar, and James Varah. A stable numerical method for inverting shape from moments. *SIAM J. Sci. Comput.*, 21(4):1222–1243, 1999.

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

[13] Björn Gustafsson, Chiyu He, Peyman Milanfar, and Mihai Putinar. Reconstructing planar domains from their moments. *Inverse Probl.*, 16(4):1053–1070, 2000.

[14] Björn Gustafsson and Makoto 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.

[15] Björn Gustafsson and Henrik Shahgholian. Existence and geometric properties of solutions of a free boundary problem in potential theory. *J. Reine Angew. Math.*, 473:137–179, 1996.

[16] Björn Gustafsson and Harold S. Shapiro. What is a quadrature domain? In *Quadrature domains and their applications*, volume 156 of *Operator Theory: Advances and Applications*, pages 1–25. Birkhäuser, 2005.

[17] Freidrich Hecht. New development in FreeFem++. *J. Numer. Math.*, 20(3-4):251–265, 2012.

[18] Frank Hettlich and William Rundell. Iterative methods for the reconstruction of an inverse potential problem. *Inverse Probl.*, 12(3):251–266, 1996.

[19] Jeffrey Hokanson. *Numerically Stable and Statistically Efficient and Algorithms for Large and Scale Exponential and Fitting*. PhD thesis, Rice University, 2013. https://hdl.handle.net/1911/77161.

[20] Victor Isakov. *Inverse source problems*, volume 34 of *Mathematical Surveys and Monographs*. American Mathematical Society, 1990.

[21] Victor Isakov, Shingyu Leung, and Jianliang Qian. A fast local level set method for inverse gravimetry. *Commun. Comput. Phys.*, 10(4):1044–1070, 2011.

[22] Victor Isakov, Shingyu Leung, and Jianliang Qian. A three-dimensional inverse gravimetry problem for ice with snow caps. *Inverse Probl. Imaging*, 7(2):523–544, 2013.

[23] Kazufumi Ito and Karl Kunisch. Semi-Smooth Newton Methods for Variational Inequalities of the First Kind. *ESAIM, Math. Model. Numer. Anal.*, 37(1):41–62, 2003.

[24] Rainer Kress. *Linear integral equations*, volume 82 of *Applied Mathematical Sciences*. Springer, third edition, 2014.
[25] Rainer Kress and William Rundell. A nonlinear integral equation and an iterative algorithm for an inverse source problem. *J. Integral Equations Appl.*, 27(2):179–197, 2015.

[26] Wenbin Li, Wangtao Lu, and Jianliang Qian. A level-set method for imaging salt structures using gravity data. *Geophysics*, 81(2):G27–G40, 2016.

[27] Wenbin Li and Jianliang Qian. Simultaneously recovering both domain and varying density in inverse gravimetry by efficient level-set methods. *Inverse Probl. Imaging*, 15(3):387–413, 2021.

[28] Wangtao Lu, Shingyu Leung, and Jianliang Qian. An improved fast local level set method for three-dimensional inverse gravimetry. *Inverse Probl. Imaging*, 9(2):479–509, 2015.

[29] William McLean. *Strongly elliptic systems and boundary integral equations*. Cambridge University Press, 2000.

[30] P. Milanfar, G. C. Verghese, W. C. Karl, and A. S. Willsky. Reconstructing polygons from moments with connections to array processing. *IEEE Trans. Signal Process.*, 43(2):432–443, 1995.

[31] V. F. Pisarenko. The Retrieval of Harmonics from a Covariance Function. *Geophys. J. Int.*, 33(3):347–366, 1973.

[32] Roland Potthast. A survey on sampling and probe methods for inverse problems. *Inverse Probl.*, 22(2):R1–R47, 2006.

[33] Nilson C. Roberty and Carlos J. S. Alves. On the identification of star-shape sources from boundary measurements using a reciprocity functional. *Inverse Probl. Sci. Eng.*, 17(2):187–202, 2009.

[34] Makoto Sakai. *Quadrature domains*, volume 934 of Lecture Notes in Mathematics. Springer, 1982.

[35] Tetsuya Sakurai, Junko Asakura, Hiroto Tadano, and Tsutomu Ikegami. Error analysis for a matrix pencil of Hankel matrices with perturbed complex moments. *JSIAM Lett.*, 1:76–79, 2009.

[36] Harold S. Shapiro. *The Schwarz function and its generalization to higher dimensions*. John Wiley & Sons, 1992.

[37] Olaf Steinbach. *Numerical approximation methods for elliptic boundary value problems*. Springer, 2008.

[38] Petre Stoica and Randolph L. Moses. *Spectral Analysis of Signals*. Pearson Prentice Hall, 2005.

[39] Ben Stray, Andrew Lamb, Aisha Kaushik, Jamie Vovrosh, Anthony Rodgers, Jonathan Winch, Farzad Hayati, Daniel Boddice, Artur Stabrawa, Alexander Niggebaum, Mehdi Langlois, Yu-Hung Lien, Samuel Lellouch, Sanaz Roshanmanesh, Kevin Ridley, Geoffrey de Villiers, Gareth Brown, Trevor Cross, George Tuckwell, Asaad Faramarzi, Nicole Metje, Kai Bongs, and Michael Holynski. Quantum sensing for gravity cartography. *Nature*, 602(7898):590–594, 2022.

[40] K. Doel, Uri M. Ascher, and Antonio Leitão. Multiple level sets for piecewise constant surface reconstruction in highly ill-posed problems. *J. Sci. Comput.*, 43(1):44–66, 2010.

[41] Dimiter Zidarov. *Inverse gravimetric problem in geoprospecting and geodesy*. *Transl. from the Bulgarian*. Elsevier, 1990.