Direct algorithm for reconstructing small absorbers in thermoacoustic tomography problem from a single data

First, we set

$$\begin{equation*}h\left(\mathbf{\text{x}},t\right)=\beta \left(\mathbf{\text{x}}\right)\frac{\partial H}{\partial t}\left(\mathbf{\text{x}},t\right),\end{equation*}$$

where H is given in (1.2), and we consider the adjoint problem

$$\begin{equation}\begin{matrix}{c}\hfill \frac{1}{{c}_{\text{s}}^{2}}{v}_{tt}-{\Delta}v=0\quad \text{in} {\Omega}{\times}\left(0,T\right).\hfill \end{matrix}\end{equation} \tag{ 2.1 }$$

Then, multiplying equation (1.3) by v, and integrating by parts, we obtain using Green's formula

$$\begin{align}\hfill {\int }_{0}^{T}{\int }_{{\Omega}}h\left(\mathbf{\text{x}},t\right)v\left(\mathbf{\text{x}},t\right)\mathrm{d}\mathbf{\text{x}}\mathrm{d}t& ={\int }_{{\Sigma}}f\left(\xi ,t\right)\frac{\partial v}{\partial \nu }\left(\xi ,t\right)\;\mathrm{d}\xi \mathrm{d}t+\frac{1}{{c}_{\text{s}}^{2}}{\int }_{{\Omega}}[{p}_{t}\left(\mathbf{\text{x}},T\right)v\left(\mathbf{\text{x}},T\right)\hfill \\ \hfill & -p\left(\mathbf{\text{x}},T\right){v}_{t}\left(\mathbf{\text{x}},T\right)]\;\mathrm{d}\mathbf{\text{x}}.\hfill \end{align} \tag{ 2.2 }$$

We note that

$$\begin{equation*}h\left(\mathbf{\text{x}},t\right)=\beta \left(\mathbf{\text{x}}\right)\frac{\partial H}{\partial t}\left(\mathbf{\text{x}},t\right)=\beta \left(\mathbf{\text{x}}\right)\sigma \left(\mathbf{\text{x}}\right)\frac{\partial }{\partial t}\vert u\left(\mathbf{\text{x}},t\right){\vert }^{2}.\end{equation*}$$

Let

$$\begin{equation*}g\left(\mathbf{\text{x}},t\right){:=}\beta \left(\mathbf{\text{x}}\right)\frac{\partial }{\partial t}\vert u\left(\mathbf{\text{x}},t\right){\vert }^{2},\end{equation*}$$

then

$$\begin{equation*}h\left(\mathbf{\text{x}},t\right)=\sigma \left(\mathbf{\text{x}}\right)g\left(\mathbf{\text{x}},t\right)=\sum _{j=0}^{m}{\sigma }_{j}\left(\mathbf{\text{x}}\right)g\left(\mathbf{\text{x}},t\right){\chi }_{Dj}.\end{equation*}$$

If we denote the so called reciprocity gap functional by

$$\begin{equation}\mathcal{R}\left(f,v\right)={\int }_{{\Sigma}}f\left(\xi ,t\right)\frac{\partial v}{\partial \nu }\left(\xi ,t\right)\;\mathrm{d}\xi \mathrm{d}t+\frac{1}{{c}_{\text{s}}^{2}}{\int }_{{\Omega}}\left[{p}_{t}\left(\mathbf{\text{x}},T\right)v\left(\mathbf{\text{x}},T\right)-p\left(\mathbf{\text{x}},T\right){v}_{t}\left(\mathbf{\text{x}},T\right)\right]\;\mathrm{d}\mathbf{\text{x}},\end{equation} \tag{ 2.3 }$$

then, from (2.2), we get

$$\begin{equation*}\mathcal{R}\left(f,v\right)=\sum _{j=0}^{m}{\int }_{0}^{T}{\int }_{{D}_{j}}{\sigma }_{j}\left(\mathbf{\text{x}}\right)g\left(\mathbf{\text{x}},t\right)v\left(\mathbf{\text{x}},t\right)\;\mathrm{d}\mathbf{\text{x}}\mathrm{d}t.\end{equation*}$$

Furthermore, the large difference in conductivity between healthy tissues and cancerous ones, see [6, 7, 22, 23], permits us to assume that σ₀ = O(^κ) for some κ > 0, assumed in our work to be κ > 4. Consequently, we have

$$\begin{equation}\mathcal{R}\left(f,v\right)=\sum _{j=1}^{m}{\int }_{0}^{T}{\int }_{{D}_{j}}{\sigma }_{j}\left(\mathbf{\text{x}}\right)g\left(\mathbf{\text{x}},t\right)v\left(\mathbf{\text{x}},t\right)\;\mathrm{d}\mathbf{\text{x}}\mathrm{d}t+O\left({{\epsilon}}^{\kappa }\right).\end{equation} \tag{ 2.4 }$$

To solve our inverse problem, we need to know the values of p(x, T) and p_t(x, T), unless we impose that the solution v of (2.1) satisfies

$$\begin{equation}v\left(\mathbf{\text{x}},T\right)={v}_{t}\left(\mathbf{\text{x}},T\right)=0\quad \text{in} {\Omega}.\end{equation} \tag{ 2.5 }$$

To overcome this difficulty, and for a source term of the form ${\sum }_{j=1}^{N}{\lambda }_{j}\left(t\right){\delta }_{{\mathbf{S}}_{j}}$, the authors in [14] opted for the determination of the values of the unknowns p_t(x, T) and p(x, T) using the fact that the intensities λ_j(t) become null after a certain time T* < T and by applying the exact controllability approach. In here, this approach is not employed where we opt for the choice of the condition (2.5). Indeed, we introduce what one calls test functions, the solutions of (2.1) verifying the condition (2.5).

In the following, we will show how an appropriate choice of test functions unveils the desired information. The general idea behind our identification method is the projection of the problem onto chosen test functions. This idea is not new, but our approach employs a specific family of functions that allows a practical solution.Let ${\alpha }_{i},{\beta }_{i},{\gamma }_{i}\in \mathbb{R}$, for i = 1, 2 such that

$$\begin{equation*}{\Omega}\subset \left[{\alpha }_{1},{\alpha }_{2}\right]{\times}\left[{\beta }_{1},{\beta }_{2}\right]{\times}\left[{\gamma }_{1},{\gamma }_{2}\right].\end{equation*}$$

Consider now the one-dimensional wave equation,

$$\begin{equation}\left\{\begin{aligned}\hfill & \frac{1}{{c}_{s}^{2}}{\rho }_{tt}\left(z,t\right)-{\rho }_{zz}\left(z,t\right)=0\text{in}\left({\gamma }_{1},{\gamma }_{2}\right){\times}\left(0,T\right)\hfill \\ \hfill & \rho \left(z,T\right)={\rho }_{t}\left(z,T\right)=0\text{in}\left({\gamma }_{1},{\gamma }_{2}\right).\hfill \end{aligned}\right.\end{equation} \tag{ 2.6 }$$

Then, employing the test functions

$$\begin{equation}{v}_{n}^{a}\left(x,y,z,t\right)=\rho \left(z,t\right){\left(x+iy\right)}^{n},\quad n\in \mathbb{N},\end{equation} \tag{ 2.7 }$$

the relation (2.4) becomes

$$\begin{equation}\mathcal{R}\left(f,{v}_{n}^{a}\right)=\sum _{j=1}^{m}{\int }_{0}^{T}{\int }_{{D}_{j}}{\sigma }_{j}\left(\mathbf{\text{x}}\right)g\left(\mathbf{\text{x}},t\right)\rho \left(z,t\right){\left(x+iy\right)}^{n}\mathrm{d}\mathbf{\text{x}}\mathrm{d}t+O\left({{\epsilon}}^{\kappa }\right)\quad \text{for} \text{all} n\in \mathbb{N}.\end{equation} \tag{ 2.8 }$$

Thus, using the change of variables x = S_j + τ with τ = (τ₁, τ₂, τ₃), one obtains

$$\begin{equation}\mathcal{R}\left(f,{v}_{n}^{a}\right)=\sum _{j=1}^{m}{\int }_{0}^{T}{\int }_{{B}_{j}}{{\epsilon}}^{3}{\gamma }_{j}\left(\boldsymbol{\tau },t\right)\rho \left({z}_{j}+{\epsilon}{\tau }_{3},t\right){\left[{x}_{j}+i{y}_{j}+{\epsilon}\left({\tau }_{1}+i{\tau }_{2}\right)\right]}^{n}\mathrm{d}\boldsymbol{\tau }\mathrm{d}t+O\left({{\epsilon}}^{\kappa }\right),\end{equation} \tag{ 2.9 }$$

where

$$\begin{equation}{\gamma }_{j}\left(\boldsymbol{\tau },t\right)={\sigma }_{j}\left({\mathbf{\text{S}}}_{j}+{\epsilon}\boldsymbol{\tau }\right)g\left({\mathbf{\text{S}}}_{j}+{\epsilon}\boldsymbol{\tau },t\right),\end{equation} \tag{ 2.10 }$$

and consequently,

$$\begin{equation}\mathcal{R}\left(f,{v}_{n}^{a}\right)=\sum _{j=1}^{m}\sum _{k=0}^{n}{\nu }_{j}^{k,a}\left(\genfrac{}{}{0.0pt}{}{n}{k}\right){\left({P}_{j}^{a}\right)}^{n-k}+O\left({{\epsilon}}^{\kappa }\right),\quad \text{for} \text{all} n\in \mathbb{N},\end{equation} \tag{ 2.11 }$$

where

$$\begin{equation}{\nu }_{j}^{k,a}={{\epsilon}}^{3+k}{\int }_{0}^{T}{\int }_{{B}_{j}}{\gamma }_{j}\left(\boldsymbol{\tau },t\right)\rho \left({z}_{j}+{\epsilon}{\tau }_{3},t\right){\left[{\tau }_{1}+i{\tau }_{2}\right]}^{k}\mathrm{d}\boldsymbol{\tau }\mathrm{d}t,\end{equation} \tag{ 2.12 }$$

$$\begin{equation*}{P}_{j}^{a}={x}_{j}+i{y}_{j},\end{equation*}$$

and

$$\begin{equation*}\left(\genfrac{}{}{0.0pt}{}{n}{k}\right)=\begin{cases}\frac{n!}{k!\left(n-k\right)!}\text{if}n{\geqslant}k\quad \hfill \\ 0\text{if}n{< }k.\hfill \end{cases}\end{equation*}$$

Our goal now is to solve the equation (2.11) for the reconstruction of the unknown coefficients ${\nu }_{j}^{k,a}$ and the centers ${P}_{j}^{a}$. We observe that $\forall {n}_{0}\in \mathbb{N}$, the system (2.11), with n = 0, ..., n₀, has n₀ + 1 equations and m(n₀ + 2) unknowns. Then, the number of unknowns is variable with n₀, and it increases with the increase of the number of equations. Thus, it is impossible to uniquely reconstruct the unknowns no matter how many equations are taken. To overcome this difficulty, we will truncate these equations from a small error. First, for a given positive < 1, we choose a fixed positive integer K such that K < κ − 4 and ^K+4 is small enough, and we set

$$\begin{equation}{\alpha }_{n}^{a}=\sum _{j=1}^{m}\sum _{k=0}^{K}{\nu }_{j}^{k,a}\left(\genfrac{}{}{0.0pt}{}{n}{k}\right){\left({P}_{j}^{a}\right)}^{n-k},\quad \text{for} \text{all}\;n\in \mathbb{N}.\end{equation} \tag{ 2.13 }$$

Then, according to (2.11), we can see that for n ⩽ K one has

$$\begin{equation*}\mathcal{R}\left(f,{v}_{n}^{a}\right)={\alpha }_{n}^{a}+O\left({{\epsilon}}^{\kappa }\right),\end{equation*}$$

and for n > K,

$$\begin{equation*}\mathcal{R}\left(f,{v}_{n}^{a}\right)={\alpha }_{n}^{a}+O\left({{\epsilon}}^{K+4}\right),\end{equation*}$$

where

$$\begin{equation*}O\left({{\epsilon}}^{K+4}\right)=\sum _{j=1}^{m}\sum _{k=K+1}^{n}{\nu }_{j}^{k,a}\left(\genfrac{}{}{0.0pt}{}{n}{k}\right){\left({P}_{j}^{a}\right)}^{n-k}.\end{equation*}$$

From (2.13), we notice that the number of unknowns is now fixed at m(K + 1). Thus, we can increase the number of equations, by taking n large enough, to be able to uniquely determine the unknowns. This is the objective of the next subsection.

Remark 2.1. Assuming that the observation time T is large enough, $T{ >}\frac{\text{diam}\,{\Omega}}{{c}_{\text{s}}}$, then a special solution of (2.6) would be

$$\begin{equation*}\rho \left(z,t\right)={\rho }_{\varepsilon }\left(t+\frac{z}{{c}_{\text{s}}}-\theta \right),\end{equation*}$$

where ${\rho }_{\varepsilon }\in {C}_{0}^{\infty }\left(\mathbb{R}\right)$ is a mollifier function supported in the interval [−ɛ, ɛ] and satisfying ${\int }_{\mathbb{R}}{\rho }_{\varepsilon }\left(s\right)\mathrm{d}s=1$, and θ ∈ [θ_min + ɛ, θ_max − ɛ] with

$$\begin{equation*}{\theta }_{\mathrm{min}}=\frac{{\gamma }_{2}}{{c}_{\text{s}}},\quad {\theta }_{\mathrm{max}}=T+\frac{{\gamma }_{1}}{{c}_{\text{s}}}.\end{equation*}$$

It is obvious that for every θ and ɛ fixed, the function ρ belongs to the functional space C^∞((γ₁, γ₂) × (0, T)). Moreover, for each z fixed, ρ(z, t) is supported in the time interval $\left[\theta -\frac{z}{{c}_{\text{s}}}-\varepsilon ,\theta -\frac{z}{{c}_{\text{s}}}+\varepsilon \right]$. Therefore, taking θ ∈ [θ_min + ɛ, θ_max − ɛ] we get ρ(., T) = ρ_t(., T) = 0.

The mollifier function ρ was used by the authors in [24] for solving the inverse problem of reconstruction of time-varying point sources for wave equation.

Before solving the equation (2.11), we need to know if the projections ${P}_{j}^{a}$ are mutually distinct, which is necessary in order to use the algebraic method proposed below. Indeed, one can remark that there is only a finite number of planes containing the origin such that at least two points among S_j, are projected onto the same point on this plane. So, if a basis is chosen randomly, one is almost sure that the points S_j are projected onto distinct points on every coordinate plane. Therefore, without loss of generality, we assume that:

(H) The projections onto the xy, yz and xz-planes of the points S_j are mutually distinct.

Before presenting our identification algorithm we introduce, for simplicity, some notations and definitions that are used throughout this subsection.

First, denote by ${P}_{j}^{b}$ and ${P}_{j}^{c}$, the projections of S_j onto the yz and xz-complex planes respectively. Then, using in (2.4), the following test functions

$$\begin{equation}{v}_{n}^{b}={\left(y+iz\right)}^{n}{\rho }_{1}\left(x,t\right)\quad \text{and}\quad {v}_{n}^{c}={\left(x+iz\right)}^{n}{\rho }_{2}\left(y,t\right)\end{equation} \tag{ 2.14 }$$

where ρ₁, and ρ₂ are the solutions of (2.6) with the space variable z being replaced by x and y respectively, leads to the following algebraic quantities

$$\begin{equation}{\alpha }_{n}^{b}=\sum _{j=1}^{m}\sum _{k=0}^{K}{\nu }_{j}^{k,b}\left(\genfrac{}{}{0.0pt}{}{n}{k}\right){\left({P}_{j}^{b}\right)}^{n-k}\quad \text{for} \text{all} n\in \mathbb{N}\end{equation} \tag{ 2.15 }$$

and

$$\begin{equation}{\alpha }_{n}^{c}=\sum _{j=1}^{m}\sum _{k=0}^{K}{\nu }_{j}^{k,c}\left(\genfrac{}{}{0.0pt}{}{n}{k}\right){\left({P}_{j}^{c}\right)}^{n-k}\quad \text{for} \text{all}\;n\in \mathbb{N},\end{equation} \tag{ 2.16 }$$

where

$$\begin{equation*}{\nu }_{j}^{k,b}={{\epsilon}}^{3+k}{\int }_{0}^{T}{\int }_{{B}_{j}}{\gamma }_{j}\left(\boldsymbol{\tau },t\right){\rho }_{1}\left({x}_{j}+{\epsilon}{\tau }_{1},t\right){\left[{\tau }_{2}+i{\tau }_{3}\right]}^{k}\mathrm{d}\boldsymbol{\tau }\mathrm{d}t\end{equation*}$$

and

$$\begin{equation*}{\nu }_{j}^{k,c}={{\epsilon}}^{3+k}{\int }_{0}^{T}{\int }_{{B}_{j}}{\gamma }_{j}\left(\boldsymbol{\tau },t\right){\rho }_{2}\left({y}_{j}+{\epsilon}{\tau }_{2},t\right){\left[{\tau }_{1}+i{\tau }_{3}\right]}^{k}\mathrm{d}\boldsymbol{\tau }\mathrm{d}t.\end{equation*}$$

Finally, bringing together the equations (2.13), (2.15) and (2.16), we can write

$$\begin{equation}{\alpha }_{n}^{r}=\sum _{j=1}^{m}\sum _{k=0}^{K}{\nu }_{j}^{k,r}\left(\genfrac{}{}{0.0pt}{}{n}{k}\right){\left({P}_{j}^{r}\right)}^{n-k},\quad \text{for} r=a,b,c \text{and} n\in \mathbb{N}.\end{equation} \tag{ 2.17 }$$

On the other hand, we assume that we know an upper bound M of the number of sources m, then $\overline{J}=M\left(K+1\right)$ is an upper bound of the number J = m(K + 1). We define for r = a, b, c, and k = 0, ..., K, the complex vectors

$$\begin{equation*}{\xi }_{n}^{r}={\left({\alpha }_{n}^{r},\dots ,\;{\alpha }_{\overline{J}+n-1}^{r}\right)}^{t},\quad {{\Lambda}}_{k}^{r}={\left({\nu }_{1}^{k,r},\dots ,\;{\nu }_{m}^{k,r}\right)}^{t}\quad \;\text{and}\;\;{{\Lambda}}^{r}={\left({\left({{\Lambda}}_{0}^{r}\right)}^{t},{\left({{\Lambda}}_{1}^{r}\right)}^{t},\dots ,\;{\left({{\Lambda}}_{K}^{r}\right)}^{t}\right)}^{t},\end{equation*}$$

of sizes $\overline{J},m$ and m(K + 1) respectively, and consider, for all $n\in \mathbb{N}$, the complex matrices ${A}_{n}^{r}$ of size $\overline{J}{\times}J$,

$$\begin{equation}{A}_{n}^{r}=\left({V}_{n}^{0,r},\dots ,\;{V}_{n}^{K,r}\right),\end{equation} \tag{ 2.18 }$$

where, for k ∈ {0, ..., K}, ${V}_{n}^{k,r}$ are the confluent $\overline{J}{\times}m$ Vandermonde matrices

$$\begin{equation*}{V}_{n}^{k,r}=\left(\begin{array}{ccc}\left(\genfrac{}{}{0.0pt}{}{n}{k}\right){\left({P}_{1}^{r}\right)}^{n-k}\hfill & \hfill \cdots \hfill & \left(\genfrac{}{}{0.0pt}{}{n}{k}\right){\left({P}_{m}^{r}\right)}^{n-k}\hfill \\ \left(\genfrac{}{}{0.0pt}{}{n+1}{k}\right){\left({P}_{1}^{r}\right)}^{n-k+1}\hfill & \hfill \cdots \hfill & \left(\genfrac{}{}{0.0pt}{}{n+1}{k}\right){\left({P}_{m}^{r}\right)}^{n-k+1}\hfill \\ {\vdots}\hfill & \hfill \ddots \hfill & {\vdots}\hfill \\ \left(\genfrac{}{}{0.0pt}{}{\overline{J}+n-1}{k}\right){\left({P}_{1}^{r}\right)}^{n-k+\overline{J}-1}\hfill & \hfill \cdots \hfill & \left(\genfrac{}{}{0.0pt}{}{\overline{J}+n-1}{k}\right){\left({P}_{m}^{r}\right)}^{n-k+\overline{J}-1}\hfill \end{array}\right).\end{equation*}$$

Now, let us introduce the Hankel matrix

$$\begin{equation}{H}_{\overline{J}}^{r}=\left(\begin{pmatrix}{cccc}\hfill {\alpha }_{0}^{r}\hfill & \hfill {\alpha }_{1}^{r}\hfill & \hfill \cdots \hfill & \hfill {\alpha }_{\overline{J}-1}^{r}\hfill \\ \hfill {\alpha }_{1}^{r}\hfill & \hfill {\alpha }_{2}^{r}\hfill & \hfill \cdots \hfill & \hfill {\alpha }_{\overline{J}}^{r}\hfill \\ \hfill {\vdots}\hfill & \hfill {\vdots}\hfill & \hfill {\vdots}\hfill & \hfill {\vdots}\hfill \\ \hfill {\alpha }_{\overline{J}-1}^{r}\hfill & \hfill {\alpha }_{\overline{J}}^{r}\hfill & \hfill \cdots \hfill & \hfill {\alpha }_{2\overline{J}-2}^{r}\hfill \end{pmatrix}\right)\end{equation} \tag{ 2.19 }$$

and the following multi-diagonal matrix

$$\begin{equation}{L}^{r}=\left(\begin{pmatrix}{cccc}\hfill {\nu }^{0,r}\hfill & \hfill {\nu }^{1,r}\hfill & \hfill \cdots \hfill & \hfill {\nu }^{K,r}\hfill \\ \hfill {\vdots}\hfill & \hfill {\vdots}\hfill & \hfill ...\hfill & \hfill {\vdots}\hfill \\ \hfill {\nu }^{K-1,r}\hfill & \hfill {\nu }^{K,r}\hfill & \hfill \cdots \hfill & \hfill 0\hfill \\ \hfill {\nu }^{K,r}\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill \end{pmatrix}\right),\end{equation} \tag{ 2.20 }$$

where

$$\begin{equation*}{\nu }^{k,r}=\mathrm{diag}\left({\nu }_{1}^{k,r},\dots ,\;{\nu }_{m}^{k,r}\right).\end{equation*}$$

We propose an identification processes in four steps. First, we determine the absorbers' number m using the Hankel matrix defined in (2.19), then we obtain the projections ${P}_{j}^{r}$ using a companion matrix. Next, we get the coefficients ${\nu }_{j}^{k,r}$ as a solution of a linear system via the Hankel matrix ${H}_{J}^{r}$, and finally we end up with the locations S_j by solving a similar linear system and using the coefficients ${\nu }_{j}^{k,r}$.

Step 1: determination of the number of absorbers

The first step consists in determining the number of sources by means of the following theorem.

Theorem 2.2. Let ${H}_{\overline{J}}^{r}$ be the Hankel matrix defined in (2.19) where $\overline{J}$ is a known upper bound of J = m(K + 1). Assume that, for r = a, b, c, the projected points ${P}_{j}^{r}$ of S_j are mutually distinct, then we have

$$\begin{equation*}\text{rank}\left({H}_{\overline{J}}^{r}\right)=m\left(K+1\right)\quad \text{if}\;\text{and} \text{only} \text{if} {\nu }_{j}^{K,r}\ne 0\quad \text{for} j=1,\dots ,\;m.\end{equation*}$$

The proof of this theorem is similar to that of theorem 2.8 in [19]. It is based on the decomposition of the Hankel matrix ${H}_{\overline{J}}^{r}$ given in lemma 2.9 in [19]. For the convenience of the reader, we recall this lemma here for which we present a sketch of the proof.

Lemma 2.3. Let ${H}_{\overline{J}}^{r}$ be the Hankel matrix defined in (2.19), L^r the multi-diagonal matrix defined in (2.20), ${A}_{0}^{r}$ the Vandermonde matrix defined in (2.18), and ${\left({A}_{0}^{r}\right)}^{t}$ its transpose matrix. Then,

$$\begin{equation}{H}_{\overline{J}}^{r}={A}_{0}^{r}{L}^{r}{\left({A}_{0}^{r}\right)}^{t}.\end{equation} \tag{ 2.21 }$$

Proof. Indeed, from (2.18), we notice that

$$\begin{equation*}{A}_{n}^{r}{{\Lambda}}^{r}=\sum _{k=0}^{K}{V}_{n}^{k,r}{{\Lambda}}_{k}^{r}.\end{equation*}$$

Moreover,

$$\begin{equation*}{V}_{n}^{k,r}{{\Lambda}}_{k}^{r}=\left(\begin{pmatrix}{c}\hfill \sum _{j=1}^{m}\left(\genfrac{}{}{0.0pt}{}{n}{k}\right){\left({P}_{j}^{r}\right)}^{n-k}{\nu }_{j}^{k,r}\hfill \\ \hfill {\vdots}\hfill \\ \hfill \sum _{j=1}^{m}\left(\genfrac{}{}{0.0pt}{}{\overline{J}+n-1}{k}\right){\left({P}_{j}^{r}\right)}^{n-k+\overline{J}-1}{\nu }_{j}^{k,r}\hfill \end{pmatrix}\right).\end{equation*}$$

Therefore, the algebraic relation (2.17) can be rewritten in a matrix form as

$$\begin{equation}{\xi }_{n}^{r}={A}_{n}^{r}{{\Lambda}}^{r},\quad \forall n\in \mathbb{N}.\end{equation} \tag{ 2.22 }$$

Furthermore, if we denote by T^r the block J × J upper triangular complex matrix

$$\begin{equation}{T}^{r}=\left(\begin{pmatrix}{cccc}\hfill {D}_{p}^{r}\hfill & \hfill I\hfill & \hfill \cdots \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill {\vdots}\hfill \\ \hfill {\vdots}\hfill & \hfill \ddots \hfill & \hfill {D}_{p}^{r}\hfill & \hfill I\hfill \\ \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill & \hfill {D}_{p}^{r}\hfill \end{pmatrix}\right),\end{equation} \tag{ 2.23 }$$

with

$$\begin{equation*}{D}_{p}^{r}=\mathrm{diag}\left({P}_{1}^{r},\dots ,\;{P}_{m}^{r}\right)\quad \text{and}\quad I=\mathrm{diag}\left(1,\dots ,\;1\right),\end{equation*}$$

then, using the Pascal formula $\left(\genfrac{}{}{0.0pt}{}{n}{j-1}\right)+\left(\genfrac{}{}{0.0pt}{}{n}{j}\right)=\left(\genfrac{}{}{0.0pt}{}{n+1}{j}\right)$, it follows that

$$\begin{equation*}{V}_{n+1}^{k,r}={V}_{n}^{k-1,r}+{V}_{n}^{k,r}{D}_{p}^{r},\quad \forall n\in \mathbb{N},\quad k=1,\dots ,\;K\end{equation*}$$

and

$$\begin{equation*}{V}_{n+1}^{0,r}={V}_{n}^{0,r}{D}_{p}^{r},\quad \forall n\in \mathbb{N}.\end{equation*}$$

From the definition of the matrix ${A}_{n}^{r}$, we obtain

$$\begin{equation*}{A}_{n+1}^{r}={A}_{n}^{r}{T}^{r}={A}_{0}^{r}{\left({T}^{r}\right)}^{n+1},\quad \forall n\in \mathbb{N},\end{equation*}$$

and therefore, from (2.22), one gets

$$\begin{equation}{\xi }_{n}^{r}={A}_{0}^{r}{\left({T}^{r}\right)}^{n}{{\Lambda}}^{r},\quad \forall n\in \mathbb{N}.\end{equation} \tag{ 2.24 }$$

Now, thanks to (2.24), one can verify by a simple calculation the following relationship

$$\begin{equation*}{H}_{\overline{J}}^{r}={A}_{0}^{r}\left[{{\Lambda}}^{r},{T}^{r}{{\Lambda}}^{r},\dots ,\;{\left({T}^{r}\right)}^{\overline{J}-1}{{\Lambda}}^{r}\right]={A}_{0}^{r}{L}^{r}{\left({A}_{0}^{r}\right)}^{t},\end{equation*}$$

which ends the proof of Lemma 2.3. □

Proof of Theorem 2.2. We can easily check that $\text{rank}\left({A}_{0}^{r}\right)=\text{rank}{\left({A}_{0}^{r}\right)}^{t}=J$ if and only if the projections ${P}_{j}^{r}$ are mutually distinct, for j = 1, ..., m. Moreover, the matrix L^r is nonsingular if and only if ${\nu }_{j}^{K,r}\ne 0\quad \forall j=1,\dots ,\;m$. Therefore, under the hypothesis (H), and assuming that the coefficients ${\nu }_{j}^{K,r}\ne 0$, we deduce that ${A}_{0}^{r}{L}^{r}$ is surjective and therefore we have $\text{rank}\left({A}_{0}^{r}{L}^{r}{\left({A}_{0}^{r}\right)}^{t}\right)=\text{rank}{\left({A}_{0}^{r}\right)}^{t}=J$, which leads to the desired result. □

Step 2: reconstruction of the projections ${P}_{j}^{r}$

The second step consists in determining the projections onto the xy-, xz- and yz-complex planes of the centers S_j.

Henceforth, we replace $\overline{J}$ by J in the quantities defined above. Thus, from (2.24), we can easily derive that

$$\begin{equation*}{\xi }_{n+1}^{r}={B}^{r}{\xi }_{n}^{r},\quad \forall n\in \mathbb{N},\end{equation*}$$

where we have set

$$\begin{equation}{B}^{r}={A}_{0}^{r}{T}^{r}{\left({A}_{0}^{r}\right)}^{-1}.\end{equation} \tag{ 2.25 }$$

Here, the matrix ${A}_{0}^{r}$ is invertible since the projected points ${P}_{j}^{r}$ are assumed distinct. Moreover, since $\text{rank}\left({H}_{J}^{r}\right)=J$, the family ${\left({\xi }_{n}^{r}\right)}_{n=0,\dots ,J-1}$ forms a basis of ${\mathbb{C}}^{J}$, so the J × J complex matrix B^r is given explicitly by

$$\begin{equation}{B}^{r}=\left(\begin{pmatrix}{ccccc}\hfill 0\hfill & \hfill 1\hfill & \hfill \cdots \hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill \cdots \hfill & \hfill 0\hfill \\ \hfill {\vdots}\hfill & \hfill {\vdots}\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill 0\hfill \\ \hfill {\vdots}\hfill & \hfill {\vdots}\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill 1\hfill \\ \hfill {c}_{0}^{r}\hfill & \hfill {c}_{1}^{r}\hfill & \hfill \cdots \hfill & \hfill \cdots \hfill & \hfill {c}_{J-1}^{r}\hfill \end{pmatrix}\right),\end{equation} \tag{ 2.26 }$$

where the vector ${C}^{r}={\left({c}_{0}^{r},\dots ,\;{c}_{J-1}^{r}\right)}^{t}$ is obtained by solving the linear system

$$\begin{equation}{H}_{J}^{r}{C}^{r}={\xi }_{J}^{r}.\end{equation} \tag{ 2.27 }$$

Indeed, as the family ${\left({\xi }_{n}^{r}\right)}_{n=0,\dots ,J-1}$ forms a basis of ${\mathbb{C}}^{J}$, given ξ_J, there exists a vector ${C}^{r}={\left({c}_{0}^{r},\dots ,\;{c}_{J-1}^{r}\right)}^{t}$ such that

$$\begin{equation*}{\xi }_{J}^{r}=\sum _{n=0}^{J-1}{c}_{n}^{r}{\xi }_{n}^{r},\end{equation*}$$

which is equivalent to (2.27) since the matrix ${H}_{J}^{r}=\left({\xi }_{0}^{r},\dots ,\;{\xi }_{J-1}^{r}\right)$.

Finally, the projections ${P}_{j}^{r}$ are given by the following theorem.

Theorem 2.4. Let B^r, for r = a, b, c, be the companion matrices defined in (2.26). Assume that the projected points ${P}_{j}^{r}$ of S_j are distinct and away from the origin and that ${\nu }_{j}^{K,r}\ne 0\;$ for j = 1, ..., m. Then,

• (a)  
  B^r admits m eigenvalues of multiplicity K + 1;
• (b)  
  The m eigenvalues of multiplicity K + 1 are the projections ${P}_{j}^{r}$ of the point sources S_j.

Proof. The proof of this theorem follows from (2.23) and (2.25). □

Remark 2.5. We know according to formulas (2.23), (2.25) and hypothesis (H), that the matrix B^r admits m distinct eigenvalues having the same multiplicity K + 1, where K is given. Then, numerically, one should obtain m distinct eigenvalues, m being the rank of the Hankel matrix ${H}_{\overline{J}}^{r}$ defined in (2.19). The rank of ${H}_{\overline{J}}^{r}$ can be obtained numerically as the number of its non-zero singular values.

Step 3: determination of the vectors Λ^r

In this step, we shall determine the quantities ${\nu }_{j}^{k,r}$, which we will use later to reconstruct the locations S_j. In fact, in order to determine ${\nu }_{j}^{k,r}$, it is sufficient to solve the linear systems

$$\begin{equation}{A}_{0}^{r}{{\Lambda}}^{r}={\xi }_{0}^{r}.\end{equation} \tag{ 2.28 }$$

Step 4: complete reconstruction of locations S_j

In step 2, we were able to reconstruct modulo ^K+4 the projections of the centers S_j on the xy-, yz- and xz-planes. However, this is not enough. In fact, the combination of different projections does not uniquely determine the centers S_j, see remark 2.6. For this purpose, we fix r = a and we reconstruct, modulo ⁴, the third components z_j that correspond to the previously calculated projections ${P}_{j}^{a}$, on which we need that the coefficients ${\nu }_{j}^{0,a}$ calculated in the previous step to satisfy

$$\begin{equation}{\nu }_{j}^{0,a}\ne 0\;\quad \text{for} \;j=1,\dots ,\;m.\end{equation} \tag{ 2.29 }$$

We propose an algorithm for the reconstruction of the components z_j using specific wave functions. For this purpose, for $n\in \mathbb{N}$, we take ${v}_{n}^{a}$ defined in (2.7) and we consider the following wave functions,

$$\begin{equation}{g}_{n}^{a}\left(\mathbf{\text{x}},t\right)=-\frac{\partial }{\partial z}{v}_{n}^{a}\left(\mathbf{\text{x}},t\right)=-{\rho }_{z}\left(z,t\right){\left(x+iy\right)}^{n},\end{equation} \tag{ 2.30 }$$

$$\begin{align}\hfill {h}_{n}^{a}\left(\mathbf{\text{x}},t\right)& =z\left(\frac{\partial }{\partial x}-i\frac{\partial }{\partial y}\right){v}_{n}^{a}\left(\mathbf{\text{x}},t\right)-\left(x-iy\right)\frac{\partial }{\partial z}{v}_{n}^{a}\left(\mathbf{\text{x}},t\right)\hfill \\ \hfill & =2nz\rho \left(z,t\right){\left(x+iy\right)}^{n-1}-{\rho }_{z}\left(z,t\right)\left(x-iy\right){\left(x+iy\right)}^{n},\hfill \end{align} \tag{ 2.31 }$$

which are solutions of (2.1) satisfying the condition (2.5). Therefore, substituting ${g}_{n}^{a}$ in (2.4) we get

$$\begin{align*}\hfill R\left(f,{g}_{n}^{a}\right)& =\sum _{j=1}^{m}{\int }_{0}^{T}{\int }_{{D}_{j}}{\sigma }_{j}\left(\mathbf{\text{x}}\right)g\left(\mathbf{\text{x}},t\right){g}_{n}^{a}\left(\mathbf{\text{x}},t\right)\mathrm{d}\mathbf{\text{x}}\mathrm{d}t+O\left({{\epsilon}}^{\kappa }\right)\hfill \\ \hfill & =-\sum _{j=1}^{m}{\int }_{0}^{T}{\int }_{{D}_{j}}{\sigma }_{j}\left(\mathbf{\text{x}}\right)g\left(\mathbf{\text{x}},t\right){\rho }_{z}\left(z,t\right){\left(x+iy\right)}^{n}\mathrm{d}\mathbf{\text{x}}\mathrm{d}t+O\left({{\epsilon}}^{\kappa }\right)\hfill \\ \hfill & =\sum _{j=1}^{m}\sum _{k=1}^{n}\left(\genfrac{}{}{0.0pt}{}{n}{k}\right){\tilde {{\nu }_{j}}}^{k,a}{\left({P}_{j}^{a}\right)}^{n-k}+O\left({{\epsilon}}^{\kappa }\right),\hfill \end{align*}$$

where for every j = 1, ..., m and k = 1, ..., K,

$$\begin{equation}{\tilde {{\nu }_{j}}}^{k,a}=-{{\epsilon}}^{3+k}{\int }_{0}^{T}{\int }_{{B}_{j}}{\sigma }_{j}\left({\mathbf{S}}_{j}+{\epsilon}\boldsymbol{\tau }\right)g\left({\mathbf{S}}_{j}+{\epsilon}\boldsymbol{\tau },t\right){\rho }_{z}\left({z}_{j}+{\epsilon}{\tau }_{3},t\right){\left({\tau }_{1}+i{\tau }_{2}\right)}^{k}\mathrm{d}\boldsymbol{\tau }\mathrm{d}t.\end{equation} \tag{ 2.32 }$$

Similar to the previous section, we have for n ⩽ K

$$\begin{equation}R\left(f,{g}_{n}^{a}\right)=\sum _{j=1}^{m}\sum _{k=0}^{K}\left(\genfrac{}{}{0.0pt}{}{n}{k}\right){\tilde {{\nu }_{j}}}^{k,a}{\left({P}_{j}^{a}\right)}^{n-k}+O\left({{\epsilon}}^{\kappa }\right),\end{equation} \tag{ 2.33 }$$

and for n > K

$$\begin{equation}R\left(f,{g}_{n}^{a}\right)=\sum _{j=1}^{m}\sum _{k=0}^{K}\left(\genfrac{}{}{0.0pt}{}{n}{k}\right){\tilde {{\nu }_{j}}}^{k,a}{\left({P}_{j}^{a}\right)}^{n-k}+O\left({{\epsilon}}^{K+4}\right).\end{equation} \tag{ 2.34 }$$

Now, if we choose for $n\in \mathbb{N}$

$$\begin{equation}{\tilde {\alpha }}_{n}^{a}=\sum _{j=1}^{m}\sum _{k=0}^{K}\left(\genfrac{}{}{0.0pt}{}{n}{k}\right){\tilde {{\nu }_{j}}}^{k,a}{\left({P}_{j}^{a}\right)}^{n-k},\end{equation} \tag{ 2.35 }$$

and we define the following quantities

$$\begin{equation*}{\tilde {\xi }}_{n}^{a}={\left({\tilde {\alpha }}_{n}^{a},\dots ,\;{\tilde {\alpha }}_{J+n-1}\right)}^{t},\end{equation*}$$

$$\begin{equation*}{\tilde {{\Lambda}}}_{k}={\left({\tilde {\nu }}_{1}^{k,a},\dots ,\;{\tilde {\nu }}_{m}^{k,a}\right)}^{t}\quad \text{for} \;k=0,\dots ,\;K\end{equation*}$$

and

$$\begin{equation*}\tilde {{\Lambda}}={\left({\left({\tilde {{\Lambda}}}_{0}\right)}^{t},\dots ,\;{\left({\tilde {{\Lambda}}}_{K}\right)}^{t}\right)}^{t},\end{equation*}$$

where J = m(K + 1), we get using the formula (2.35) the following system of equations

$$\begin{equation}{A}_{0}^{a}\tilde {{\Lambda}}={\tilde {\xi }}_{0}^{a},\end{equation} \tag{ 2.36 }$$

where ${A}_{0}^{a}$ is defined in (2.18), and ${\tilde {\xi }}_{0}^{a}$ can be approximated using the equations (2.34) and (2.35). Under the assumption that the projections ${P}_{j}^{a}$ are mutually distinct, then ${A}_{0}^{a}$ is invertible, and system (2.36) admits a unique solution. By solving this system, we can reconstruct the quantities ${\tilde {{\nu }_{j}}}^{k,a}$.On the other hand, if we substitute ${h}_{n}^{a}$ in (2.4) we get

$$\begin{align*}\hfill R\left(f,{h}_{n}^{a}\right)& =\sum _{j=1}^{m}{\int }_{0}^{T}{\int }_{{D}_{j}}{\sigma }_{j}\left(\mathbf{\text{x}}\right)g\left(\mathbf{\text{x}},t\right){h}_{n}^{a}\left(\mathbf{\text{x}},t\right)\mathrm{d}\mathbf{\text{x}}\mathrm{d}t+O\left({{\epsilon}}^{\kappa }\right)\hfill \\ \hfill & =2n\sum _{j=1}^{m}{\int }_{0}^{T}{\int }_{{D}_{j}}{\sigma }_{j}\left(\mathbf{\text{x}}\right)g\left(\mathbf{\text{x}},t\right)\rho \left(z,t\right)z{\left(x+iy\right)}^{n-1}\mathrm{d}\mathbf{\text{x}}\mathrm{d}t\hfill \\ \hfill & -\sum _{j=1}^{m}{\int }_{0}^{T}{\int }_{{D}_{j}}{\sigma }_{j}\left(\mathbf{\text{x}}\right)g\left(\mathbf{\text{x}},t\right){\rho }_{z}\left(z,t\right)\left(x-iy\right){\left(x+iy\right)}^{n}\mathrm{d}\mathbf{\text{x}}\mathrm{d}t+O\left({{\epsilon}}^{\kappa }\right).\hfill \end{align*}$$

Applying the change of variables D_j = S_j + B_j one gets,

$$\begin{align*}\hfill R\left(f,{h}_{n}^{a}\right)& =2n\sum _{j=1}^{m}{{\epsilon}}^{3}{\int }_{0}^{T}{\int }_{{B}_{j}}{\gamma }_{j}\left(\boldsymbol{\tau },t\right)\rho \left({z}_{j}+{\epsilon}{\tau }_{3},t\right)\left({z}_{j}+{\epsilon}{\tau }_{3}\right){\left({P}_{j}^{a}+{\epsilon}\left({\tau }_{1}+i{\tau }_{2}\right)\right)}^{n-1}\mathrm{d}\boldsymbol{\tau }\;\mathrm{d}t\hfill \\ \hfill & -\;\sum _{j=1}^{m}{{\epsilon}}^{3}{\int }_{0}^{T}{\int }_{{B}_{j}}{\gamma }_{j}\left(\boldsymbol{\tau },t\right){\rho }_{z}\left({z}_{j}+{\epsilon}{\tau }_{3},t\right)\left({x}_{j}-i{y}_{j}+{\epsilon}\left({\tau }_{1}-i{\tau }_{2}\right)\right)\hfill \\ \hfill & {\times}{\left({P}_{j}^{a}+{\epsilon}\left({\tau }_{1}+i{\tau }_{2}\right)\right)}^{n}\mathrm{d}\boldsymbol{\tau }\mathrm{d}t+O\left({{\epsilon}}^{\kappa }\right)\hfill \\ \hfill & =2n\sum _{j=1}^{m}{{\epsilon}}^{3}{\int }_{0}^{T}{\int }_{{B}_{j}}{\gamma }_{j}\left(\boldsymbol{\tau },t\right)\rho \left({z}_{j}+{\epsilon}{\tau }_{3},t\right){z}_{j}{\left({P}_{j}^{a}\right)}^{n-1}\mathrm{d}\boldsymbol{\tau }\mathrm{d}t\hfill \\ \hfill & -\sum _{j=1}^{m}{{\epsilon}}^{3}{\int }_{0}^{T}{\int }_{{B}_{j}}{\gamma }_{j}\left(\boldsymbol{\tau },t\right){\rho }_{z}\left({z}_{j}+{\epsilon}{\tau }_{3},t\right)\left({x}_{j}-i{y}_{j}\right){\left({P}_{j}^{a}\right)}^{n}\mathrm{d}\boldsymbol{\tau }\mathrm{d}t+O\left({{\epsilon}}^{4}\right),\hfill \end{align*}$$

where γ_j(τ, t) is defined in (2.10). Therefore, one has

$$\begin{equation*}R\left(f,{h}_{n}^{a}\right)=2n\sum _{j=1}^{m}{\left({P}_{j}^{a}\right)}^{n-1}{\nu }_{j}^{0,a}{z}_{j}-\sum _{j=1}^{m}\left({x}_{j}-i{y}_{j}\right){\left({P}_{j}^{a}\right)}^{n}{\tilde {{\nu }_{j}}}^{0,a}+O\left({{\epsilon}}^{4}\right).\end{equation*}$$

Hence,

$$\begin{equation*}\frac{1}{2n}\left(R\left(f,{h}_{n}^{a}\right)+\sum _{j=1}^{m}\left({x}_{j}-i{y}_{j}\right){\left({P}_{j}^{a}\right)}^{n}{\tilde {{\nu }_{j}}}^{0,a}\right)=\sum _{j=1}^{m}{\left({P}_{j}^{a}\right)}^{n-1}{\nu }_{j}^{0,a}{z}_{j}+O\left({{\epsilon}}^{4}\right)\quad \forall n\in \mathbb{N}.\end{equation*}$$

Now, if we define for $n\in \mathbb{N}$

$$\begin{equation}{\tilde {R}}_{n}=\frac{1}{2n}\left(R\left(f,{h}_{n}^{a}\right)+\sum _{j=1}^{m}\left({x}_{j}-i{y}_{j}\right){\left({P}_{j}^{a}\right)}^{n}{\tilde {{\nu }_{j}}}^{0,a}\right),\end{equation} \tag{ 2.37 }$$

and we take

$$\begin{equation*}\tilde {R}={\left({\tilde {R}}_{1},\dots ,\;{\tilde {R}}_{m}\right)}^{t}\end{equation*}$$

and

$$\begin{equation*}{\tilde {{\Lambda}}}_{z}={\left({\nu }_{1}^{0,a}{z}_{1},\dots ,\;{\nu }_{m}^{0,a}{z}_{m}\right)}^{t},\end{equation*}$$

we get, modulo O(⁴), the following system

$$\begin{equation}{P}_{m}{\tilde {{\Lambda}}}_{z}=\tilde {R},\end{equation} \tag{ 2.38 }$$

where P_m is the m × m Vandermonde matrix defined by

$$\begin{equation*}{P}_{m}=\left(\begin{pmatrix}{cccc}\hfill 1\hfill & \hfill 1\hfill & \hfill \dots \hfill & \hfill 1\hfill \\ \hfill {P}_{1}^{a}\hfill & \hfill {P}_{2}^{a}\hfill & \hfill \dots \hfill & \hfill {P}_{m}^{a}\hfill \\ \hfill {\vdots}\hfill & \hfill \hfill & \hfill \hfill & \hfill {\vdots}\hfill \\ \hfill {\left({P}_{1}^{a}\right)}^{m-1}\hfill & \hfill {\left({P}_{2}^{a}\right)}^{m-1}\hfill & \hfill \dots \hfill & \hfill {\left({P}_{m}^{a}\right)}^{m-1}\hfill \end{pmatrix}\right).\end{equation*}$$

Since the projections ${P}_{j}^{a}$ are assumed to be mutually distinct, thus P_m is invertible and equation (2.38) admits a unique solution. By solving this equation we can reconstruct ${\nu }_{j}^{0,a}{z}_{j}$ up to O(⁴), and since ${\nu }_{j}^{0,a}$ are previously calculated and chosen to be different from zero, thus we can calculate the values of z_j. Finally, combining this result with that of section 2.3, we get the locations of the centers S_j.

Remark 2.6. 

• (a)  
  The projections on the xy-, yz-, and xz-planes do not uniquely determine the centers S_j. For example, the two sets of points

  $$\begin{equation*}{\mathcal{F}}_{1}=\left\{\left(1,0,2\right),\left(1,-1,3\right),\left(2,0,3\right),\left(2,-1,2\right)\right\},\end{equation*}$$

  and

  $$\begin{equation*}{\mathcal{F}}_{2}=\left\{\left(1,0,3\right),\left(1,-1,2\right),\left(2,0,2\right),\left(2,-1,3\right)\right\}\end{equation*}$$

  have the same projections

  $$\begin{equation*}{P}^{a}=\left\{\left(1,0\right),\left(1,-1\right),\left(2,0\right),\left(2,-1\right)\right\},\end{equation*}$$

  $$\begin{equation*}{P}^{b}=\left\{\left(0,2\right),\left(-1,3\right),\left(0,3\right),\left(-1,2\right)\right\},\end{equation*}$$

  and

  $$\begin{equation*}{P}^{c}=\left\{\left(1,2\right),\left(1,3\right),\left(2,3\right),\left(2,2\right)\right\}\end{equation*}$$

  on the xy-, yz- and xz-planes respectively. Thus, knowing the projections P^a, P^b and P^c is not enough to uniquely determine the corresponding locations S_j.
• (b)  
  The last step allows us to uniquely determine the centers S_j. However, we lose some accuracy as the obtained values of z_j are given modulo ⁴.

Summary: Identification algorithm

To sum up, our identification algorithm can be summarized as follows. First, we choose a fixed integer K > 0 such that ^K+4 is small enough and we suppose that we know an upper bound $\overline{J}$ of the number J = m(K + 1). Then, we proceed in the following steps:

• Step 1. We fix r = a, b or c, and we estimate for $n=0,\dots ,\;2\overline{J}-2$ the coefficients ${\alpha }_{n}^{r}$ by $R\left(f,{v}_{n}^{r}\right)$ defined in (2.3), which introduces an accuracy error O(^K+4) in our identification algorithm. Then we construct the Hankel matrix ${H}_{\overline{J}}^{r}$ defined in (2.19).
• Step 2. We calculate J as the rank of the Hankel matrix ${H}_{\overline{J}}^{r}$, which can be numerically estimated using, for example, the singular value decomposition (SVD) method, and we deduce the number of absorbers m using the relation J = m(K + 1).
• Step 3. We construct ${C}^{r}={\left({c}_{0}^{r},\dots ,\;{c}_{J-1}^{r}\right)}^{t}$ by solving the system ${H}_{J}^{r}{C}^{r}={\xi }_{J}^{r}$, where ${\xi }_{J}^{r}={\left({\alpha }_{J}^{r},\dots ,\;{\alpha }_{2J-1}^{r}\right)}^{t}$. The projections ${P}_{j}^{r}$ of the centers S_j then can be determined as the eigenvalues of the matrix B^r defined in (2.26), which can be constructed using the previously calculated coefficients ${\left\{{c}_{j}^{r}\right\}}_{0{\leqslant}j{\leqslant}J-1}$.
• Step 4. The coefficients ${\nu }_{j}^{k,a}$ can be obtained by solving the system (2.28).
• Step 5. We fix r = a and we estimate for n = 0, 1, ..., J − 1 the coefficients ${\tilde {\alpha }}_{n}^{a}$ by $R\left(f,{g}_{n}^{a}\right)$ defined in (2.3), where ${g}_{n}^{a}$ is given by (2.30), then we calculate the coefficients ${\tilde {\nu }}_{j}^{k,a}$ defined in (2.32) by solving system (2.36).
• Step 6. We calculate $R\left(f,{h}_{n}^{a}\right)$ for n = 1, ..., m, where ${h}_{n}^{a}$ is defined by (2.31). Then using these values and the previously calculated coefficients ${\tilde {\nu }}_{j}^{k,a}$, we may calculate the coefficients ${\tilde {R}}_{n}$ defined in (2.37).
• Step 7. Finally, by solving (2.38) we get the values of ${\nu }_{j}^{0,a}{z}_{j}$ for j = 1, ..., m. Then under the assumption that the coefficients ${\nu }_{j}^{0,a}$ calculated in step 4 are different from zero, we may calculate the values of z_j that correspond to the projections ${P}_{j}^{a}$, and thus we deduce the complete location of the centers S_j.

Remark 2.7. From the previous calculations, one can see that the resolution of the problem depends on the conductivity of the background medium σ₀. More precisely, on how much σ₀ is small. In other words, from equations (2.11)–(2.14), we notice that to solve the problem, it is necessary to have κ > 4 + K, for some positive integer K, where κ is given by the relation σ₀ = O(^κ). On the other hand, the conductivity coefficient σ₀ plays an important role in the resolution of the problem in higher dimensions, on which we notice that our reconstruction algorithm can be generalized to higher dimensions ${\mathbb{R}}^{N}$, if we can find a positive integer K such that

$$\begin{equation*}\quad N+1+K{< }\kappa .\end{equation*}$$

Detailed explanations are given in the appendix.

Before giving our stability results, we first introduce several notations and definitions needed throughout the rest of the paper.

First, let S = (x, y, z) ∈ Ω and d(Γ, S) be the Euclidean distance between the boundary Γ and S. Let D_j, j = 1, ..., m, be m absorbers with centers

$$\begin{equation*}{\mathbf{S}}_{j}=\left({x}_{j},{y}_{j},{z}_{j}\right),\quad \text{for} j=1,\dots ,\;m.\end{equation*}$$

We set

$$\begin{equation*}\alpha =\underset{1{\leqslant}j{\leqslant}m}{\mathrm{min}}d\left({\Gamma},{\mathbf{S}}_{j}\right),\end{equation*}$$

which is strictly positive, since S_j ∈ Ω. Then, we define the set

$$\begin{equation*}{{\Omega}}_{\alpha }=\left\{S\in {\Omega}:d\left({\Gamma},S\right){\geqslant}\alpha \right\},\end{equation*}$$

and let

$$\begin{equation*}d=\text{diam}\left({\Omega}\right)-\alpha .\end{equation*}$$

As denoted before, we take

$$\begin{equation*}{P}_{j}^{a}={x}_{j}+i{y}_{j}\quad \text{and}\;\quad {P}_{j}^{b}={y}_{j}+i{z}_{j}\end{equation*}$$

to be the projections of the centers S_j onto the xy and yz-planes respectively, and set

$$\begin{equation*}{P}^{a}={\left({P}_{j}^{a}\right)}_{1{\leqslant}j{\leqslant}m}\quad \text{and}\;\quad {P}^{b}={\left({P}_{j}^{b}\right)}_{1{\leqslant}j{\leqslant}m}.\end{equation*}$$

Then, we introduce the following real coefficients

$$\begin{equation*}{\delta }_{r}=\underset{\underset{j\ne n}{1{\leqslant}j{\leqslant}m,\;\;1{\leqslant}n{\leqslant}m}}{\mathrm{min}}{\Vert}{P}_{j}^{r}-{P}_{n}^{r}{\Vert},\quad r=a,b,\end{equation*}$$

and take

$$\begin{equation}\delta =\mathrm{min}\left({\delta }_{a},{\delta }_{b}\right).\end{equation} \tag{ 3.1 }$$

Thus, δ is strictly positive according to hypothesis (H). Finally, for two points configurations ${R}^{\ell }={\left({R}_{j}^{\ell }\right)}_{1{\leqslant}j{\leqslant}{N}_{\ell }}$ for ℓ = 1, 2, we recall the Hausdorff distance between R¹ and R², defined as follows

$$\begin{equation*}{d}_{H}\left({R}^{1},{R}^{2}\right)=\mathrm{max}\left[\underset{1{\leqslant}n{\leqslant}{N}_{2};}{\mathrm{max}}\underset{1{\leqslant}j{\leqslant}{N}_{1}}{\mathrm{min}}{\Vert}{R}_{n}^{2}-{R}_{j}^{1}{\Vert},\underset{1{\leqslant}n{\leqslant}{N}_{1};}{\mathrm{max}}\underset{1{\leqslant}j{\leqslant}{N}_{2}}{\mathrm{min}}{\Vert}{R}_{n}^{1}-{R}_{j}^{2}{\Vert}\right].\end{equation*}$$

Let

$$\begin{equation*}{\left({\gamma }_{j}^{\ell },{\mathbf{S}}_{j}^{\ell }\right)}_{1{\leqslant}j{\leqslant}{m}_{\ell }}\quad \text{for} \ell =1,2,\end{equation*}$$

be two absorber configurations with ${\gamma }_{j}^{\ell }$ defined in (2.10), and consider the following constants

$$\begin{equation}{c}_{0}=\underset{\underset{1{\leqslant}j{\leqslant}{m}_{2}}{1{\leqslant}i{\leqslant}{m}_{1}}}{\mathrm{min}}\left(\left\vert {\nu }_{i,1}^{0,a}\right\vert ,\left\vert {\nu }_{j,2}^{0,a}\right\vert \right)\end{equation} \tag{ 3.2 }$$

and

$$\begin{equation}{c}_{1}={\left(2\frac{{\left({m}_{1}+{m}_{2}-1\right)}^{2}}{{d}^{2}}+1\right)}^{\frac{1}{2}},\end{equation} \tag{ 3.3 }$$

where ${\nu }_{i,1}^{0,a}$ and ${\nu }_{j,2}^{0,a}$ are as defined in (2.12):

$$\begin{equation*}{\nu }_{i,1}^{0,a}={{\epsilon}}^{3}{\int }_{0}^{T}{\int }_{{B}_{i}^{1}}{\gamma }_{i}^{1}\left(\boldsymbol{\tau },t\right)\rho \left({z}_{i}^{1}+{\epsilon}{\tau }_{3},t\right)\mathrm{d}\boldsymbol{\tau }\mathrm{d}t\end{equation*}$$

and

$$\begin{equation*}{\nu }_{j,2}^{0,a}={{\epsilon}}^{3}{\int }_{0}^{T}{\int }_{{B}_{j}^{2}}{\gamma }_{j}^{2}\left(\boldsymbol{\tau },t\right)\rho \left({z}_{j}^{2}+{\epsilon}{\tau }_{3},t\right)\mathrm{d}\boldsymbol{\tau }\mathrm{d}t,\end{equation*}$$

and the function ρ is the solution of (2.6). Then take ϱ > 0 such that

$$\begin{equation}\underset{t\in \left[0,T\right]}{\mathrm{sup}}{\Vert}\rho \left(.,t\right){\parallel }_{{H}^{1}\left({\gamma }_{1},{\gamma }_{2}\right)}{\leqslant}\varrho ,\end{equation} \tag{ 3.4 }$$

where ${\Vert}\cdot {{\Vert}}_{{H}^{1}\left({\gamma }_{1},{\gamma }_{2}\right)}$ is the classical H¹ norm defined by

$$\begin{equation*}{\Vert}\rho \left(.,t\right){{\Vert}}_{{H}^{1}\left({\gamma }_{1},{\gamma }_{2}\right)}={\left(\vert \rho \left(.,t\right){\vert }_{{L}^{2}\left({\gamma }_{1},{\gamma }_{2}\right)}^{2}+\vert {\rho }_{z}\left(.,t\right){\vert }_{{L}^{2}\left({\gamma }_{1},{\gamma }_{2}\right)}^{2}\right)}^{\frac{1}{2}}.\end{equation*}$$

Now we are able to state our stability result for the absorbers' centers.

Theorem 3.1. Let p^k, for k = 1, 2, be the solutions of (1.3) corresponding to two different sets of absorption domains, ${D}_{j}^{k}={\mathbf{S}}_{j}^{k}+{\epsilon}{B}_{j}^{k},1{\leqslant}j{\leqslant}{m}_{k}$, and characterized by the configurations ${\left({\gamma }_{j}^{k},{\mathbf{S}}_{j}^{k}\right)}_{1{\leqslant}j{\leqslant}{m}_{k}}$ with ${\nu }_{j,k}^{0,a}\ne 0$. For k = 1, 2, we take ${f}^{k}{:=}{p}_{{\Gamma}{\times}\left(0,T\right)}^{k}$ to be the corresponding measurements on the boundary Σ = Γ × (0, T), and ${S}^{k}={\left\{{\mathbf{S}}_{j}^{k}\right\}}_{1{\leqslant}j{\leqslant}{m}_{k}}$. Assuming that ${\mathbf{S}}_{j}^{k}\in {{\Omega}}_{\alpha }$, then the following estimate holds,

$$\begin{equation*}{d}_{H}\left({S}^{1},{S}^{2}\right){\leqslant}2\underset{k=1,2}{\mathrm{max}}{\left[\frac{{c}_{1}\varrho \sqrt{T\vert {\Gamma}\vert }}{{c}_{0}{\delta }^{{m}_{3-k}-1}}{d}^{{m}_{1}+{m}_{2}-1}{\Vert}\;{f}^{2}-{f}^{1}{\parallel }_{{L}^{2}\left({\Sigma}\right)}+O\left({{\epsilon}}^{4}\right)\right]}^{\frac{1}{{m}_{k}}}.\end{equation*}$$

The proof of theorem 3.1 is based on the following lemma.

Lemma 3.2. Let ${P}^{r,k}={\left\{{P}_{j}^{r,k}\right\}}_{1{\leqslant}j{\leqslant}{m}_{k}}$, for r = a, b, be respectively the corresponding projections on the xy- and yz-planes of S^k. Under the assumptions of theorem 3.1, we have, for r = a, b and k = 1, 2,

$$\begin{equation}{d}_{H}\left({P}^{r,1},{P}^{r,2}\right){\leqslant}\underset{k=1,2}{\mathrm{max}}{\left[\frac{{c}_{1}\varrho \sqrt{T\vert {\Gamma}\vert }}{{c}_{0}{\delta }^{{m}_{3-k}-1}}{d}^{{m}_{1}+{m}_{2}-1}{\Vert}{f}^{2}-{f}^{1}{\parallel }_{{L}^{2}\left({\Sigma}\right)}+O\left({{\epsilon}}^{4}\right)\right]}^{\frac{1}{{m}_{k}}}.\end{equation} \tag{ 3.5 }$$

Proof. The proof is based on the idea of choosing specific wave functions. In fact, if we take for 1 ⩽ ℓ ⩽ m₂,

$$\begin{equation*}{{\Phi}}_{\ell }\left(x,y\right)=\prod _{m=1}^{{m}_{1}}\left(x+iy-{P}_{m}^{a,1}\right)\prod _{m\ne \ell }^{{m}_{2}}\left(x+iy-{P}_{m}^{a,2}\right)\end{equation*}$$

and

$$\begin{equation*}{{\Psi}}_{\ell }\left(\mathbf{\text{x}},t\right)=\rho \left(z,t\right){{\Phi}}_{\ell }\left(x,y\right),\end{equation*}$$

where ρ is defined in (2.6). Then, we observe that Ψ_ℓ is a solution of (2.1) satisfying (2.5). Thus,

$$\begin{equation}\mathcal{R}\left(f,{{\Psi}}_{\ell }\right)={\int }_{0}^{T}{\int }_{{\Gamma}}f\left(\xi ,t\right)\frac{\partial {{\Psi}}_{\ell }}{\partial \nu }\left(\xi ,t\right)\;\mathrm{d}\xi \mathrm{d}t,\end{equation} \tag{ 3.6 }$$

where $\mathcal{R}\left(f,{{\Psi}}_{\ell }\right)$ is defined in (2.3). Also, substituting Ψ_ℓ in (2.4) we get

$$\begin{align*}\hfill \mathcal{R}\left({f}^{2},{{\Psi}}_{\ell }\right)=& \sum _{j=1}^{{m}_{2}}{\int }_{0}^{T}{\int }_{{D}_{j}^{2}}{\sigma }_{j}^{2}\left(\mathbf{\text{x}}\right){g}^{2}\left(\mathbf{\text{x}},t\right)\rho \left(z,t\right)\prod _{m=1}^{{m}_{1}}\left(x+iy-{P}_{m}^{a,1}\right)\prod _{m\ne \ell }^{{m}_{2}}\left(x+iy-{P}_{m}^{a,2}\right)\mathrm{d}\mathbf{\text{x}}\mathrm{d}t\hfill \\ \hfill & +O\left({{\epsilon}}^{\kappa }\right).\hfill \end{align*}$$

Applying the change of variable ${D}_{j}^{2}={\mathbf{S}}_{j}^{2}+{\epsilon}{B}_{j}^{2}$, we obtain

$$\begin{align*}\hfill \mathcal{R}\left({f}^{2},{{\Psi}}_{\ell }\right)=& {{\epsilon}}^{3}\sum _{j=1}^{{m}_{2}}{\int }_{0}^{T}{\int }_{{B}_{j}^{2}}{\sigma }_{j}^{2}\left({\mathbf{S}}_{j}^{2}+{\epsilon}\boldsymbol{\tau }\right){g}^{2}\left({\mathbf{S}}_{j}^{2}+{\epsilon}\boldsymbol{\tau },t\right)\rho \left({z}_{j}^{2}+{\epsilon}{\tau }_{3},t\right)\hfill \\ \hfill & {\times}\prod _{m=1}^{{m}_{1}}\left({P}_{j}^{a,2}-{P}_{m}^{a,1}+{\epsilon}\left({\tau }_{1}+i{\tau }_{2}\right)\right)\prod _{m\ne \ell }^{{m}_{2}}\left({P}_{j}^{a,2}-{P}_{m}^{a,2}+{\epsilon}\left({\tau }_{1}+i{\tau }_{2}\right)\right)\mathrm{d}\boldsymbol{\tau }\mathrm{d}t+O\left({{\epsilon}}^{\kappa }\right).\hfill \end{align*}$$

Since κ > 4, we have

$$\begin{align*}\hfill \mathcal{R}\left({f}^{2},{{\Psi}}_{\ell }\right)& ={{\epsilon}}^{3}\sum _{j=1}^{{m}_{2}}{\int }_{0}^{T}{\int }_{{B}_{j}^{2}}{\sigma }_{j}^{2}\left({\mathbf{S}}_{j}^{2}+{\epsilon}\boldsymbol{\tau }\right){g}^{2}\left({\mathbf{S}}_{j}^{2}+{\epsilon}\boldsymbol{\tau },t\right)\rho \left({z}_{j}^{2}+{\epsilon}{\tau }_{3},t\right)\prod _{m=1}^{{m}_{1}}\left({P}_{j}^{a,2}-{P}_{m}^{a,1}\right)\hfill \\ \hfill & {\times}\prod _{m\ne \ell }^{{m}_{2}}\left({P}_{j}^{a,2}-{P}_{m}^{a,2}\right)\mathrm{d}\boldsymbol{\tau }\mathrm{d}t+O\left({{\epsilon}}^{4}\right)\hfill \\ \hfill & =\prod _{m=1}^{{m}_{1}}\left({P}_{\ell }^{a,2}-{P}_{m}^{a,1}\right)\prod _{m\ne \ell }^{{m}_{2}}\left({P}_{\ell }^{a,2}-{P}_{m}^{a,2}\right){\nu }_{\ell ,2}^{0,a}+O\left({{\epsilon}}^{4}\right).\hfill \end{align*}$$

Similarly, one has

$$\begin{align}\hfill \mathcal{R}\left({f}^{1},{{\Psi}}_{\ell }\right)& =\sum _{j=1}^{{m}_{1}}{{\epsilon}}^{3}{\int }_{0}^{T}{\int }_{{B}_{j}^{1}}{\sigma }_{j}^{1}\left({\mathbf{S}}_{j}^{1}+{\epsilon}\boldsymbol{\tau }\right){g}^{1}\left({\mathbf{S}}_{j}^{1}+{\epsilon}\boldsymbol{\tau },t\right)\rho \left({z}_{j}^{1}+{\epsilon}{\tau }_{3},t\right)\prod _{m=1}^{{m}_{1}}\left({P}_{j}^{a,1}-{P}_{m}^{a,1}\right)\hfill \\ \hfill & {\times}\prod _{m\ne \ell }^{{m}_{2}}\left({P}_{j}^{a,1}-{P}_{m}^{a,2}\right)\mathrm{d}\boldsymbol{\tau }\mathrm{d}t+O\left({{\epsilon}}^{4}\right)=O\left({{\epsilon}}^{4}\right).\hfill \end{align}$$

Subtracting the two previous sums, we get

$$\begin{equation*}\mathcal{R}\left({f}^{2}-{f}^{1},{{\Psi}}_{\ell }\right)=\prod _{m=1}^{{m}_{1}}\left({P}_{\ell }^{a,2}-{P}_{m}^{a,1}\right)\prod _{m\ne \ell }^{{m}_{2}}\left({P}_{\ell }^{a,2}-{P}_{m}^{a,2}\right){\nu }_{\ell ,2}^{0,a}+O\left({{\epsilon}}^{4}\right).\end{equation*}$$

Therefore, from (3.1), (3.2) and (3.6), we obtain

$$\begin{equation*}{c}_{0}\underset{1{\leqslant}m{\leqslant}{m}_{1}}{\mathrm{min}}{\Vert}{P}_{\ell }^{a,2}-{P}_{m}^{a,1}{{\Vert}}^{{m}_{1}}{\leqslant}\frac{1}{{\delta }^{{m}_{2}-1}}{\Vert}{f}^{2}-{f}^{1}{{\Vert}}_{{L}^{2}\left({\Sigma}\right)}{\Vert}\frac{\partial {{\Psi}}_{\ell }}{\partial \nu }{\Vert}{L}^{2}\left({\Sigma}\right)+O\left({{\epsilon}}^{4}\right).\end{equation*}$$

Moreover, from (3.3) and (3.4), one can easily check that

$$\begin{equation*}{\Vert}\frac{\partial {{\Psi}}_{\ell }}{\partial \nu }{\Vert}{L}^{2}\left({\Sigma}\right){\leqslant}{c}_{1}\varrho \sqrt{T\vert {\Gamma}\vert }{d}^{{m}_{1}+{m}_{2}-1}.\end{equation*}$$

Finally, we get

$$\begin{equation}\underset{1{\leqslant}\ell {\leqslant}{m}_{2};}{\mathrm{max}}\underset{1{\leqslant}m{\leqslant}{m}_{1}}{\mathrm{min}}{\Vert}{P}_{\ell }^{a,2}-{P}_{m}^{a,1}{\Vert}{\leqslant}{\left[\frac{{c}_{1}\varrho \sqrt{T\vert {\Gamma}\vert }}{{c}_{0}{\delta }^{{m}_{2}-1}}{d}^{{m}_{1}+{m}_{2}-1}{\Vert}{f}^{2}-{f}^{1}{{\Vert}}_{{L}^{2}\left({\Sigma}\right)}+O\left({{\epsilon}}^{4}\right)\right]}^{\frac{1}{{m}_{1}}}.\end{equation} \tag{ 3.7 }$$

Similarly, by considering, first,

$$\begin{equation*}{\tilde {{\Phi}}}_{\ell }\left(x,y\right)=\prod _{m\ne \ell }^{{m}_{1}}\left(x+iy-{P}_{m}^{a,1}\right)\prod _{m=1}^{{m}_{2}}\left(x+iy-{P}_{m}^{a,2}\right),\end{equation*}$$

and then,

$$\begin{equation*}{\tilde {{\Psi}}}_{\ell }\left(\mathbf{\text{x}},t\right)=\rho \left(z,t\right){\tilde {{\Phi}}}_{\ell }\left(x,y\right),\end{equation*}$$

we get

$$\begin{equation}\underset{1{\leqslant}\ell {\leqslant}{m}_{1};}{\mathrm{max}}\underset{1{\leqslant}m{\leqslant}{m}_{2}}{\mathrm{min}}{\Vert}{P}_{\ell }^{a,1}-{P}_{m}^{a,2}{\Vert}{\leqslant}{\left[\frac{{c}_{1}\varrho \sqrt{T\vert {\Gamma}\vert }}{{c}_{0}{\delta }^{{m}_{1}-1}}{d}^{{m}_{1}+{m}_{2}-1}{\Vert}{f}^{2}-{f}^{1}{{\Vert}}_{{L}^{2}\left({\Sigma}\right)}+O\left({{\epsilon}}^{4}\right)\right]}^{\frac{1}{{m}_{2}}}.\end{equation} \tag{ 3.8 }$$

Taking the maximum between (3.7) and (3.8), we get (3.5) for r = a. To prove the second inequality, we consider the following functions

$$\begin{equation*}{{\Psi}}_{\ell }\left(\mathbf{\text{x}},t\right)={\rho }_{1}\left(x,t\right){{\Phi}}_{\ell }\left(y,z\right),\quad {\tilde {{\Psi}}}_{\ell }\left(\mathbf{\text{x}},t\right)={\rho }_{1}\left(x,t\right){\tilde {{\Phi}}}_{\ell }\left(y,z\right),\end{equation*}$$

where ρ₁ is the solution of (2.6) with the space variable z replaced by x. □

Now we are able to give the proof of theorem 3.1.

Proof. Using the fact that

$$\begin{equation*}\underset{1{\leqslant}n{\leqslant}{m}_{1};}{\mathrm{max}}\underset{1{\leqslant}j{\leqslant}{m}_{2}}{\mathrm{min}}{\Vert}{\mathbf{S}}_{n}^{1}-{\mathbf{S}}_{j}^{2}{\Vert}{\leqslant}\underset{1{\leqslant}n{\leqslant}{m}_{1};}{\mathrm{max}}\underset{1{\leqslant}j{\leqslant}{m}_{2}}{\mathrm{min}}{\Vert}{P}_{n}^{a,1}-{P}_{j}^{a,2}{\Vert}+\underset{1{\leqslant}n{\leqslant}{m}_{1};}{\mathrm{max}}\underset{1{\leqslant}j{\leqslant}{m}_{2}}{\mathrm{min}}{\Vert}{P}_{n}^{b,1}-{P}_{j}^{b,2}{\Vert}\end{equation*}$$

and

$$\begin{equation*}\underset{1{\leqslant}n{\leqslant}{m}_{2};}{\mathrm{max}}\underset{1{\leqslant}j{\leqslant}{m}_{1}}{\mathrm{min}}{\Vert}{\mathbf{S}}_{n}^{2}-{\mathbf{S}}_{j}^{1}{\Vert}{\leqslant}\underset{1{\leqslant}n{\leqslant}{m}_{2};}{\mathrm{max}}\underset{1{\leqslant}j{\leqslant}{m}_{1}}{\mathrm{min}}{\Vert}{P}_{n}^{a,2}-{P}_{j}^{a,1}{\Vert}+\underset{1{\leqslant}n{\leqslant}{m}_{2};}{\mathrm{max}}\underset{1{\leqslant}j{\leqslant}{m}_{1}}{\mathrm{min}}{\Vert}{P}_{n}^{b,2}-{P}_{j}^{b,1}{\Vert},\end{equation*}$$

we have

$$\begin{equation*}{d}_{H}\left({S}^{1},{S}^{2}\right){\leqslant}{d}_{H}\left({P}_{1}^{a},{P}_{2}^{a}\right)+{d}_{H}\left({P}_{1}^{b},{P}_{2}^{b}\right).\end{equation*}$$

Using the previous lemma, we get the desired result. □

Remark 3.3 (Numerical simulations). The numerical implementation of our algorithm and some extensions of this work, in particular, the investigation of stability in case of variable sound speed, will be considered in some future work. For this reason, we leave the numerical analysis of this inverse problem for some future investigations. Moreover, we would like to compare our results to those obtained in [1] even if they were in the context of photoacoustics.