Identifying an acoustic source in a two-layered medium from multi-frequency phased or phaseless far-field patterns

The field of inverse scattering problems, which gains momentum from significant scientific and industrial applications, has witnessed remarkable growth in recent years [8, 3, 4]. Being an essential topic in the field of inverse problems, the research background for the layered medium is rooted in the observation that numerous natural and artificial materials possess a two-layered structure. For example, the skin has a dermis and an epidermis [2], and printed circuit boards have a copper layer and a dielectric layer [1]. In addition, many technical processes involve using two-layered media, such as high-speed printing [13], high-precision cutting, and thin-film solar cells. The exploration of two-layered media is crucial for comprehending their physical attributes and enhancing their performance across various applications. However, compared with single-layered media, the study of two-layered media is more challenging due to the intricate interactions between the two layers. Furthermore, reflection, refraction, and transmission of waves occur at these interfaces, which complicates the analysis and computation of the wave propagation in a layered medium.

The research background for the layered medium is extensive and intertwines multiple disciplines, including materials science, physics, chemistry, and engineering. Mathematically, the study of two-layered media has garnered enduring attention. The authors of [21, 22] have established increasing stability results for the inverse source problem (ISP) in the one-dimensional Helmholtz equation in a multi-layered medium. In [9], the authors have developed increasing stability estimates for the inverse source scattering problem in two-dimensional space and introduced a recursive Kaczmarz-Landweber iteration scheme with incomplete data. To identify point sources from multi-frequency sparse far-field patterns, a multi-step numerical scheme is proposed in [17]. For a deeper understanding of the direct and inverse scattering problems in a two-layered medium, we refer to [10, 14, 11, 16, 12] and the associated reference materials.

Although the ISP in two-layered media has been the focus of research interest thus far, most existing methods are either sampling-based (direct but qualitative) or iterative-type (quantitative but time-consuming). For the sake of rapid reconstruction with high precision, it is necessary to investigate novel imaging algorithms.

In this work, we aim to reconstruct the source function with compact support in a two-layered medium by utilizing the multi-frequency far-field data collected on the upper half-space. The problem under consideration is more sophisticated than that in an unbounded space. A reason accounting for this is that the medium consists of two distinct layers, each with its own specific properties, and the propagation of waves through this medium is influenced by both layers. Additionally, the availability of far-field models is limited to the upper half space, forming a major challenge due to the practical physical model. Moreover, another obstruction to solving the problem lies in the fact that the observation aperture is closely related to the refractive index (defined as the ratio of the speed above and below the interface). Thus, when the refractive index decreases, the aperture shrinks sharply. It is well known that the small observation aperture will significantly increase the ill-posedness. All the aforementioned difficulties constitute the challenges of the underlying problem.

Motivated by the Fourier method proposed for the ISP in the full space for the Helmholtz equation [7, 18] and the biharmonic equation [20], we shall investigate the feasibility of extending this method to the two-layer media case. A notable distinction to previous full-space problems is that the far-field measurements in our two-layered model are only accessible in the upper half-space. A computational advantage of the Fourier method is that we can directly determine the Fourier coefficients from the far-field data without prior knowledge of the source function. Therefore, it is crucial to establish a correspondence between the admissible wavenumbers and the Fourier coefficients. Nevertheless, we highlight that the distinct sound speeds above and below the interface result in different wavenumbers, which may render the existing Fourier method inapplicable. Hence, extending the full-space approach to its two-layer version is non-trivial, thus novel modifications deserve investigation.

In many practical applications, only the intensity information of the radiated field can be measured, whereas the phase information is difficult to access or completely unavailable. To address this issue, we introduce an easy-to-implement phase retrieval technique by incorporating auxiliary known reference source points into the inverse source system. This technique allows us to recover the phase information, thereby reformulating the phaseless problem into the standard ISP with phase information. Compared to the full-space case, a notable difference is that the configuration of reference source points plays an important role in the specific details of phase recovery. In a nutshell, although the deployment of the auxiliary point in the lower space offers computational convenience, it comes at the expense of practical engineering costs. In comparison, it is more convenient and realistic to place the auxiliary point sources in the upper half-space, but this naturally requires an immense endeavor to derive novel computational schemes.

Once the phase information is retrieved, the phaseless inverse source problem (PLISP) can be easily reduced to its phased version, which can be solved using the Fourier method introduced here. Instead of delving into mathematical intricacies such as the stability and solvability of the equation system in the phase-retrieval process, we shall focus on feasible implementations and demonstrate its effectiveness through extensive numerical results conducted in two and three dimensions.

The remaining part of the work is organized as follows: In the next section, we present the forward and inverse source model with a two-layer media. The Fourier method for solving the phased problem is given in Section 3. Then, the phase retrieval algorithm for the phaseless inverse source problem is discussed in Section 4. Next, several numerical examples are provided in Section 5 to illustrate the effectiveness of the proposed methods, and the conclusions follow in Section 6.

Let us start this section with a mathematical description of the model setup. As shown in Figure 1, the flat plane $\Gamma_{0}:=\{x\in\mathbb{R}^{n}:x_{n}=0\}$ delimits the whole space $\mathbb{R}^{n}\,(n=2,3)$ into two half-spaces: the upper half space $\mathbb{R}_{+}^{n}:=\{x=(x_{1},x_{2},\cdots,x_{n})\in\mathbb{R}^{n}:x_{n}>0\}$ and the lower half space $\mathbb{R}_{-}^{n}:=\{x\in\mathbb{R}^{n}:x_{n}<0\}.$ The unbounded domains $\mathbb{R}_{\pm}^{n}$ are filled with homogeneous and isotropic acoustic material with wave speed $c_{\pm}>0$, respectively.

Given source function $S$, the radiated field $u$ is governed by the Helmholtz equation

$$\Delta u+k^{2}(x)u=S,\quad x\in\mathbb{R}^{n},$$

(1)

with $k(x)=k_{\pm}=\omega/c_{\pm}$ for $x\in\mathbb{R}^{n}_{\pm}$ and $\omega>0$ is the wave frequency. Throughout this article, we assume that $S$ is independent of $k_{\pm}$ and located in $\mathbb{R}_{-}^{n},$ with $\text{supp}\,S\subset V_{0}.$ Here, $V_{0}\subset\mathbb{R}_{-}^{n},n=2,3$ is some square ($n=2$) or cuboidal ($n=3$) domain. Under the assumption that $\Gamma_{0}$ is flat, the radiated field $u$ is supposed to satisfy the following Sommerfeld radiation condition

$$\lim\limits_{\lvert x\rvert\to\infty}\lvert x\rvert^{\frac{n-1}{2}}\left(\frac{\partial u}{\partial\lvert x\rvert}-\mathrm{i}k_{\pm}u\right)=0.$$

(2)

We want to point out that, for a general geometry setting on $\Gamma_{0}$, the Sommerfeld radiation condition (2) may no longer be valid to describe the asymptotic behavior of the radiated wave $u$ as $\lvert x\rvert\to\infty$ (see [12]).

In Figure 1, the observation aperture $\Theta_{c}^{n}$ is defined by

$$\Theta_{c}^{n}:=\begin{cases}(\theta_{c},\pi-\theta_{c}),&n=2,\\ \left(\theta_{c},{\pi}/{2}\right],&n=3,\end{cases}$$

with the critical angle given by

$$\theta_{c}=\begin{cases}\arccos\left(\dfrac{c_{+}}{c_{-}}\right),&\text{if}\ c_{-}>c_{+},\\ 0,&\text{if}\ c_{-}\leq c_{+}.\end{cases}$$

Correspondingly, the observation directions of interest are defined through:

	$\displaystyle\mathbb{S}_{c}^{1}$	$\displaystyle:=\left\{\hat{x}\in\mathbb{S}^{1}:\hat{x}=(\cos\theta,\sin\theta),\ \theta\in\Theta_{c}^{2}\right\},$
	$\displaystyle\mathbb{S}_{c}^{2}$	$\displaystyle:=\left\{\hat{x}\in\mathbb{S}^{2}:\hat{x}=(\cos\phi\cos\theta,\sin\phi\cos\theta,\sin\theta),\ \theta\in\Theta_{c}^{3},\ \phi\in[0,2\pi)\right\},$

respectively. For later use, the transmitted direction $\hat{x}^{t}\in\mathbb{S}^{n-1}\cap\mathbb{R}_{+}^{n}$ for $\hat{x}\in\mathbb{S}_{c}^{n}$ is denoted as

$\displaystyle\hat{x}^{t}:=\left\{\begin{aligned} &\left(\frac{c_{-}}{c_{+}}\cos\theta,\sqrt{1-\frac{c_{-}^{2}}{c_{+}^{2}}\cos^{2}\theta}\right),&&n=2,\\ &\left(\frac{c_{-}}{c_{+}}\cos\theta\cos\phi,\frac{c_{-}}{c_{+}}\cos\theta\sin\phi,\sqrt{1-\frac{c_{-}^{2}}{c_{+}^{2}}\cos^{2}\theta}\right),&&n=3.\end{aligned}\right.$

For $x=(x_{1},x_{2},\cdots,x_{n})\in\mathbb{R}^{n}(n=2,3),$ we also define $\tilde{x}=(x_{1},x_{2},\cdots,x_{n-1})\in\mathbb{R}^{n-1}$ and $x^{s}:=(x_{1},x_{2},\cdots,-x_{n})\in\mathbb{R}^{n}$.

It is well known that the solution to (1)–(2) is given by

$\displaystyle u(x,\omega)=-\int_{V_{0}}\Phi_{\omega}(x,y)S(y)\mathrm{d}y,$

(3)

where $\Phi_{\omega}(x,y)$ is the outgoing Green’s function to the following transmission problem:

$\displaystyle\left\{\begin{aligned} &\Delta_{x}\Phi_{\omega}(x,y)+k(x)^{2}\Phi_{\omega}(x,y)=-\delta_{y}(x),\quad x\in\mathbb{R}_{+}^{n}\cup\mathbb{R}_{-}^{n},\\ &\left[\Phi_{\omega}(x,y)\right]=\left[\partial_{x_{n}}\Phi_{\omega}(x,y)\right]=0,\quad x\in\Gamma_{0}.\end{aligned}\right.$

(4)

with, in each domain, the Rellich-Sommerfeld radiation condition satisfied:

$\displaystyle\lim\limits_{r=\lvert x\rvert\to\infty}r^{\frac{n-1}{2}}\left(\partial_{r}{\Phi}_{\omega}-\mathrm{i}k_{\pm}\Phi_{\omega}\right)=0.$

In (4), $[\,\cdot\,]$ denotes the jump across the interface $\Gamma_{0},$ and $\delta_{y}(x):=\delta(x-y)$ is the Dirac function in $\mathbb{R}^{n}.$ As $\Gamma_{0}$ considered here is a flat plane, the Green function has the following explicit expressions [17, 15, 12]:

$\displaystyle\Phi_{\omega}(x,y)=\left\{\begin{aligned} &\frac{1}{2\pi}\int_{\mathbb{R}^{n-1}}\frac{\mathrm{e}^{-\gamma_{1}x_{n}+\gamma_{2}y_{n}}}{\gamma_{1}+\gamma_{2}}\mathrm{e}^{\mathrm{i}\xi\cdot(\tilde{x}-\tilde{y})}\mathrm{d}\xi,&&y\in\mathbb{R}_{-}^{n},\\ &\frac{1}{2\pi}\int_{\mathbb{R}^{n-1}}\frac{\mathrm{e}^{-\gamma_{1}\lvert x_{n}-y_{n}\rvert}+\gamma\mathrm{e}^{-\gamma_{1}(x_{n}+y_{n})}}{2\gamma_{1}}\mathrm{e}^{\mathrm{i}\xi\cdot(\tilde{x}-\tilde{y})}\mathrm{d}\xi,&&y\in\mathbb{R}_{-}^{n},\end{aligned}\right.$

(5)

where $\gamma_{i}:=\left(\lvert\xi\rvert^{2}-k_{i}^{2}\right)^{1/2}$ has non-negative real part and non-positive imaginary part, and $\gamma:=(\gamma_{1}-\gamma_{2})/(\gamma_{1}+\gamma_{2})$. Further, the layered Green’s function has the following asymptotic behavior [17, 15, 12]:

$$\Phi_{\omega}(x,y)=\beta_{n}\frac{\mathrm{e}^{\mathrm{i}k_{+}\lvert x\rvert}}{\lvert x\rvert^{\frac{n-1}{2}}}\left\{\Phi_{\omega}^{\infty}(\hat{x},y)+o(1)\right\},\quad\text{as}\quad|x\rvert\to\infty,$$

(6)

which holds uniformly for all directions $\hat{x}\in\mathbb{S}_{c}^{n-1},$ with

$$\beta_{n}=\frac{\mathrm{e}^{\mathrm{i}\pi/4}}{\sqrt{8k_{+}\pi}}\left(\mathrm{e}^{-\mathrm{i}\frac{\pi}{4}}\sqrt{\frac{k_{+}}{2\pi}}\right)^{n-2},$$

and

$\displaystyle\Phi_{\omega}^{\infty}(\hat{x},y)=\left\{\begin{aligned} &T(\theta)\mathrm{e}^{-\mathrm{i}k_{-}\hat{x}^{t}\cdot y},&&y\in\mathbb{R}_{-}^{n},\\ &H(\theta)\mathrm{e}^{-\mathrm{i}k_{+}\hat{x}\cdot y^{s}}+\mathrm{e}^{-\mathrm{i}k_{+}\hat{x}\cdot y},&&y\in\mathbb{R}_{+}^{n},\end{aligned}\right.\quad\hat{x}\in\mathbb{S}_{c}^{n-1}.$

(7)

Here,

	$\displaystyle T(\theta)$	$\displaystyle=\frac{2\sin\theta}{\sin\theta+\sqrt{c_{+}^{2}/c_{-}^{2}-\cos^{2}\theta}},\quad\theta\in\Theta_{c}^{n},$
	$\displaystyle H(\theta)$	$\displaystyle=\frac{\sin\theta-\sqrt{c_{+}^{2}/c_{-}^{2}-\cos^{2}\theta}}{\sin\theta+\sqrt{c_{+}^{2}/c_{-}^{2}-\cos^{2}\theta}},\quad\theta\in\Theta_{c}^{n}.$

Accordingly, for $x\in\mathbb{R}_{+}^{n},$ $u(x,\omega)$ admits the following asymptotic expansion:

$\displaystyle u(x,\omega)=\beta_{n}\frac{\mathrm{e}^{\mathrm{i}k_{+}\lvert x\rvert}}{\lvert x\rvert^{\frac{n-1}{2}}}\left\{u^{\infty}(\hat{x},\omega)+o(1)\right\},\quad\text{as}\ |x\rvert\to\infty,$

(8)

where $u^{\infty}(\hat{x},\omega)$ defined on the upper half unit sphere/circle $\mathbb{S}_{+}^{n-1}$ is known as the far-field pattern (scattering amplitude) with $\hat{x}$ being the observation direction:

$\displaystyle u^{\infty}(\hat{x},\omega)=T(\theta)\int_{V_{0}}\mathrm{e}^{-\mathrm{i}k_{-}\hat{x}^{t}\cdot y}S(y)\mathrm{d}y.$

(9)

Let $\Omega_{N}=\{\omega_{j}\}_{j=1}^{N},\,(N\in\mathbb{N})$ be an admissible set consisting of a finite number of frequencies, and we are now to propose the two inverse problems under consideration as follows:

To tackle 2, we adopt a two-stage strategy. Initially, we introduce several appropriate auxiliary source points into the inverse system to recover the phase information. Subsequently, we transform 2 into 1.

Before heading for the inverse problem, we introduce several notations and the relevant Sobolev spaces in the rest of this section. Without loss of generality, we take $a>0,\,L>0$ such that the cubic domain

$$V_{0}=\left(-\frac{a}{2},\frac{a}{2}\right)^{n-1}\times\left(-L,0\right),$$

contains the support of the source function $S,$ i.e., $\text{supp }S\subset V_{0}.$ For $\sigma>0,$ the periodic Sobolev space $H^{\sigma}(V_{0})$ is defined by

$$H^{\sigma}(V_{0})=\{v\in L^{2}(V_{0}):\|v\|_{\sigma}<\infty\},$$

equipped with the norm

$$\|v\|_{\sigma}=\left(\sum_{\vec{l}\in\mathbb{Z}^{n}}\left(1+{\lvert{\vec{l}}\rvert}^{2}\right)\lvert\hat{v}_{\vec{l}}\,\rvert^{2}\right)^{1/2},$$

where

$$\hat{v}_{\vec{l}}=\frac{1}{a^{n}}\int_{V_{0}}v(x)\overline{\phi_{\vec{l}}\,(x)}\mathrm{d}x$$

is the Fourier coefficient and the Fourier basis function $\phi_{\vec{l}}\,(x)$ is given by

$$\phi_{\vec{l}}\,(x)=\mathrm{e}^{\mathrm{i}\frac{2\pi}{a}\vec{l}\cdot x},\quad\vec{l}=(l_{1},l_{2},\cdots,l_{n})\in\mathbb{Z}^{n}.$$

Within this preparation, the main idea of this article is to approximate the source function $S$ by

$\displaystyle S_{N}(x)=\sum_{\lvert\vec{l}\rvert_{\infty}\leq N}\hat{s}_{\vec{l}}\,\phi_{\vec{l}}\,(x),$

(10)

where $\hat{s}_{\vec{l}}$ is the Fourier coefficient defined by

$\displaystyle\hat{s}_{\vec{l}}=\frac{1}{a^{n}}\int_{V_{0}}S(x)\overline{\phi_{\vec{l}}\,(x)}\mathrm{d}x.$

(11)

For $S_{N}$ defined by (10), we have the following approximation result ([5]):

In this section, we shall develop a Fourier method to deal with 1. To appropriately utilize the multi-frequencies data, we are supposed to introduce the admissible frequencies.

Moreover, the admissible frequencies set $\Omega_{N}$, the admissible wave number set $\mathbb{K}_{+,N}$ in the upper half-space $\mathbb{R}_{+}^{n}$, and the observation angles are uniquely determined, i.e.,

$$\Omega_{N}=c_{-}\mathbb{K}_{-,N},\quad\mathbb{K}_{+,N}=\Omega_{N}/c_{+},\quad\frac{c_{-}}{c_{+}}\cos\theta=\frac{l_{1}}{\lvert{l}\rvert}.$$

We now deliver the following uniqueness theorem:

Combining (10), (17), and (18) gives the truncated Fourier series $S_{N}$ of the form

$\displaystyle S_{N}(x)=\hat{s}_{\vec{0}}+\sum_{1\leq\lvert\vec{l}\rvert_{\infty}\leq N}\hat{s}_{\vec{l}}\,\phi_{\vec{l}}\,(x),$

(19)

which can be viewed as an approximation to the source function.

In the following, instead of investigating the mathematical justification with regard to the uniqueness and stability of the proposed scheme, we shall head for an intriguing numerical exploration, and further consider the case where the measured data is phaseless.

The strategy to reconstruct the source from the phaseless data (2) can be divided into two steps: In the first step, we retrieve the phase information by adding several artificial source points to the setup, which transform the PLISP into its phased problem; In the second step, the Fourier method developed in Section 3 can be utilized to produce the final source reconstruction. We focus our concentration on the first step in this section.

Motivated by the phase retrieval technique developed in [6, 7, 19], we shall propose a novel formula to retrieve the phase information by introducing appropriate auxiliary source points to the model setup. Noticing (7), the locations of the auxiliary artificial source points would have an important influence on this process.

We now present the method to recover the phase information. For each observation direction $\hat{x}\in\mathbb{S}^{n-1},\,n=2,3,$ we take two points $z_{j}:=\alpha_{j}\hat{x}$ with $j=1,2,\,\alpha_{j}\in\mathbb{R}.$ Then the far-field pattern of the fundamental solution at points $z_{j}$ is given by

$\displaystyle\Phi_{\omega}^{\infty}(\hat{x},z_{j})=\left\{\begin{aligned} &T(\theta)\mathrm{e}^{-\mathrm{i}k_{-}\hat{x}^{t}\cdot z_{j}},&&z_{j}\in\mathbb{R}_{-}^{n},\\ &H(\theta)\mathrm{e}^{-\mathrm{i}k_{+}\hat{x}\cdot z_{j}^{s}}+\mathrm{e}^{-\mathrm{i}k_{+}\hat{x}\cdot z_{j}},&&z_{j}\in\mathbb{R}_{+}^{n},\end{aligned}\right.\quad\hat{x}\in\mathbb{S}_{c}^{n-1}.$

(20)

By defining the scaling parameters

$\displaystyle c_{j}=\frac{\|u^{\infty}(\cdot,\omega)\|_{\infty}}{\|\Phi_{\omega}^{\infty}(\cdot,z_{j})\|_{\infty}},\quad j=1,2,$

(21)

with $\|\cdot\|_{\infty}=\|\cdot\|_{L^{\infty}(\mathbb{S}^{n-1})}$, the auxiliary function $\phi_{j}(x,\omega)=-c_{j}\Phi_{\omega}(x,z_{j})$ satisfies the following inhomogeneous Helmholtz equation

$$\Delta\phi_{j}(x;\omega)+k(x)^{2}\phi_{j}(x)=c_{j}\delta_{z_{j}}(x)\quad\text{in }\ \mathbb{R}^{n}_{+}\cup\mathbb{R}^{n}_{-}.$$

By denoting $v_{j}=u+\phi_{j}\,(j=1,2)$, we derive that $v_{j}$ is the unique solution to

$$\begin{cases}\left(\Delta+k(x)^{2}\right)v_{j}(x,\omega)=\left(S+c_{j}\right)(x),\ \ {\rm in}\ \mathbb{R}^{n},\\ \lim\limits_{r=\lvert x\rvert\to\infty}r^{\frac{n-1}{2}}\left(\partial_{r}{v_{j}}(x,\omega)-\mathrm{i}k_{\pm}v_{j}(x,\omega)\right)=0.\end{cases}$$

Furthermore, the far-field pattern of $v_{j}$ is given by

	$\displaystyle v_{j}^{\infty}$	$\displaystyle=u^{\infty}+\phi_{j}^{\infty}$
		$\displaystyle=u^{\infty}-c_{j}\Phi_{\omega}^{\infty}$
		$\displaystyle=u^{\infty}-c_{j}\begin{cases}T(\theta)\mathrm{e}^{-\mathrm{i}k_{-}\hat{x}^{t}\cdot z_{j}},&z_{j}\in\mathbb{R}_{-}^{n},\\ H(\theta)\mathrm{e}^{-\mathrm{i}k_{+}\hat{x}\cdot z_{j}^{s}}+\mathrm{e}^{-\mathrm{i}k_{+}\hat{x}\cdot z_{j}},&z_{j}\in\mathbb{R}_{+}^{n},\end{cases}\quad j=1,2.$		(22)

We now turn to the following phase retrieval problem:

For simplicity, we denote

$$u^{\infty}=u^{\infty}(\hat{x},\omega),\quad v_{j}^{\infty}=v_{j}^{\infty}(\hat{x},\omega).$$

Then we derive from (22) that

$\displaystyle v_{j}^{\infty}=\Re u^{\infty}-c_{j}\Re\Phi_{\omega,j}^{\infty}+\mathrm{i}\left(\Im u^{\infty}-c_{j}\Im\Phi_{\omega,j}^{\infty}\right),$

with

	$\displaystyle\Re\Phi_{\omega,j}^{\infty}=\begin{cases}T(\theta)\cos(k_{-}\hat{x}^{t}\cdot z_{j}),&z_{j}\in\mathbb{R}_{-}^{n},\\ H(\theta)\cos(k_{+}\hat{x}\cdot z_{j}^{s})+\cos(k_{+}\hat{x}\cdot z_{j}),&z_{j}\in\mathbb{R}_{+}^{n},\end{cases}$
	$\displaystyle\Im\Phi_{\omega,j}^{\infty}=\begin{cases}-T(\theta)\sin(k_{-}\hat{x}^{t}\cdot z_{j}),&z_{j}\in\mathbb{R}_{-}^{n},\\ -H(\theta)\sin(k_{+}\hat{x}\cdot z_{j}^{s})-\sin(k_{+}\hat{x}\cdot z_{j}),&z_{j}\in\mathbb{R}_{+}^{n},\end{cases}$

Furthermore, we derive that

		$\displaystyle\lvert u^{\infty}\rvert^{2}=\left(\Re u^{\infty}\right)^{2}+\left(\Im u^{\infty}\right)^{2},$		(23)
		$\displaystyle\lvert v_{j}^{\infty}\rvert^{2}=\left(\Re u^{\infty}-c_{j}\Re\Phi_{\omega,j}^{\infty}\right)^{2}+\left(\Im u^{\infty}-c_{j}\Im\Phi_{\omega,j}^{\infty}\right)^{2}.$		(24)

Subtracting (23) from (24) gives that

$\displaystyle\Re\Phi_{\omega,j}^{\infty}\Re u^{\infty}+\Im\Phi_{\omega,j}^{\infty}\Im u^{\infty}=f_{j},\quad j=1,2,$

(25)

with

$$f_{j}=-\frac{1}{2c_{j}}\left(\lvert v_{j}^{\infty}\rvert^{2}-\lvert u^{\infty}\rvert^{2}-c_{j}^{2}\lvert\Phi_{\omega,j}^{\infty}\rvert^{2}\right).$$

Once (25) is solvable, the phase information of $u^{\infty}$ can be retrieved correspondingly, which further leads to the phased far field:

$\displaystyle u^{\infty}=\frac{\det D^{R}}{\det D}+\mathrm{i}\frac{\det D^{I}}{\det D},$

(26)

with

$$D=\begin{bmatrix}\Re\Phi_{\omega,1}^{\infty}&\Im\Phi_{\omega,1}^{\infty}\\ \Re\Phi_{\omega,2}^{\infty}&\Im\Phi_{\omega,2}^{\infty}\end{bmatrix},\quad D^{R}=\begin{bmatrix}f_{1}&\Im\Phi_{\omega,1}^{\infty}\\ f_{2}&\Im\Phi_{\omega,2}^{\infty}\end{bmatrix},\quad D^{I}=\begin{bmatrix}\Re\Phi_{\omega,1}^{\infty}&f_{1}\\ \Re\Phi_{\omega,2}^{\infty}&f_{2}\end{bmatrix}.$$

To ensure the unique solvability of (25), the determination of $D$ should be non-zero, which can be fulfilled by choosing $\alpha_{j},\,j=1,2,$ properly. Instead of heading for the mathematical argument to clarify the stability of the phase retrieval technique, we shall illustrate the performance of the proposed through several numerical examples. Nevertheless, the mathematical analysis to prove the stability of phase retrieval can be obtained similarly to that in [7].

By solving (25), we can obtain the phase information of the radiated field. We can then employ the Fourier method developed in Section 3 to reconstruct the source function.

In this section, we shall conduct several numerical experiments to verify the algorithms proposed in this paper. Both two- and three-dimensional cases are considered. The forward problem is solved by direct integration. The admissible wavenumber set in the lower half-space $\mathbb{R}_{-}^{n}$ is given by

$$\mathbb{K}_{-,N}=\left\{k:k=2\pi\lvert\vec{l}\rvert,\ \vec{l}\in\mathbb{L}_{N}\right\}\cup\left\{2\pi\lambda\right\},\ \lambda=10^{-3},$$

with

$$\mathbb{L}_{N}=\left\{\vec{l}:\vec{l}=(l_{1},l_{2},\cdots,l_{n})\in\mathbb{Z}^{n},\ 1\leq\lvert\vec{l}\rvert_{\infty}\leq N,l_{n}>0,\ \theta_{\vec{l}}\in\Theta_{c}^{n}\right\},$$

and $N$ is chosen to be 50 for convenience. Once $\mathbb{K}_{-,N}$ is chosen properly, the admissible wave frequency $\Omega_{N}$ and the admissible wave number $\mathbb{K}_{+,N}$ in $\mathbb{R}_{+}^{n}$ are set to be

$\displaystyle\Omega_{N}=\left\{\omega:\omega=k_{-}c_{-},\ k_{-}\in\mathbb{K}_{-,N}\right\},\quad\mathbb{K}_{+,N}=\Omega_{N}/c_{+}.$

In the following examples, the wave speed $c_{\pm}$ is taken to be $c_{-}=2,$ $c_{+}=c_{-}-\pi/1000.$ With the admissible wave number set chosen above, we obtain the synthesis far-field data from (9) for $\text{supp}f\subset V_{0}$. We use the Gauss quadrature to calculate the volume integrals over the $100\times 100$ or $50^{3}$ Gauss-Legendre points for $V_{0}\subset\mathbb{R}^{n},\,n=2,3.$ Specifically, $V_{0}$ is chosen to be $[-0.5,0.5]^{n-1}\times[-0.5,0]$.

To better exhibit the reconstruction, we divide this section into two subsections. In Section 5.1, we show the performance of the phase retrieval technique, both in two- and three-dimensional space. In Section 5.2, we shall verify the effectiveness of the Fourier method. The exact source functions are chosen to be

	$\displaystyle S_{\rm 2D}(x)$	$\displaystyle=1.1\exp\left(-200((x_{1}-0.01)^{2}+(x_{2}+0.38)^{2})\right)$
		$\displaystyle\quad-100\left((x_{2}+0.5)^{2}-x_{1}^{2}\right)\exp\left(-90\left(x_{1}^{2}+(x_{2}+0.5)^{2}\right)\right)$

for $n=2,$ and

	$\displaystyle S_{\rm 3D}(x)$	$\displaystyle=1.1\exp\left(-200((x_{1}-0.01)^{2}+(x_{2}-0.12)^{2}+(x_{3}+0.5)^{2})\right)$
		$\displaystyle\quad-100\left(x_{2}^{2}-x_{1}^{2}\right)\exp\left(-90\left(x_{1}^{2}+x_{2}^{2}+(x_{3}+0.5)^{2}\right)\right)$

for $n=3.$