¹¹footnotetext: Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China. [email protected]

Single Mode Multi-frequency Factorization Method for the Inverse Source Problem in Acoustic Waveguides

Shixu Meng¹

Abstract.

This paper investigates the inverse source problem with a single propagating mode at multiple frequencies in an acoustic waveguide. The goal is to provide both theoretical justifications and efficient algorithms for imaging extended sources using the sampling methods. In contrast to the existing far/near field operator based on the integral over the space variable in the sampling methods, a multi-frequency far-field operator is introduced based on the integral over the frequency variable. This far-field operator is defined in a way to incorporate the possibly non-linear dispersion relation, a unique feature in waveguides. The factorization method is deployed to establish a rigorous characterization of the range support which is the support of source in the direction of wave propagation. A related factorization-based sampling method is also discussed. These sampling methods are shown to be capable of imaging the range support of the source. Numerical examples are provided to illustrate the performance of the sampling methods, including an example to image a complete sound-soft block.

Key Words. inverse source problem, waveguide, multi-frequency, sampling method, Helmholtz equation.

1. Introduction

Inverse scattering is of great importance in non-destructive testing, medical imaging, geophysical exploration and numerous problems associated with target identification. In the last thirty years, sampling methods such as the linear sampling method [18], the factorization method [23] and their extensions have attracted a lot of interests. These sampling methods do not require a priori information about the scattering objects, and provide both theoretical justifications and robust numerical algorithms. There have been recent interests in inverse scattering for waveguides, mainly motivated by their applications in ocean acoustics, non-destructive testing of slender structures, imaging in and of tunnels [3, 22, 31]. One of the early works is the generalized dual space indicator method for underwater imaging [35]. The linear sampling method and factorization method were studied in acoustic waveguides [2, 10, 11, 27] and in elastic waveguides [9, 12]. Theories and applications of sampling methods have been extended to periodic waveguides [8, 33], a time domain linear sampling method [28], sampling methods in electromagnetic waveguide [26, 29]. It is also worth mentioning related imaging methods such as the time migration imaging method in an acoustic terminating waveguide [34], in an electromagnetic waveguide [16] and the time reversal imaging in an electromagnetic terminating waveguide [7]. Relations between the time migration imaging and the factorization method as well as a factorization-based imaging method is discussed in [6]. We also refer to [4] for related discussions on the connections among different imaging methods in the homogeneous space case. It is worth mentioning that a scatterer may be invisible in a waveguide that supports one propagating mode at a single frequency [20].

One difference between the waveguide and the homogeneous whole space is that the waveguide supports finitely many propagating modes and infinitely many evanescent modes at a fixed frequency. The theories of the linear sampling method and factorization method in waveguides [2, 5, 6, 10, 11, 27] have been established using both propagating modes and evanescent modes. Particularly in the far-field, the measurements stem from the propagating modes so that the linear sampling method and factorization method are usually implemented using only the propagating modes. The imaging results are good as far as there are sufficiently many propagating modes. However, an interesting observation made in [27] is that the linear sampling method may not be capable of imaging a complete block ( where a complete block could be an obstacle spanning the entire cross-section) using far-field measurements at a single frequency. This seems not a problem due to the sufficiency or deficiency of propagating modes. This motivates us to investigate further the role of propagating modes. One fact about the propagating modes is that each propagating mode propagates with different group velocity, indicated by the dispersion relation which is a unique feature in waveguide. Intuitively, a single propagating mode at multi-frequencies may help understand further the problem [27] of imaging a complete block with far-field measurements, since the phase (travel time) provides information on the bulk location. The efficiency of using a single propagating mode may be of practical interests in long-distance communication optical devices or tunnels. The above discussions lead to the following question in waveguide imaging: how to make use of a single propagating mode at multiple frequencies in a mathematically rigorous way?

While most of the existing works focus on the sampling methods at a single frequency with full-aperture measurements, the time migration imaging method [34] and the linear sampling method [5] investigated the limited-aperture problem with multiple frequencies. The time domain linear sampling method [28] amounts to use multiple frequency measurements. However, there seems no work to treat the problem with less measurements (such as backscattering measurements obtained from a single propagating mode) at multiple frequencies to lay the theoretical foundation for extended scatterers. This is one of the motivations of this paper: to provide theoretical foundation for extended sources using the sampling methods with a single propagating mode at multiple frequencies. This necessitates the study of the factorization method which provides theoretical justifications for extended sources. Another motivation is to help understand further the problem of imaging a complete block of the waveguide as pointed out in [27]. Though the single frequency linear sampling method may not be capable of imaging a complete block using far-field measurements, measurements with a single mode at multiple frequencies may help locate such a block efficiently. In a certain frequency range, the waveguide supports a single propagating mode so that the far-field measurements stem from such single propagating mode; such a singe mode waveguide was also studied [17] for applications such as invisibility and perfect reflectivity. Based on these two motivations, this paper considers the factorization method for the inverse source problem with single mode multi-frequency measurements for a waveguide that supports a single propagating mode; the measurements are made at one opening of the waveguide, and we call these “backscattering” measurements (where “” is to distinguish the inverse source problem from the inverse scattering problem). A related factorization-based sampling method would also be discussed. These sampling methods are shown to be capable of imaging the range support of the source, but not the support in the cross-section direction. With a single propagating mode, it seems to us that only the support in one direction could be imaged. It was similarly reported in [21] that the support in one direction can be imaged with far-field measurements at a single direction (and its opposite direction).

The methods in this paper inherit the tradition of the classical sampling methods to define a data operator in the far-field. However, in contrast to the existing far/near field operator based on the integral over the space variable, a multi-frequency far-field operator is introduced based on the integral over the frequency variable. This operator has a symmetric factorization which allows us to investigate the factorization method and the factorization-based sampling method. Sampling methods to treat the inverse source problems in the whole space were studied in [1, 21, 25], including the factorization method and direct/orthogonal sampling method. One difference between the waveguide imaging in this paper and the imaging in the whole space arises from the dispersion relation, a unique feature in waveguide. This brings an additional difficulty in defining or factorizing the multi-frequency far-field operator. This paper overcomes this difficulty by incorporating the dispersion relation into the definition of the far-field operator. Another difference is that the measurements under consideration are “backscattering” measurements (in contrast to measurements obtained in two opposite directions in the homogeneous space case [21]); we shall show how to use these less measurements to design appropriate multi-frequency operators to deploy the factorization method. As demonstrated in [1, 21, 25] and the references therein, one may hope to image an approximate support of the source with only one measurement point/direction for imaging in the whole space. This paper reveals a similar result in the waveguide imaging problem: it turns out that the range support of the source can be reconstructed with a single propagating mode at multiple frequencies.

This paper is further organized as follows. Section 2 provides the mathematical model for the inverse source problem in a waveguide. We also briefly discuss the concept of propagating/evanescent modes and dispersion relation, as well as the Green function and the forward problem. The multi-frequency far-field operator is proposed in Section 3, and it is shown that it has a symmetric factorization where the factorized middle operator is coercive. Section 4.1 is devoted to the factorization method which provides a mathematically rigorous characterization of the source support. Motivated by the factorization of the multi-frequency far-field operator, a factorization-based sampling method is then discussed in Section 4.2. Finally numerical examples are provided in Section 5 to illustrate the performance of the two sampling methods to image extended sources. As an extension and application, a numerical example is also provided to illustrate the potential of the multi-frequency sampling methods to image a complete sound-soft block. Section 6 discusses another two multi-frequency operators for sampling methods which impose less restrictive condition on the source. Finally we conclude this paper with a conclusion in Section 7.

2. Mathematical model and problem setup

Let us now introduce the mathematical model of the inverse source problem. The waveguide in $\mathbb{R}^{d}$ is given by $W=(-\infty,\infty)\times\Sigma$ , where $\Sigma$ is the cross-section of the waveguide. In this paper we restrict ourself to the two dimensional case $d=2$ where $\Sigma=(0,|\Sigma|)$ with $|\Sigma|$ denoting the length of $\Sigma$ , and remark that the extension to the three dimensional case can follow similarly. We consider the model that the boundary of the waveguide is either sound-hard, sound-soft, or mixed (Dirichlet on part of the boundary and Neumann on the remaining boundary). Let us denote the support of the source $f$ by $D$ with the assumption that $f\in L^{\infty}(D)$ and $D$ is a domain with Lipschitz boundary. See Figure 1 for an illustration of the problem. For any $x\in\mathbb{R}^{d}$ , let $x=(x_{1},x_{\perp})$ be the coordinate of $x$ where $x_{1}$ is the variable in the range and $x_{\perp}$ is the variable in the cross-section.

Figure 1. Source

D

in the acoustic waveguide

W

Let the interrogating wave number be a positive number $k$ . Let $u^{s}(\cdot;k)$ be the wave field due to the source $f$ . The wave field $u^{s}(\cdot;k)$ belongs to $\widetilde{H}^{1}_{loc}(W)$ and satisfies

$\displaystyle\Delta u^{s}(\cdot;k)+k^{2}u^{s}(\cdot;k)={-}f\quad$	$\displaystyle\mbox{in}\quad W,$	(1)
$\displaystyle\mathcal{B}u^{s}(\cdot;k)=0\quad$	$\displaystyle\mbox{on}\quad\partial W,$	(2)
$\displaystyle u^{s}(\cdot;k)\quad$	$\displaystyle\mbox{satisfies a radiation condition},$	(3)

where the radiation condition is specified in Section 2.1, and $\mathcal{B}$ denotes either the Dirichlet, Neumann, or a mixed boundary condition, i.e.

\mathcal{B}u^{s}:=\Bigg{\{}\begin{array}[]{ccc}u^{s}=0&\mbox{Dirichlet}\\ \frac{\partial u^{s}}{\partial\nu}=0&\mbox{Neumann}\\ \mbox{Dirichlet on }\Gamma_{1},\mbox{ Neumann on }\Gamma_{2}&\mbox{mixed}\end{array},

(4)

with $\Gamma_{1}$ (resp. $\Gamma_{2}$ ) denoting either the top (resp. bottom) or the bottom (resp. top) boundary, $\nu$ is the unit outward normal, and $\widetilde{H}^{1}_{loc}(W)$ is the set

\{u:u|_{\Omega}\in H^{1}(\Omega),\forall\,{\Omega=(-b,b)\times\Sigma}\mbox{ s.t. }D\subset\Omega\}

with $H^{1}$ denoting the standard Sobolev space consisting of bounded $L^{2}$ functions whose derivatives belong to $L^{2}$ .

The goal is to reconstruct the support $D$ of the source with a single propagating mode at multiple frequencies. In the case that the waveguide supports one propagating mode, the use of one propagating mode at multiple frequencies amounts to the use of multi-frequency measurements at one measurement point. In particular, the inverse problem is to characterize the support $D$ of the source from the following non-vanishing measurements

•

“Backscattering” measurements:

$\{u^{s}_{p}({x^{*}};k):\quad{x^{*}}\in\{-a\}\times\Sigma,\,k\in K\}$ (5)

where $u^{s}_{p}$ represents the propagating part (far-field) of $u^{s}$ , ${x^{*}}$ is a measurement point at a measurement cross section $\{-a\}\times\Sigma$ , and $K$ is the range of wavenumber (i.e. frequency).

Here $K$ is the range of wavenumber such that the waveguide supports only one propagating mode (which will be specified later in details in Section 2.1). In practice, $a>0$ is such that the measurements are at the far-field of the waveguide (where the contribution from the evanescent part is negligible), this is feasible as the frequencies are always discrete numbers and the evanescent part will be negligible for some $a$ . In this paper we fix the notation that waves measured at $\{-a\}\times\Sigma$ are called “backscattering” measurements (to distinguish the inverse source problem from the inverse scattering problem) and waves measured at $\{a\}\times\Sigma$ are called “forward scattering” measurements. It is also possible to consider both the “backward and forward scattering” measurements, i.e. measurements at $\{\pm a\}\times\Sigma$ , however this type of measurements may not be applicable in applications where measurements can only be made at one opening of the waveguide. Nevertheless, a detailed remark is to be made in Section 6 concerning this type of measurements.

2.1. Propagating modes

Let us introduce the eigensystem $\{\psi_{n},\lambda_{n}^{2}\}_{n=1}^{\infty}$ (where we sort $0\leq\lambda_{n}\leq\lambda_{n+1}$ ) by

	$\displaystyle\Delta_{\perp}\psi_{n}+\lambda_{n}^{2}\psi_{n}=0\quad$		$\displaystyle\mbox{in}\quad\Sigma,$		(6)
	$\displaystyle\mathcal{B}_{\perp}\psi_{n}=0\quad$		$\displaystyle\mbox{on}\quad\partial\Sigma,$		(7)

where $\mathcal{B}_{\perp}$ (defined according to $\mathcal{B}$ in (4)) denotes either the Dirichlet, Neumann, or mixed boundary condition, i.e.

\mathcal{B}_{\perp}\psi_{n}=\Bigg{\{}\begin{array}[]{ccc}\psi_{n}=0&\mbox{Dirichlet}\\ \frac{\partial\psi_{n}}{\partial\nu_{\perp}}=0&\mbox{Neumann}\\ \mbox{Dirichlet on }\Gamma_{1\Sigma},\mbox{ Neumann on }\Gamma_{2\Sigma}&\mbox{mixed}\end{array},

where $\nu_{\perp}$ denotes the unit outward normal to $\partial\Sigma$ . Here $\Gamma_{1\Sigma}$ is either the top $\{x_{\perp}:x_{\perp}=|\Sigma|\}$ or the bottom $\{x_{\perp}:x_{\perp}=0\}$ part of $\partial\Sigma$ according to (4), and $\Gamma_{2\Sigma}$ is the remaining part of $\partial\Sigma$ .

For every fixed $k$ with $k^{2}\not\in\{\lambda_{n}^{2}\}_{n=1}^{\infty}$ , the following physical wave functions (can be found via separation of variables) are solutions to the Helmholtz equation (1) away from the support of the source

\Bigg{\{}\begin{array}[]{ccc}e^{\pm i\sqrt{k^{2}-\lambda_{n}^{2}}x_{1}}\psi_{n}(x_{\perp}),&k>\lambda_{n}\\ e^{-\sqrt{\lambda_{n}^{2}-k^{2}}|x_{1}|}\psi_{n}(x_{\perp}),&k<\lambda_{n}\end{array},\quad n=1,2,\cdots,\quad{|x_{1}|\gg 1},

(8)

here we have chosen the branch cut of the square root with positive imaginary part, and we call it a propagating mode if $k>\lambda_{n}$ and an evanescent mode if $k<\lambda_{n}$ .

For our problem, we consider the frequency range $(\lambda_{1},\lambda_{2})$ such that the waveguide supports only one propagating mode and we set $K:=(\lambda_{1},\lambda_{2})$ . For later purposes we note that the eigenfunction $\psi_{1}$ can be chosen as real-valued with the following explicit expression

\psi_{1}(x_{\perp})=\left\{\begin{array}[]{cccc}\sqrt{2/|\Sigma|}\sin(\pi x_{\perp}/|\Sigma|)&\mbox{Dirichlet}\\ \sqrt{1/|\Sigma|}&\mbox{Neumann}\\ \sqrt{2/|\Sigma|}\cos(\pi x_{\perp}/(2|\Sigma|))&\qquad\Gamma_{1\Sigma}=\{x_{\perp}:x_{\perp}=|\Sigma|\}\\ \sqrt{2/|\Sigma|}\sin(\pi x_{\perp}/(2|\Sigma|))&\qquad\Gamma_{1\Sigma}=\{x_{\perp}:x_{\perp}=0\}\end{array}\right..

We now give the definition of the radiation condition.

Definition 1.

We say a solution to the Helmholtz equation (1) satisfies the radiation condition if

u^{s}(x;k)=\sum_{n=1}^{N(k)}a_{n}e^{i\sqrt{k^{2}-\lambda_{n}^{2}}|x_{1}|}\psi_{n}(x_{\perp})+\sum_{n=N(k)+1}^{\infty}a_{n}e^{-\sqrt{\lambda_{n}^{2}-k^{2}}|x_{1}|}\psi_{n}(x_{\perp}),\quad|x_{1}|\gg 1,

(9)

where $N(k)$ is the number of the propagating modes, and $a_{n}$ are coefficients determined by $u^{s}(x;k)$ .

This definition would ensure that the forward problem (1)–(3) has at most one solution, see for instance [5, 10].

2.2. Dispersion relation

For each propagating mode $e^{\pm i\sqrt{k^{2}-\lambda_{n}^{2}}x_{1}}\psi_{n}(x_{\perp})$ , we identify the concept of group wave number by

\mu_{n}(k):=\sqrt{k^{2}-\lambda_{n}^{2}},

and we call the set $\{(k,\mu_{n}(k)):k\in(\lambda_{n},+\infty)\}$ the $n$ -th “dispersion relation”, which describes how this propagating mode propagates with respect to a given wave number $k$ . To illustrate the dispersion relation, in Figure 2 we plot $\{(k,\mu_{1}(k)):k\in(1,2)\}$ for a two-dimensional waveguide $(-\infty,\infty)\times(0,\pi)$ with Dirichlet boundary condition. For sake of completeness, we also define $\mu_{n}(k):=\sqrt{k^{2}-\lambda_{n}^{2}}$ (where the chosen branch cut of the square root is with positive imaginary part) for each evanescent modes when $k<\lambda_{n}$ .

This dispersion relation is a unique feature of waveguide (which is in contrast to the linear dispersion appeared in the far-field patterns in an homogeneous background) and poses difficulty in the single mode multi-frequency sampling method.

Refer to caption — Figure 2. The plot of $\{(k,\mu_{1}(k)):k\in(1,2)\}$ for a two-dimensional waveguide $(-\infty,\infty)\times(0,\pi)$ with Dirichlet boundary condition.

2.3. Green function and forward problem

Now we introduce the Green function $G(x,y;k)$ for the waveguide $W$ that solves

$\displaystyle\Delta_{x}G(\cdot,y;k)+k^{2}G(\cdot,y;k)=-\delta_{y}\quad$	$\displaystyle\mbox{in}\quad W,$	(10)
$\displaystyle\mathcal{B}G(\cdot,y;k)=0\quad$	$\displaystyle\mbox{on}\quad\partial W,$	(11)
$\displaystyle G(\cdot,y;k)\quad$	$\displaystyle\mbox{satisfies the radiation condition},$	(12)

here the derivatives are with respect to variable $x$ . The Green function has the following series expansion [5, 10]

G(x,y;k)=\sum_{n=1}^{\infty}\frac{i}{2\mu_{n}(k)}\psi_{n}(x_{\perp})\psi_{n}(y_{\perp})e^{i\mu_{n}(k)|x_{1}-y_{1}|}.

(13)

Since the forward problem has at most one unique solution (see for instance [5, 10]), then the following lemma follows directly by verifying $u^{s}$ (14) is a solution to the forward problem via the volume integral, see for instance [19, Theorem 2.1].

Lemma 1.

There exists a unique solution $u^{s}\in\widetilde{H}^{1}_{loc}(W)$ to the forward problem (1)–(3). Furthermore it holds that

u^{s}(\cdot;k)=\int_{D}G(\cdot,y;k)f(y)\,\mbox{d}y.

(14)

3. The multi-frequency far-field operator and its factorization

The imaging methods under consideration are the qualitative methods or the sampling methods [13, 14, 15, 19], and the goal is to lay the theoretical foundation for extended sources and design efficient imaging algorithms with these measurements. To achieve this, we introduce the following multi-frequency far-field operator and study its properties in this section.

To begin with, set $\omega_{\sigma\gamma}:=\sqrt{\lambda_{1}^{2}+(\sigma-\gamma)^{2}}$ ( and details will be given in later context to motivate its definition). We introduce the single mode multi-frequency far-field operator $F:L^{2}(k_{-},k_{+})\to L^{2}(k_{-},k_{+})$ by

(Fg)(\sigma;x^{*}):=\int_{k_{-}}^{k_{+}}-i\mu_{1}(\omega_{\sigma\gamma})\Big{[}H(\sigma-\gamma)e^{i\theta}u_{p}^{s}\left(x^{*};\omega_{\sigma\gamma}\right)-H(\gamma-\sigma)e^{-i\theta}\overline{u_{p}^{s}\left(x^{*};\omega_{\sigma\gamma}\right)}\Big{]}g(\gamma)\,\mbox{d}\gamma,

(15)

$\sigma\in(k_{-},k_{+})$ , where $u_{p}^{s}$ represents the propagating part of the wave field $u^{s}$ , $H$ is the Heaviside function

H(s)=\Big{\{}\begin{array}[]{ccc}1,&s\geq 0\\ 0,&s<0\end{array},

$k_{-}=0$ and $k_{+}=\sqrt{\lambda_{2}^{2}-\lambda_{1}^{2}}$ is such that $\omega_{\sigma\gamma}\in(\lambda_{1},\lambda_{2}),\forall\sigma\not=\gamma$ . The assumption made on $f$ is that $e^{i\theta}f>c_{1}>0$ or $e^{i\theta}f<c_{2}<0$ with some $\theta\in[0,2\pi)$ . We remark that this assumption in the “backscattering” case can be relaxed if both “backward and forward scattering” measurements are available (see Section 6.1), or we can design another less obvious multi-frequency operator (see Section 6.2). We consider this case to fully illustrate the multi-frequency sampling methods in waveguide. Another note is that $\psi_{1}({{x}^{*}_{\perp}})\not=0$ since otherwise the measurement set (5) would be a vanishing one.

As can be seen, the possibly non-linear dispersion relation is incorporated into the definition of the far-field operator via $\omega_{\sigma\gamma}$ . In particular, one finds via a direct calculation that $\mu_{1}(\omega_{\sigma\gamma})={|\sigma-\gamma|}$ which would allow us to investigate the factorization method further. This is in contrast to the inverse source problem in the homogeneous whole space where the dispersion is already linear.

To deploy the factorization method based on the multi-frequency operator (15), we first introduce the operator $S:L^{2}(k_{-},k_{+})\to L^{2}(D)$ by

(S\varphi)(y):=\int_{k_{-}}^{k_{+}}\sqrt{\psi_{1}(y_{\perp})}e^{-i\sigma(y_{1}-x^{*}_{1})}\varphi(\sigma)\,\mbox{d}\sigma,\quad y\in D,

(16)

and the operator $T:L^{2}(D)\to L^{2}(D)$ by

(Th)(y):=\frac{e^{i\theta}\psi_{1}(x^{*}_{\perp})f(y)h(y)}{2},\quad y\in D.

(17)

We can directly derive the adjoint operator of $S$ , namely $S^{*}:L^{2}(D)\to L^{2}(k_{-},k_{+})$ by

(S^{*}h)(\sigma):=\int_{D}\overline{\sqrt{\psi_{1}(y_{\perp})}}e^{i\sigma(y_{1}-x^{*}_{1})}h(y)\,\mbox{d}y,\quad\sigma\in(k_{-},k_{+}).

(18)

We have chosen to keep the conjugate in (18) to formally retain the structure of the adjoint. Note that $\psi_{1}>0$ a.e. on the cross section, there is no confusion in the use of square root in (16) and (18). In the following, we give a factorization of the far-field operator. Note that $\psi_{1}>0$ is used in the proof of the following theorem as well.

Theorem 1.

The far-field operator $F:L^{2}(k_{-},k_{+})\to L^{2}(k_{-},k_{+})$ given by (15) can be factorized by

F=S^{*}TS,

(19)

where $S:L^{2}(k_{-},k_{+})\to L^{2}(D)$ , $T:L^{2}(D)\to L^{2}(D)$ , and $S^{*}:L^{2}(D)\to L^{2}(k_{-},k_{+})$ are given by (16), (17), and (18) respectively.

Proof.

Set $G_{p}$ the propagating part of the Green function $G$ . From the expression (14) of $u^{s}$ and (13), it follows that, for $\sigma\not=\gamma$ , the propagating part $u_{p}^{s}$ reads

	$\displaystyle u^{s}_{p}({x^{*}};\omega_{\sigma\gamma})$	$\displaystyle=$	$\displaystyle\int_{D}G_{p}\left({x^{*}},y;\omega_{\sigma\gamma}\right)f(y)\,\mbox{d}y$
		$\displaystyle=$	$\displaystyle\frac{i}{2\mu_{1}(\omega_{\sigma\gamma})}\int_{D}\psi_{1}({x^{}_{\perp}})\psi_{1}(y_{\perp})e^{i\mu_{1}(\omega_{\sigma\gamma})\|x^{}_{1}-y_{1}\|}f(y)\,\mbox{d}y,\quad\sigma\not=\gamma,$

by noting that $\mu_{1}(\omega_{\sigma\gamma})=|\sigma-\gamma|$ , we then continue to derive that

\displaystyle u^{s}_{p}({x^{*}};\omega_{\sigma\gamma})

\displaystyle=

\displaystyle\frac{i}{2|\sigma-\gamma|}\int_{D}\psi_{1}(x^{*}_{\perp})\psi_{1}(y_{\perp})e^{i|\sigma-\gamma||x^{*}_{1}-y_{1}|}f(y)\,\mbox{d}y,\quad\sigma\not=\gamma,

and consequently (note that $\psi_{1}$ and $e^{i\theta}f$ are real-valued)

			$\displaystyle\Big{[}H(\sigma-\gamma)e^{i\theta}u_{p}^{s}\left(x^{};\omega_{\sigma\gamma}\right)-H(\gamma-\sigma)e^{-i\theta}\overline{u_{p}^{s}\left(x^{};\omega_{\sigma\gamma}\right)}\Big{]}$
	$\displaystyle=$		$\displaystyle\frac{ie^{i\theta}}{2\|\sigma-\gamma\|}\int_{D}\psi_{1}(x^{}_{\perp})\psi_{1}(y_{\perp})e^{i(\sigma-\gamma)\|x^{}_{1}-y_{1}\|}f(y)\,\mbox{d}y,\quad\sigma\not=\gamma.$

Now from the definition of $F$ , it follows that

			$\displaystyle(Fg)(\sigma;x^{})=\int_{k_{-}}^{k_{+}}{-i}\mu_{1}(\omega_{\sigma\gamma})\Big{[}H(\sigma-\gamma)e^{i\theta}u_{p}^{s}\left(x^{};\omega_{\sigma\gamma}\right)-H(\gamma-\sigma)e^{-i\theta}\overline{u_{p}^{s}\left(x^{*};\omega_{\sigma\gamma}\right)}\Big{]}g(\gamma)\,\mbox{d}\gamma$
		$\displaystyle=$	$\displaystyle 1/2\int_{k_{-}}^{k_{+}}\int_{D}\psi_{1}(y_{\perp})e^{i(\sigma-\gamma)\|x^{}_{1}-y_{1}\|}e^{i\theta}{\psi_{1}(x^{}_{\perp})}f(y)g(\gamma)\,\mbox{d}y\,\mbox{d}\gamma$
		$\displaystyle=$	$\displaystyle\int_{D}e^{i\sigma\|x^{}_{1}-y_{1}\|}\overline{\sqrt{\psi_{1}(y_{\perp})}}\left[\left(\int_{k_{-}}^{k_{+}}\sqrt{\psi_{1}(y_{\perp})}e^{-i\gamma\|x^{}_{1}-y_{1}\|}g(\gamma)\,\mbox{d}\gamma\right)\frac{e^{i\theta}{\psi_{1}(x^{*}_{\perp})}f(y)}{2}\right]\,\mbox{d}y,$

where in the last step, we used the property that $\psi_{1}>0$ a.e. on the cross section so that $\sqrt{\psi_{1}(y_{\perp})}=\overline{\sqrt{\psi_{1}(y_{\perp})}}$ . The proof is then completed by recalling the definition of $S$ , $T$ , and $S^{*}$ in (16), (17), and (18) respectively. ∎

In the following, we study the properties of the factorized operators which are needed to investigate the factorization method in the next section.

Theorem 2.

The operator $T:L^{2}(D)\to L^{2}(D)$ given by (17) is self adjoint and coercive, and

|\langle Th,h\rangle_{L^{2}(D)}|\geq c\|h\|^{2}_{L^{2}(D)},

for some positive constant $c$ .

Proof.

Recall that $f\in L^{\infty}(D)$ with $e^{i\theta}f>c_{1}>0$ or $e^{i\theta}f<c_{2}<0$ with some $\theta\in[0,2\pi)$ , then $T$ is self adjoint; moreover from the definition of $T$ in (17)

\displaystyle\left|\langle Th,h\rangle_{L^{2}(D)}\right|

\displaystyle=

\displaystyle\left|\int_{D}\frac{e^{i\theta}f(y){\psi_{1}(x^{*}_{\perp})}|h(y)|^{2}}{2}\,\mbox{d}y\right|,

then we can directly prove the theorem by noting that $e^{i\theta}f>c_{1}>0$ or $e^{i\theta}f<c_{2}<0$ (and note that $\psi_{1}(x^{*}_{\perp})\not=0$ as the measurements are non-vanishing). This completes the proof. ∎

Lemma 2.

The operator $S^{*}:L^{2}(D)\to L^{2}(k_{-},k_{+})$ given by (18) is compact and has dense range in $L^{2}(k_{-},k_{+})$ .

Proof.

From the definition (18) of $S^{*}$ , for any $h\in L^{2}(D)$ we have that

(S^{*}h)(\sigma):=\int_{D}\overline{\sqrt{\psi_{1}(y_{\perp})}}e^{i\sigma(y_{1}-x^{*}_{1})}h(y)\,\mbox{d}y,\quad\sigma\in(k_{-},k_{+}),

and since $h\in L^{2}(D)$ , one can directly derive that $(S^{*}h)(\sigma)\in H^{1}(k_{-},k_{+})$ by taking the derivative of $(S^{*}h)(\sigma)$ with respect to $\sigma$ . Then we can conclude the compactness of $S^{*}:L^{2}(D)\to L^{2}(k_{-},k_{+})$ from the compact embedding from $H^{1}(k_{-},k_{+})$ into $L^{2}(k_{-},k_{+})$ .

To show the dense range of $S^{*}$ in $L^{2}(k_{-},k_{+})$ , it is sufficient to show that $S$ is injective. We prove this by contradiction. Indeed if there exists $\varphi\in L^{2}(k_{-},k_{+})$ such that $S\varphi=0$ , i.e.,

\int_{k_{-}}^{k_{+}}\sqrt{\psi_{1}(y_{\perp})}e^{-i\sigma(y_{1}-x^{*}_{1})}\varphi(\sigma)\,\mbox{d}\sigma=0,\quad\forall y\in D,

note that $\psi_{1}\not=0$ , then we can derive that

\mathcal{F}\{\varphi\chi_{(k_{-},k_{+})}\}(y_{1}-x^{*}_{1})=0,\quad\forall y_{1}\in D_{1}^{+},

where $\mathcal{F}$ is the Fourier transform $\mathcal{F}:\varphi\chi_{(k_{-},k_{+})}\to\int_{\mathbb{R}}\varphi(y)\chi_{(k_{-},k_{+})}(y)e^{-isy}\,\mbox{d}y$ , $D^{+}_{1}:=\{y_{1}:\exists\,y_{\perp}^{*}\mbox{ s.t. }(y_{1},y_{\perp}^{*})\in D\}$ and $\chi_{(k_{-},k_{+})}$ is the characteristic function such that $\chi_{(k_{-},k_{+})}(x)=1$ if $x\in(k_{-},k_{+})$ and $\chi_{(k_{-},k_{+})}(x)=0$ if $x\not\in(k_{-},k_{+})$ . This implies that $\varphi\chi_{(k_{-},k_{+})}$ vanishes and thereby $\varphi$ vanishes in $L^{2}(k_{-},k_{+})$ . This proves that $S$ is injective and hence $S^{*}$ has dense range in $L^{2}(k_{-},k_{+})$ . This completes the proof. ∎

4. The sampling methods

In this section, we study the sampling methods [13, 14, 15, 19, 24] where the goal is to lay the theoretical foundation for extended sources and design efficient imaging algorithms with the single mode multi-frequency measurements. The factorization method [23] has established a mathematically rigorous characterization of extended scatterers in the classical setting and has been applied to the inverse source problem in the homogeneous whole space [21] with sparse multi-frequency measurements. There have been relevant work on the factorization method in waveguide with a single frequency, yet we aim to deploy the factorization method with a single propagating mode at multiple frequencies. Based on the multi-frequency operator and its properties derived in Section 3, the factorization method is to be established in Section 4.1. Moreover, the factorization of the multi-frequency operator may also lead to a factorization-based sampling method, as demonstrated in waveguide imaging in the single frequency case [6] and the inverse source problem in the homogeneous space [25] with sparse multi-frequency measurements. In Section 4.2 the factorization-based sampling method is to be analyzed.

4.1. The factorization method

In this section, we shall investigate the factorization method which gives a mathematically rigorous characterization of the range support of the source. To begin with, denote by

\psi_{z}^{\epsilon}(\sigma):=\frac{1}{|B(z,\epsilon)|}\int_{B(z,\epsilon)}e^{i\sigma(y_{1}-x^{*}_{1})}\,\mbox{d}y,

(20)

where $B(z,\epsilon)$ is a circle centered at $z$ with radius $\epsilon$ , and $|B(z,\epsilon)|$ denotes its area.

We first give a characterization of the range of $S^{*}$ . Let $D^{+}:=\{(y_{1},y_{\perp})\in W:\exists\,y_{\perp}^{*}\mbox{ s.t. }(y_{1},y_{\perp}^{*})\in D\}$ , and recall that $D^{+}_{1}=\{y_{1}:\exists\,y_{\perp}^{*}\mbox{ s.t. }(y_{1},y_{\perp}^{*})\in D\}$ . For any domain $\Omega$ , denote by $\chi_{\Omega}$ the characteristic function such that $\chi_{\Omega}(x)=1$ if $x\in\Omega$ and $\chi_{\Omega}(x)=0$ if $x\not\in\Omega$ . For any point $z=(z_{1},z_{\perp})\in W$ , we call $z_{1}$ the range coordinate and $z_{\perp}$ the cross-section coordinate.

Lemma 3.

Let $S^{*}:L^{2}(D)\to L^{2}(k_{-},k_{+})$ be given by (18) and $\psi_{z}^{\epsilon}$ be given by (20) respectively. For any sampling point in $W$ with range coordinate $z_{1}$ , it holds that

(1)

If $z_{1}\in D^{+}_{1}$ , then there exists a positive $\epsilon$ depending on $z_{1}$ such that for any $z$ with range coordinate $z_{1}$ , $\psi_{z}^{\epsilon}\in Range(S^{*})$ .
(2)

If $z_{1}\not\in D^{+}_{1}$ , then for any sufficiently small $\epsilon>0$ with $[z_{1}-\epsilon,z_{1}+\epsilon]\cap\overline{D_{1}^{+}}=\emptyset$ and any $z$ with range coordinate $z_{1}$ , $\psi_{z}^{\epsilon}\not\in Range(S^{*})$ .

Proof.

(1). If ${z_{1}}\in D^{+}_{1}$ , then there must exist $\epsilon>0$ and $z^{*}=(z_{1},z^{*}_{\perp})\in D$ such that $B(z^{*},\epsilon)\in D$ . Let $h(y)=\frac{\chi_{B(z^{*},\epsilon)}(y)}{|B(z^{*},\epsilon)|}\frac{1}{\overline{\sqrt{\psi_{1}(y_{\perp})}}}$ , then $h\in L^{2}(D)$ and it is directly verified that

(S^{*}h)(\sigma)=\int_{D}\frac{\chi_{B(z^{*},\epsilon)}(y)}{|B(z^{*},\epsilon)|}e^{i\sigma(y_{1}-x^{*}_{1})}\,\mbox{d}y=\psi_{z^{*}}^{\epsilon}(\sigma).

This proves that $\psi_{z^{*}}^{\epsilon}\in Range(S^{*})$ . Now a direct calculation

	$\displaystyle\psi_{z^{*}}^{\epsilon}(\sigma)$	$\displaystyle=$	$\displaystyle\frac{1}{\|B(z^{},\epsilon)\|}\int_{B(z^{},\epsilon)}e^{i\sigma(y_{1}-x^{}_{1})}\,\mbox{d}y=\frac{1}{\epsilon^{2}}\int_{z_{1}-\epsilon}^{z_{1}+\epsilon}\int_{z^{}_{\perp}-\sqrt{\epsilon^{2}-(y_{1}-z_{1})^{2}}}^{z^{}_{\perp}+\sqrt{\epsilon^{2}-(y_{1}-z_{1})^{2}}}e^{i\sigma(y_{1}-x^{}_{1})}\,\mbox{d}y_{\perp}\,\mbox{d}y_{1}$
		$\displaystyle=$	$\displaystyle\frac{1}{\epsilon^{2}}\int_{z_{1}-\epsilon}^{z_{1}+\epsilon}2\sqrt{\epsilon^{2}-(y_{1}-z_{1})^{2}}e^{i\sigma(y_{1}-x^{*}_{1})}\,\mbox{d}y_{1}$

yields that $\psi_{z^{*}}^{\epsilon}$ is indeed independent of the cross-section coordinate $z^{*}_{\perp}$ , i.e. for any $z=(z_{1},z_{\perp})$ and $z^{*}=(z_{1},z^{*}_{\perp})$ , $\psi_{z}^{\epsilon}=\psi_{z^{*}}^{\epsilon}\in Range(S^{*})$ . This proves the first part.

(2). If $z_{1}\not\in{D_{1}^{+}}$ , for any sufficiently small $\epsilon>0$ with $[z_{1}-\epsilon,z_{1}+\epsilon]\cap\overline{D_{1}^{+}}=\emptyset$ , assume on the contrary that $\psi_{z}^{\epsilon}\in Range(S^{*})$ for a $z=(z_{1},z_{\perp})$ . Then there exists a $g\in L^{2}(D)$ such that

\int_{D}\overline{\sqrt{\psi_{1}(y_{\perp})}}e^{i\sigma(y_{1}-x^{*}_{1})}g(y)\,\mbox{d}y=(S^{*}g)(\sigma)=\psi_{z}^{\epsilon}(\sigma)=\frac{1}{|B(z,\epsilon)|}\int_{B(z,\epsilon)}e^{i\sigma(y_{1}-x^{*}_{1})}\,\mbox{d}y,

for any $\sigma\in(k_{-},k_{+})$ . Noting the common term $e^{i\sigma x^{*}_{1}}$ , the above equation is reduced to

\int_{\mathbb{R}^{d}}e^{i\sigma y_{1}}\chi_{D}(y)\widetilde{g}(y)\,\mbox{d}y=\int_{\mathbb{R}^{d}}e^{i\sigma y_{1}}\frac{\chi_{B(z,\epsilon)}(y)}{|B(z,\epsilon)|}\,\mbox{d}y,\qquad\forall\sigma\in(k_{-},k_{+})

(21)

where $\widetilde{g}(y):=\overline{\sqrt{\psi_{1}(y_{\perp})}}g(y)$ in $D$ .

By recalling the multi-dimensional Fourier transform $\mathcal{F}:L^{2}(\mathbb{R}^{d})\to L^{2}(\mathbb{R}^{d})$

\mathcal{F}\{u\}(s)=\int_{\mathbb{R}^{d}}u(y)e^{-is\cdot y}\,\mbox{d}y,\quad s\in\mathbb{R}^{d},\quad\forall u\in L^{2}(\mathbb{R}^{d}),

(22)

and the Radon transform $\mathcal{R}^{e}:L^{2}(\mathbb{R}^{d})\to L^{1}(\mathbb{R})$ for a unit direction $e\in\mathbb{R}^{d}$ that

\mathcal{R}^{e}\{u\}(t^{\prime})=\int_{e^{\perp}}u(t^{\prime}e+{y^{\prime}})\,\mbox{d}y^{\prime},\quad t^{\prime}\in\mathbb{R},\quad\forall u\in L^{2}(\mathbb{R}^{d}),

(23)

where $e^{\perp}$ represents the hyperplane orthogonal to $e$ . From the Fourier slice theorem [30, Theorem 1.1]

\mathcal{F}\{u\}(te)=\widehat{\mathcal{R}^{e}}\{u\}(t),

(24)

where $~{}\widehat{}~{}$ denotes the one-dimensional Fourier transform.

Now we can derive from (21) that

\mathcal{F}\left\{\chi_{D}\widetilde{g}\right\}(-\sigma e_{1})=\mathcal{F}\left\{\frac{\chi_{B(z,\epsilon)}}{|B(z,\epsilon)|}\right\}(-\sigma e_{1}),\quad e_{1}=(1,0,\cdots),\qquad\forall\sigma\in(k_{-},k_{+}).

and thereby

\widehat{\mathcal{R}^{e_{1}}}\left\{\chi_{D}\widetilde{g}\right\}(-\sigma)=\widehat{\mathcal{R}^{e_{1}}}\left\{\frac{\chi_{B(z,\epsilon)}}{|B(z,\epsilon)|}\right\}(-\sigma),\qquad\forall\sigma\in(k_{-},k_{+}).

(25)

Then we can derive that

\mathcal{R}^{e_{1}}\left\{\chi_{D}\widetilde{g}\right\}=\mathcal{R}^{e_{1}}\left\{\frac{\chi_{B(z,\epsilon)}}{|B(z,\epsilon)|}\right\}.

However $~{}\mbox{supp}~{}\mathcal{R}^{e_{1}}\left\{\chi_{D}\widetilde{g}\right\}\subset e_{1}\cdot D=D^{+}_{1}$ and $~{}\mbox{supp}~{}\mathcal{R}^{e_{1}}\left\{\frac{\chi_{B(z,\epsilon)}}{|B(z,\epsilon)|}\right\}\subset e_{1}\cdot B(z,\epsilon)=(z_{1}-\epsilon,z_{1}+\epsilon)$ yields a contradiction (since $[z_{1}-\epsilon,z_{1}+\epsilon]\cap\overline{D_{1}^{+}}=\emptyset$ ). This prove the second part. These prove the Lemma. ∎

As can be seen from Lemma 3 part $(1)$ , the $\epsilon$ depends on the sampling range coordinate $z_{1}$ . The following corollary removes this dependence with an extra geometric assumption on $D$ .

Corollary 1.

Assume in addition that the minimal width of $D$ is bounded below by a positive constant, i.e., $\min\{L(y_{1}):\forall y_{1}\}>c>0$ where $L(y_{1}):={\min_{y_{\perp}}\max_{y_{\perp}^{*}}}\{|y_{\perp}-y_{\perp}^{*}|:(y_{1},\beta y_{\perp}+(1-\beta)y_{\perp}^{*})\in D,\forall\beta\in[0,1]\}$ . For any sampling point in $W$ with range coordinate $z_{1}$ , it holds that

(1)

If $z_{1}\in D^{+}_{1}$ , let $\epsilon_{1}>0$ be any given small tolerance level, then there exists a positive $\epsilon_{0}$ such that for any $z$ with $(z_{1}-\epsilon_{1},z_{1}+\epsilon_{1})\subset D_{1}^{+}$ and any positive $\epsilon<\epsilon_{0}$ , $\psi_{z}^{\epsilon}\in Range(S^{*})$ .
(2)

If $z_{1}\not\in D_{1}^{+}$ , then for any sufficiently small $\epsilon>0$ with $[z_{1}-\epsilon,z_{1}+\epsilon]\cap\overline{D_{1}^{+}}=\emptyset$ and any $z$ with range coordinate $z_{1}$ , $\psi_{z}^{\epsilon}\not\in Range(S^{*})$ .

Proof.

The only difference is in part (1). We highlight the difference. Assume that $z_{1}\in D^{+}_{1}$ . Since the minimal width of $D$ is bounded below by a positive constant, i.e., ${\min_{y_{\perp}}\max_{y_{\perp}^{*}}}\{|y_{\perp}-y_{\perp}^{*}|:(y_{1},\beta y_{\perp}+(1-\beta)y_{\perp}^{*})\in D,\forall\beta\in[0,1]\}>c>0$ , then for any range coordinate $z_{1}\in D^{+}_{1}$ with $(z_{1}-\epsilon_{1},z_{1}+\epsilon_{1})\subset D$ , there must exist $z^{*}=(z_{1},z^{*}_{\perp})\in D$ such that $B(z^{*},\epsilon)\in D$ for all $0<\epsilon<{\min\{c/4,\epsilon_{1}/4\}}$ . Now let $h(y)=\frac{\chi_{B(z^{*},\epsilon)}(y)}{|B(z^{*},\epsilon)|}\frac{1}{\overline{\sqrt{\psi_{1}(y_{\perp})}}}$ , then $h\in L^{2}(D)$ and it is directly verified that

(S^{*}h)(\sigma)=\int_{D}\frac{\chi_{B(z^{*},\epsilon)}(y)}{|B(z^{*},\epsilon)|}e^{i\sigma(y_{1}-x^{*}_{1})}\,\mbox{d}y=\psi_{z^{*}}^{\epsilon}(\sigma).

This proves that $\psi_{z^{*}}^{\epsilon}\in Range(S^{*})$ . Using the same argument before, we prove the corollary. ∎

Applying Theorem 2 and Lemma 2, the following Lemma is a direct application of [24, Corollary 1.22].

Lemma 4.

Let $F:L^{2}(k_{-},k_{+})\to L^{2}(k_{-},k_{+})$ be given by (15) and $S^{*}:L^{2}(D)\to L^{2}(k_{-},k_{+})$ be given by (18) respectively. It holds that $Range(S^{*})=Range(F^{1/2})$ .

Now we can prove the following theorem. Let $(\phi_{j},\alpha_{j})_{j=1}^{\infty}$ be the eigensystem of the self-adjoint operator $F^{1/2}$ and let $\langle\cdot,\cdot\rangle$ denote the inner product in $L^{2}(k_{-},k_{+})$ .

Theorem 3.

For any sampling point in $W$ with range coordinate $z_{1}$ , it holds that

(1)

If $z_{1}\in D^{+}_{1}$ , then there exists a positive $\epsilon$ depending on $z_{1}$ such that for any $z$ with range coordinate $z_{1}$ , $\psi_{z}^{\epsilon}\in Range(F^{1/2})$ and moreover

$\sum_{j=1}^{\infty}\frac{|\langle\psi_{z}^{\epsilon},\phi_{j}\rangle|^{2}}{\alpha_{j}^{2}}<\infty.$ (26)
(2)

If $z_{1}\not\in D^{+}_{1}$ , then for any sufficiently small $\epsilon>0$ with $[z_{1}-\epsilon,z_{1}+\epsilon]\cap\overline{D_{1}^{+}}=\emptyset$ and any $z$ with range coordinate $z_{1}$ , $\psi_{z}^{\epsilon}\not\in Range(F^{1/2})$ and moreover

$\sum_{j=1}^{\infty}\frac{|\langle{\psi_{z}^{\epsilon}},\phi_{j}\rangle|^{2}}{\alpha_{j}^{2}}=\infty.$ (27)

Proof.

Since $Range(S^{*})=Range(F^{1/2})$ . Then the proof is completed by applying Lemma 3 and Picard theorem [19, Theorem 4.8]. ∎

Corollary 2.

Assume in addition that the minimal width of $D$ is bounded below by a positive constant, i.e., $\min\{L(y_{1}):\forall y_{1}\}>c>0$ where $L(y_{1})={\min_{y_{\perp}}\max_{y_{\perp}^{*}}}\{|y_{\perp}-y_{\perp}^{*}|:(y_{1},\beta y_{\perp}+(1-\beta)y_{\perp}^{*})\in D,\forall\beta\in[0,1]\}$ . Then we have that

(1)

If $z_{1}\in D^{+}_{1}$ , let $\epsilon_{1}>0$ be any given small tolerance level, then there exists a positive $\epsilon_{0}$ such that for any $z$ with $(z_{1}-\epsilon_{1},z_{1}+\epsilon_{1})\subset D_{1}^{+}$ and any positive $\epsilon<\epsilon_{0}$ , $\psi_{z}^{\epsilon}\in Range(F^{1/2})$ and equation (26) holds.
(2)

If $z_{1}\not\in D^{+}_{1}$ , then for any sufficiently small $\epsilon>0$ with $[z_{1}-\epsilon,z_{1}+\epsilon]\cap\overline{D_{1}^{+}}=\emptyset$ and any $z$ with range coordinate $z_{1}$ , $\psi_{z}^{\epsilon}\not\in Range(F^{1/2})$ and equation (27) holds.

Remark 1.

Corollary 2 implies that the $\epsilon$ can be chosen as a fixed constant for all the sampling points $z$ when implementing the factorization method. In particular, we can choose $\epsilon$ as a small number such that $\epsilon$ is less than a tolerance level that one is seeking for. If there is no such assumption as in Corollary 2, we can still chose $\epsilon$ as a fixed tolerance level, and in this case, the part whose width is very small may not be reconstructed as indicated by Theorem 3 and Lemma 3 (and its proof). However since $\epsilon$ is the tolerance level, the factorization method theoretically would still yield a sufficiently good image.

4.2. The factorization-based sampling method

In this section, we investigate a factorization-based multi-frequency sampling method with imaging function given by

I(z):=|\langle F\psi_{z},\psi_{z}\rangle|,\quad\psi_{z}(\sigma)=e^{i\sigma(z_{1}-x^{*}_{1})}.

(28)

Its property is given by the following theorem.

Theorem 4.

There exists positive constants $c_{1}$ and $c_{2}$ such that

c_{1}\|S\psi_{z}\|^{2}_{L^{2}(D)}\leq|I(z)|\leq c_{2}\|S\psi_{z}\|^{2}_{L^{2}(D)}.

(29)

Proof.

From the factorization of the far-field operator, it first follows that

\displaystyle I(z)=|\langle F\psi_{z},\psi_{z}\rangle_{L^{2}(k_{-},k_{+})}|=|\langle S^{*}TS\psi_{z},\psi_{z}\rangle_{L^{2}(k_{-},k_{+})}|=|\langle TS\psi_{z},S\psi_{z}\rangle_{L^{2}(D)}|.

The lower bound then follows from the coercivity of $T$ in Theorem 2, the upper bound follows directly from the boundedness of $T$ . This completes the proof of the theorem. ∎

Theorem 4 implies that $I(z)$ is qualitatively the same as $\|S\psi_{z}\|^{2}_{L^{2}(D)}$ without the knowledge of $D$ . The following Lemma gives the explicit expression of $S\psi_{z}$ .

Lemma 5.

It holds that

(S\psi_{z})(y)=\sqrt{\psi_{1}(y_{\perp})}\frac{e^{ik_{+}(z_{1}-y_{1})}-e^{ik_{-}(z_{1}-y_{1})}}{i(z_{1}-y_{1})},\quad y\in D,

(30)

where the above quotient is understood in the limit sense when $z_{1}=y_{1}$ .

Proof.

From the definition of $S$ , we directly have that

	$\displaystyle(S\psi_{z})(y)$	$\displaystyle=$	$\displaystyle\sqrt{\psi_{1}(y_{\perp})}\int_{k_{-}}^{k_{+}}e^{i\sigma(z_{1}-y_{1})}\,\mbox{d}\sigma$
		$\displaystyle=$	$\displaystyle\sqrt{\psi_{1}(y_{\perp})}\frac{e^{ik_{+}(z_{1}-y_{1})}-e^{ik_{-}(z_{1}-y_{1})}}{i(z_{1}-y_{1})},$

where the above quotient is understood in the limit sense when $z_{1}=y_{1}$ . This completes the proof. ∎

From the explicit formula of $S\psi_{z}$ in Lemma 5, we find that its maximum is obtained at $z_{1}=y_{1}$ and that $S\psi_{z}$ is oscillatory and decaying to zero when $|z_{1}-y_{1}|$ becomes large. We illustrate this property in Figure 3 with a fixed $y_{1}=0$ and $y_{\perp}$ , where we plot $|S\psi_{z}|/\max\{|S\psi_{z}|\}$ for a two-dimensional waveguide $(-\infty,\infty)\times(0,\pi)$ with Dirichlet boundary condition where $k_{-}=0$ and $k_{+}=\sqrt{3}$ . We further remark that the function $S\psi_{z}$ is only sensitive to the range sampling coordinate $z_{1}$ which can be directly observed from (30). As a summary, it is heuristically observed that $\|S\psi_{z}\|^{2}_{L^{2}(D)}$ peak in $D^{+}=\{(y_{1},y_{\perp})\in W:\exists\,y_{\perp}^{*}\mbox{ s.t. }(y_{1},y_{\perp}^{*})\in D\}$ as the sampling point $z$ samples in a sampling region since there always exists $y=(y_{1},y_{\perp})\in D$ such that $y_{1}=z_{1}$ for any $z=(z_{1},z_{\perp})\in D^{+}$ . Together with Theorem 4 it is inferred that $I(z)$ is qualitatively the same as $\|S\psi_{z}\|^{2}_{L^{2}(D)}$ without the knowledge of $D$ . Though the factorization-based sampling method is not as mathematically rigorous as the factorization method, it is heuristically inferred that $I(z)$ peaks in $D^{+}$ which is the range support of $D$ in the waveguide.

5. Numerical examples

In this section, we present a variety of numerical examples in two dimensions to illustrate the performance of the single mode multi-frequency sampling methods.

We generate the synthetic data $u^{s}$ using the finite element computational software Netgen/NGSolve [32]. To be more precise, the computational domain is ${(-4,4)}\times\Sigma$ and the measurements are at a random location on the cross section ${\{-2\}}\times\Sigma$ . We apply a brick Perfectly Matched Layer (PML) in the domain $\{x_{1}:2.5<|x_{1}|<4\}\times\Sigma$ and choose PML absorbing coefficient $5+5i$ . In all of the numerical examples, we apply the second order finite element to solve for the wave field where the source is constant $1$ ; the mesh size is chosen as $|\Sigma|/10$ where we recall that $|\Sigma|$ denotes the length of the cross section. We further add $5\%$ Gaussian noise to the synthetic data $u^{s}$ to implement the imaging function in Matlab.

We first illustrate the implementation of the factorization method and the factorization-based sampling method. Let $\{k_{1},k_{2},\cdots,k_{N}\}$ be the set of equally discretized wave numbers, this yields the discretized $N\times N$ multi-frequency far-field matrix $F_{N}$ of its continuous counterpart $F$ in (15).

(1)

The factorization method. Let $(\widehat{\phi}_{j},\widehat{\alpha}_{j})_{j=1}^{N}$ be the eigensystem of the self-adjoint matrix $(F_{N})^{1/2}$ . The imaging function indicated by the factorization method is

I_{\mbox{\tiny FM}}(z):=\left(\sum_{j=1}^{N}\frac{|\langle\widehat{\psi}_{z}^{\epsilon},\widehat{\phi}_{j}\rangle|^{2}}{\widehat{\alpha}_{j}^{2}}\right)^{-1},

where $\epsilon$ is a fixed tolerance level, and $\widehat{\psi}_{z}^{\epsilon}$ is the discretization of ${\psi}_{z}^{\epsilon}$ (20) at the multiple wavenumbers (i.e. multi-frequencies) $\{k_{1},k_{2},\cdots,k_{N}\}$ . The $\langle\cdot,\cdot\rangle$ represents the inner product of two vectors that arise from the discretization of their continuous counterparts. The main Theorem 2 of the factorization method implies that $I_{FM}(z)$ is bounded below by $0$ for $z\in D^{+}_{1,\epsilon}$ and almost vanishes for $z_{1}\not\in D^{+}_{1,\epsilon}$ where $D^{+}_{1,\epsilon}$ is a neighborhood of the support $D_{1}^{+}$ within the tolerance level $\epsilon$ . We choose tolerance level $\epsilon=0.01$ in all of the numerical examples. Though we have fixed this tolerance level, we have tested other small tolerance levels and haven’t found any visual difference on the imaging results.

(2)

The factorization-based sampling method. Let $\widehat{\psi}_{z}$ be the discretization of ${\psi}_{z}$ at the multiple wavenumbers (i.e. multi-frequencies) $\{k_{1},k_{2},\cdots,k_{N}\}$ . The imaging function indicated by the factorization-based sampling method is

$I_{\mbox{\tiny FBSM}}(z):=|\langle F_{N}\widehat{\psi}_{z},\widehat{\psi}_{z}\rangle|.$

The $\langle\cdot,\cdot\rangle$ represents the inner product of two vectors that arise from the discretization of their continuous counterparts.

For the visualization, we plot both $I_{\mbox{\tiny FM}}(z)$ and $I_{\mbox{\tiny FBSM}}(z)$ over the sampling region $\{-2,2\}\times\Sigma$ and we always normalize them such that their maximum value are $1$ respectively.

5.1. Number of frequencies

The first set of numerical examples is to illustrate how the number of frequencies affects the image. We consider a rectangular source in a two dimensional waveguide $\{-\infty,\infty\}\times(0,h)$ with Neumann boundary condition and $h=\pi/12$ . In this case, the first propagating mode is $e^{ik|x_{1}|}\frac{1}{\sqrt{h}}$ . The range reconstruction of the rectangular source is shown in Figure 4 where we have considered three different cases:

•

Case I: 47 frequencies in the set $\{0.25,0.5,\cdots,11.5,11.75\}$ .
•

Case II: 23 frequencies in the set $\{0.5,1.0,\cdots,11.0,11.5\}$ .
•

Case III: 11 frequencies in the set $\{1.0,2.0,\cdots,10.0,11.0\}$ .

It is observed that the factorization method performs better with more frequencies, while the factorization-based sampling method performs almost the same in these three different cases. Such a difference may be due to that the factorization method explicitly involves solving an ill-posed equation so that a more sophisticated discretization is perhaps needed in the implementation. Another observation is that the imaging function of the factorization method sharply vanish outside the range support (which agrees with the factorization method main Theorem 3), and the imaging function of the factorization-based sampling method gradually becomes small as the sampling point samples away from the range support (which also agrees with the main theorem of the factorization-based sampling method Theorem 4 and the property of the point spread function $S\psi_{z}$ ).

5.2. A “complete block” source

Using sampling methods to image of a complete block of the waveguide is an interesting problem [27]. The L-shape case under consideration in this subsection is similar to a complete block case considered in [27]. Results in Figure 5 imply that a “complete block” source (i.e. a source which spans the entire cross-section) can be reconstructed with our sampling methods with multi-frequency measurements. This is due to that the multi-frequency measurements make use of the phase (travel time) which provides information on the bulk location.

5.3. Different boundary conditions

We further illustrate the performance of the sampling methods with different boundary conditions of the waveguide. In this regard, a two dimensional waveguide $\{-\infty,\infty\}\times(0,h)$ with a mixed boundary condition and $h=\pi/12$ is considered. The mixed boundary condition is a Dirichlet boundary condition on the top boundary and a Neumann boundary condition on the bottom boundary. The first propagating mode in this case is $e^{i\sqrt{k^{2}-(\pi/2h)^{2}}|x_{1}|}\sqrt{\frac{2}{h}}\cos(\pi x_{\perp}/2h)$ , and the corresponding dispersion relation reads $\mu_{1}(k)=\sqrt{k^{2}-(\pi/2h)^{2}}$ . Note that this mixed boundary condition yields a non-linear dispersion relation. The first passband corresponds to $k\in(6,12)$ , and the following numerical example is produced with 41 wavenumbers (i.e. frequencies) in the set $\{k=\sqrt{(\pi/2h)^{2}+\sigma^{2}}:\sigma=0.25,0.5,\cdots,10.0,10.25\}$ .

The first row of Figure 6 is the range reconstruction of a rectangular source. The support of the source is the same as the one in Figure 4 in order to compare performance with respect to different boundary conditions, and it is observed that the performances of the sampling methods in these two different boundary condition cases are similar. The difficulty arising from the non-linear dispersion relation in the mixed boundary condition case is overcame via the carefully defined far-field operator.

The second row of Figure 6 is the range reconstruction of a rhombus shape source. Note that the rhombus shape source does not satisfy the assumption in Corollary 2. In particular the left and right corner of the rhombus shape source would be reconstructed with certain tolerance for the factorization method as indicated by Theorem 3 and Lemma 3 (and its proof). Since the left and right corners are sharp, the factorization-based sampling method may also reconstruct the support with certain tolerance. These observations are illustrated by the the second row of Figure 6.

5.4. A complete sound-soft block

In the last subsection, we illustrate how to image a complete block (with sound soft boundary condition) by a preliminary numerical example. It is observed in [27] that, a complete block of a waveguide can be reconstructed with near field measurements, however it may not be reconstructed with far-field measurements using the linear sampling method at a single frequency. With the multi-frequency measurements in our framework, it is possible to image such a complete block using far-field measurements. Though we have demonstrated this via the inverse source problem, to be more convincing we show a preliminary example for imaging a sound soft block. The waveguide is $(-\infty,\infty)\times(0,h)$ with Neumann boundary condition and $h=\pi/12$ . The complete block is at $-0.5\times(0,h)$ . In this case, the scattered wave field at far-field $x_{r}\in\{-10\}\times(0,2)$ due to a point source at transmitter location $x_{s}\in\{-10\}\times(0,2)$ can be written explicitly and we further add $5\%$ noise to this wave field to generate our synthetic data.

The difference between the block case and the source case can be interpreted though the travel time: the travel time in the block case roughly doubles the one in the source case (where the block and the source are at the same location); this is because the wave travel from the transmitter to the block then back to the receiver in the block case while the wave just travel from the source to the receiver in the source case. Therefore the functions $\psi_{z}^{\epsilon}$ and $\psi_{z}$ have to be modified by doubling their phases, i.e. $\psi_{z}^{\epsilon}(\sigma)\to\frac{1}{|B(z,\epsilon)|}\int_{B(z,\epsilon)}e^{i2\sigma(y_{1}-x_{1}^{*})}\,\mbox{d}y$ and $\psi_{z}(\sigma)\to e^{i2\sigma(z_{1}-x_{1}^{*})}\sqrt{\psi_{1}(z_{\perp})}$ . This is the only change in our implementation to image a complete block with the sampling methods. Figure 7 illustrate the performance of the sampling methods to image a complete sound soft block. $47$ frequencies in the set $\{0.25,0.5,\cdots,11.5,11.75\}$ are used. It is observed that the block is correctly located and a sharper result is given by the factorization method. The obstacle/block case is more complicated than the source case and a complete theoretical justification is part of our future work.

6. Remark on designs of multi-frequency operator and assumptions on the source

In the above discussions of the multi-frequency sampling methods, we have made the assumption on the source $f$ that $e^{i\theta}f>c_{1}>0$ or $e^{i\theta}f<c_{2}<0$ with some $\theta\in[0,2\pi)$ . Indeed it is possible to relax this assumption by designing other multi-frequency operators. We introduce in this section another two multi-frequency operators to remark that the sampling methods can require less restrictive conditions on the source; this may be of interest to readers who is more interested in minimal assumption on the source. The first operator in Section 6.1 is to consider both “backward and forward scattering” measurements; the second operator in Section 6.2 still uses only the “backscattering” measurements (which is less obvious to design) and it usually uses information in part of the frequency range. Their advantages are that the sampling methods with these operators would require a less restrictive assumption that $\Re(f(x)e^{i\tau})>c_{3}>0$ almost everywhere for $x\in D$ with some constant phase $\tau$ ; their drawbacks are that either more measurements have to be used (Section 6.1) or information in part of the frequency range can usually be utilized (Section 6.2). The first technique in Section 6.1 is in the same spirit of the inverse source problem in the whole space [21] (that uses measurements from two opposite directions); the second technique in Section 6.2 seems new which may facilitate further advancement of the multi-frequency factorization method using only the backscattering measurements. In the following we briefly show how to generalize the results in Section 3–4 using these two multi-frequency operators.

6.1. A multi-frequency operator based on “backward and forward scattering” measurements

One way to relax the assumption on the source is to consider both the “backward and forward scattering” measurements, i.e. measurements at $x_{l}:=(x^{*}_{1},x^{*}_{\perp})$ and $x_{r}:=(-x^{*}_{1},x^{*}_{\perp})$ . In this case, we introduce a similar multi-frequency far-field operator $\widetilde{F}:L^{2}(k_{-},k_{+})\to L^{2}(k_{-},k_{+})$ by

(\widetilde{F}g)(\sigma;x_{l},x_{r}):=\int_{k_{-}}^{k_{+}}{-i}\mu_{1}(\omega_{\sigma\gamma})\Big{[}H(\sigma-\gamma)u_{p}^{s}\left({x_{l}};\omega_{\sigma\gamma}\right){+}H(\gamma-\sigma)e^{2i{|\sigma-\gamma|}x^{*}_{1}}{u_{p}^{s}\left({x_{r}};\omega_{\sigma\gamma}\right)}\Big{]}g(\gamma)\,\mbox{d}\gamma,

(31)

$\sigma\in(k_{-},k_{+})$ . One can then verify that

	$\displaystyle(\widetilde{F}g)(\sigma;x_{l},x_{r})$	$\displaystyle=$	$\displaystyle 1/2\int_{k_{-}}^{k_{+}}\int_{D}\psi_{1}(y_{\perp})e^{i(\sigma-\gamma)\|x^{}_{1}-y_{1}\|}f(y){\psi_{1}(x^{}_{\perp})}g(\gamma)\,\mbox{d}y\,\mbox{d}\gamma$
		$\displaystyle=$	$\displaystyle\int_{D}e^{i\sigma\|x^{}_{1}-y_{1}\|}\overline{\sqrt{\psi_{1}(y_{\perp})}}\left[\left(\int_{k_{-}}^{k_{+}}\sqrt{\psi_{1}(y_{\perp})}e^{-i\gamma\|x^{}_{1}-y_{1}\|}g(\gamma)\,\mbox{d}\gamma\right)\frac{f(y){\psi_{1}(x^{*}_{\perp})}}{2}\right]\,\mbox{d}y,$

and we can proceed as in Theorem 1 to conclude that $\widetilde{F}=S^{*}\widetilde{T}S$ with $\widetilde{T}:L^{2}(D)\to L^{2}(D)$

(\widetilde{T}h)(y):=\frac{\psi_{1}(x^{*}_{\perp})f(y)h(y)}{2},\quad y\in D.

(32)

Though there is no assumption on $f$ to derive the factorization $\widetilde{F}=S^{*}\widetilde{T}S$ in the “backward and forward scattering” measurement case, an assumption that $\Re(f(x)e^{i\tau})>{c_{3}>0}$ almost everywhere for $x\in D$ with some $\tau$ would then be added to ensure the coercivity of an operator related to $\widetilde{T}$ . To show the coercivity, we introduce for a generic bounded operator $A:L^{2}\to L^{2}$ that

\Re A=\frac{A+A^{*}}{2},\quad\Im A=\frac{A-A^{*}}{2i}.

Now we have

Lemma 6.

Suppose that $\Re(f(x)e^{i\tau})\geq c_{3}>0,\forall x\in D$ for some constant phase $\tau$ , then we have that $\Re(e^{i\tau}\widetilde{T})$ is self-adjoint and coercive.

Proof.

From the assumption that $\Re(fe^{i\tau})\geq c_{3}>0$ for some constant phase $\tau$ , we have that

\Big{(}\Re(e^{i\tau}\widetilde{T})h\Big{)}(y):=\frac{\psi_{1}(x^{*}_{\perp})\Re(fe^{i\tau})h(y)}{2},\quad y\in D,

(33)

and that

\displaystyle\langle\Re(e^{i\tau}\widetilde{T})h,h\rangle_{L^{2}(D)}

\displaystyle=

\displaystyle\int_{D}\frac{\Re(e^{i\tau}f(y)){\psi_{1}(x^{*}_{\perp})}|h(y)|^{2}}{2}\,\mbox{d}y\geq c^{\prime}\|h\|^{2}_{L^{2}(D)}

(34)

for some positive constant $c^{\prime}$ . This completes the proof. ∎

Lemma 6 would then allow us to derive a $F\#$ version of the factorization method (that requires less restrictive conditions on the source $f$ ) which is in the same spirit of the inverse source problem in the whole space [21] (that uses measurements from two opposite directions). However we are mainly concerned with the most interesting “backscattering” case and we shall discuss another possible multi-frequency operator $F_{\alpha}$ in the next section. We then finalize by generalizing the main results of the factorization method and the factorization-based sampling method using $\widetilde{F}$ and $F_{\alpha}$ .

6.2. Another multi-frequency operator based on “backscattering” measurements

Indeed we can still impose less restrictive conditions on the source $f$ in the case of “backscattering” by designing another multi-frequency operator as follows. Here we introduce a multi-frequency far-field operator $F_{\alpha}:L^{2}(k_{-},k_{+}(\alpha))\to L^{2}(k_{-},k_{+}(\alpha))$ by

(F_{\alpha}g)(\sigma;x^{*}):=\int_{k_{-}}^{k_{+}(\alpha)}-i\mu_{1}(\omega_{\sigma\gamma\alpha})u_{p}^{s}\left(x^{*};\omega_{\sigma\gamma\alpha}\right)g(\gamma)\,\mbox{d}\gamma,\quad\sigma\in(k_{-},k_{+}(\alpha))

(35)

where $\omega_{\sigma\gamma\alpha}:=\sqrt{\lambda_{1}^{2}+\big{(}\sigma-\gamma+\frac{1}{\alpha}\sqrt{\lambda_{2}^{2}-\lambda_{1}^{2}}\big{)}^{2}}$ with $\alpha\geq 2$ being a positive constant that is determined later (by the coercivity assumption of a certain operator), $k_{-}=0$ and $k_{+}(\alpha):=\frac{1}{\alpha}\sqrt{\lambda_{2}^{2}-\lambda_{1}^{2}}$ . Note that for $\sigma,\gamma\in(k_{-},k_{+}(\alpha))$ , $\omega_{\sigma\gamma\alpha}\in\Big{(}\lambda_{1},\sqrt{\lambda^{2}_{1}+\frac{4}{\alpha^{2}}(\lambda^{2}_{2}-\lambda^{2}_{1})}\Big{)}\subseteq(\lambda_{1},\lambda_{2})$ .

In this case, similar to the role of $\omega_{\sigma\gamma}$ in Section 3, one finds via a direct calculation that $\mu_{1}(\omega_{\sigma\gamma\alpha})=\sigma-\gamma+\frac{1}{\alpha}\sqrt{\lambda_{2}^{2}-\lambda_{1}^{2}}$ (for $\sigma,\gamma\in(k_{-},k_{+}(\alpha))$ ) which linearizes the dispersion in another (though less obvious) way. In principle, one may find other ways to design multi-frequency operators. The possible drawback is seen as that $\omega_{\sigma\gamma\alpha}$ has range in $\Big{(}\lambda_{1},\sqrt{\lambda^{2}_{1}+\frac{4}{\alpha^{2}}(\lambda^{2}_{2}-\lambda^{2}_{1})}\Big{)}$ which is a subset of $(\lambda_{1},\lambda_{2})$ if $\alpha>2$ , this implies that information in part of the frequency range is only used if $\alpha>2$ . Nevertheless, with such a design of the multi-frequency operator, we are able to derive a factorization method with only the “backscattering” measurements with less restrictive conditions on the source $f$ . We also refer to [21] for a discussion on an extension of the factorization method to “backscattering” measurements.

The operators $S_{\alpha}:L^{2}(k_{-},k_{+}(\alpha))\to L^{2}(D)$ and $S^{*}_{\alpha}:L^{2}(D)\to L^{2}(k_{-},k_{+}(\alpha))$ remain formally the same as in (16) and (18), respectively (by only replacing $k_{+}$ by $k_{+}(\alpha)$ ). The difference is the operator $T_{\alpha}:L^{2}(D)\to L^{2}(D)$ defined via

(T_{\alpha}h)(y):=\frac{e^{i\frac{1}{\alpha}\sqrt{\lambda_{2}^{2}-\lambda_{1}^{2}}|x_{1}^{*}-y_{1}|}\psi_{1}(x^{*}_{\perp})f(y)h(y)}{2},\quad y\in D.

(36)

Using the definition of $F_{\alpha}$ (35), we can prove similarly that

			$\displaystyle(F_{\alpha}g)(\sigma;x^{})=1/2\int_{k_{-}}^{k_{+}(\alpha)}\int_{D}\psi_{1}(y_{\perp}){e^{i\left(\sigma-\gamma+\frac{1}{\alpha}\sqrt{\lambda_{2}^{2}-\lambda_{1}^{2}}\right)\|x^{}_{1}-y_{1}\|}\psi_{1}(x^{*}_{\perp})}f(y)g(\gamma)\,\mbox{d}y\,\mbox{d}\gamma$
		$\displaystyle=$	$\displaystyle\int_{D}e^{i\sigma\|x^{}_{1}-y_{1}\|}\overline{\sqrt{\psi_{1}(y_{\perp})}}\Bigg{[}\left(\int_{k_{-}}^{k_{+}(\alpha)}\sqrt{\psi_{1}(y_{\perp})}e^{-i\gamma\|x^{}_{1}-y_{1}\|}g(\gamma)\,\mbox{d}\gamma\right)\frac{{e^{i\frac{1}{\alpha}\sqrt{\lambda_{2}^{2}-\lambda_{1}^{2}}\|x^{}_{1}-y_{1}\|}\psi_{1}(x^{}_{\perp})}f(y)}{2}\Bigg{]}\,\mbox{d}y$
		$\displaystyle=$	$\displaystyle\big{(}S_{\alpha}^{*}T_{\alpha}S_{\alpha}g\big{)}(\sigma).$

Next we investigate the coercivity of the operator $\Re(e^{i\tau}T_{\alpha})$ . Recall that $\widetilde{T}:L^{2}(D)\to L^{2}(D)$ is given by (32) in Section 6.1.

Theorem 5.

Suppose that $\Re(e^{i\tau}\widetilde{T})$ is positive and coercive for some constant phase $\tau$ . Let $T_{\alpha}$ be given by (36), then there exists a positive constant $\alpha\geq 2$ such that the operator $\Re(e^{i\tau}T_{\alpha}):L^{2}(D)\to L^{2}(D)$ is self-adjoint, coercive and

\langle\Re(e^{i\tau}T_{\alpha})h,h\rangle_{L^{2}(D)}\geq c\|h\|^{2}_{L^{2}(D)},

for some positive constant $c$ .

Proof.

We first introduce $\zeta:=\frac{1}{\alpha}\sqrt{\lambda_{2}^{2}-\lambda_{1}^{2}}|x_{1}^{*}-y_{1}|$ in the proof to avoid writing long equations so that the presentation is clean and compact. From the definitions of $T_{\alpha}$ in (36), we have that

\Big{(}(e^{i\tau}T_{\alpha})h\Big{)}(y)=\frac{e^{i\zeta}\psi_{1}(x^{*}_{\perp})e^{i\tau}f(y)h(y)}{2},\quad\Big{(}(e^{i\tau}T_{\alpha})^{*}h\Big{)}(y)=\frac{e^{-i\zeta}\psi_{1}(x^{*}_{\perp})\overline{e^{i\tau}f(y)}h(y)}{2},\quad y\in D,

whereby

\displaystyle\Big{(}\Re(e^{i\tau}T_{\alpha})h\Big{)}(y)=\Big{(}\cos(\zeta)\Re(e^{i\tau}f(y))-\sin(\zeta)\Im(e^{i\tau}f(y))\Big{)}\frac{\psi_{1}(x^{*}_{\perp})h(y)}{2},\quad y\in D.

(37)

Now from the expressions of $\Re(e^{i\tau}T_{\alpha})$ in (37) and of $\Re(e^{i\tau}\widetilde{T})$ in (33),

			$\displaystyle\left\|\langle\Re(e^{i\tau}T_{\alpha})h,h\rangle_{L^{2}(D)}-\langle\Re(e^{i\tau}\widetilde{T})h,h\rangle_{L^{2}(D)}\right\|$		(38)
		$\displaystyle=$	$\displaystyle\Big{\|}\int_{D}\Big{(}\cos(\zeta)\Re(e^{i\tau}f(y))-\Re(e^{i\tau}f(y))-\sin(\zeta)\Im(e^{i\tau}f(y))\Big{)}\frac{\psi_{1}(x^{*}_{\perp})\|h(y)\|^{2}}{2}\,\mbox{d}y\Big{\|}.\qquad$		(38)

Having $y\in D$ in a fixed sampling region, then $|\sqrt{\lambda_{2}^{2}-\lambda_{1}^{2}}(y_{1}-x_{1}^{*})|$ is bounded above so that we can always choose $\alpha$ such that $|\zeta|=|\frac{1}{\alpha}\sqrt{\lambda_{2}^{2}-\lambda_{1}^{2}}(y_{1}-x_{1}^{*})|<\epsilon^{\prime}$ for a sufficiently small $\epsilon^{\prime}>0$ , thereby

\Big{|}\cos(\zeta)\Re(e^{i\tau}f(y))-\Re(e^{i\tau}f(y))-\sin(\zeta)\Im(e^{i\tau}f(y))\Big{|}<\epsilon^{\prime\prime}|f|

for some sufficiently small $\epsilon^{\prime\prime}$ . Note that $f\in L^{\infty}(D)$ , therefore we can derive from (38) that

\displaystyle\left|\langle\Re(e^{i\tau}T_{\alpha})h,h\rangle_{L^{2}(D)}-\langle\Re(e^{i\tau}\widetilde{T})h,h\rangle_{L^{2}(D)}\right|\leq\epsilon^{\prime\prime\prime}\|h\|^{2}_{L^{2}(D)}

(39)

for some sufficiently small $\epsilon^{\prime\prime\prime}$ . Together with the assumption that $\Re(e^{i\tau}\widetilde{T})$ is coercive with estimate (34), we have that

\langle\Re(e^{i\tau}T_{\alpha})h,h\rangle_{L^{2}(D)}\geq c\|h\|^{2}_{L^{2}(D)},

for some positive constant $c=c^{\prime}-\epsilon^{\prime\prime\prime}>0$ . This completes the proof. ∎

Based on Theorem 5 and Lemma 6, a sufficient and less restrictive condition (that $\Re(e^{i\tau}f(x))>c_{3}>0$ almost everywhere for $x\in D$ with some $\tau$ ) on $f$ can be found to ensure that $\Re(e^{i\tau}T_{\alpha})$ is coercive.

6.3. Generalization of the main results

Based on the discussion of the multi-frequency operators $\widetilde{F}$ and $F_{\alpha}$ , we discuss below how to generalize the main Theorems 3 and 4. In particular, we can prove exactly in the same way that Theorem 3 holds when replacing $F^{1/2}$ by the self-adjoint operator $(F_{\alpha\#})^{1/2}$ or $(\widetilde{F}_{\#})^{1/2}$ , where $F_{\alpha\#}=\Re(e^{i\tau}F_{\alpha})$ and $\widetilde{F}_{\#}=\Re(e^{i\tau}\widetilde{F})$ . On the other hand, Theorem 4 for the factorization-based sampling method holds when replacing $F$ by $\Re(e^{i\tau}F_{\alpha})$ or $\Re(e^{i\tau}\widetilde{F})$ .

Acknowledgement

The author greatly thanks the anonymous referees for their helpful comments that have improved the quality of this manuscript.

7. Conclusion

The factorization method is investigated to treat the inverse source problem in waveguides with a single propagating mode at multiple frequencies. A factorization-based sampling method is also discussed. The main result of the factorization method is to provide both theoretical justifications and efficient imaging algorithms to image extended sources. One key step is to use the multi-frequency measurements to introduce the multi-frequency far-field operator, which is carefully designed to incorporate the possibly non-linear dispersion relation. We also emphasize that the measurements under consideration are “backscattering” measurements (in contrast to measurements obtained in two opposite directions), and we have shown how to use these less measurements to design appropriate multi-frequency operators to deploy the factorization method. Each multi-frequency operator can be factorized as $S^{*}TS$ which plays a fundamental role in the mathematically rigorous justification of the factorization method. A factorization-based sampling method is also discussed, where the factorization $S^{*}TS$ and a relevant point spread function play important roles. The sampling methods are capable to image the range support of the extended sources.

Another interesting observation is the potential to image a complete block, since the linear sampling method may not be capable of imaging such a complete block using far-field measurements at a single frequency. We provide examples to image both a complete source block and a complete obstacle block using the multi-frequency sampling methods. It is shown that the use of one propagating mode at multiple frequencies can efficiently locate a complete block where such efficiency may be of practical interests in long-distance communication optical devices or tunnels.

This work establishes the fundamental result with a single propagating mode at multiple frequencies. One direction of future work will be to investigate the factorization method in waveguide imaging using multiple propagating modes to establish theoretical justifications for extended scatterers in inverse obstacle/medium problems.

References

[1] A Alzaalig, G Hu, X Liu, and J Sun. Fast acoustic source imaging using multi-frequency sparse data, Inverse Problems 36, 025009, 2020.
[2] T Arens, D Gintides, and A Lechleiter. Direct and inverse medium scattering in a three-dimensional homogeneous planar waveguide. SIAM J. Appl. Math. 71(3), 753–772, 2011.
[3] AB Baggeroer, WA Kuperman, and PN Mikhalevsky. An overview of matched field methods in ocean acoustics. IEEE Journal of Oceanic Engineering 18(4), 401–424, 1993.
[4] C Bellis, M Bonnet, and F Cakoni. Acoustic inverse scattering using topological derivative of far-field measurements-based l2 cost functionals. Inverse Problems 29, 075012, 2013.
[5] L Borcea, F Cakoni, and S Meng. A direct approach to imaging in a waveguide with perturbed geometry, J. Comput. Phys. 392, 556–577, 2019.
[6] L Borcea and S Meng. Factorization method versus migration imaging in a waveguide. Inverse Problems 35(12),124006, 2019.
[7] L Borcea and DL Nguyen. Imaging with electromagnetic waves in terminating waveguides. Inverse Probl. Imaging 10, 915–941, 2016.
[8] L Bourgeois and S Fliss. On the identification of defects in a periodic waveguide from far field data. Inverse Problems 30(9), 095004, 2014.
[9] L Bourgeois, F Le Louër, and E Lunéville. On the use of Lamb modes in the linear sampling method for elastic waveguides. Inverse Problems 27(5), 055001, 2011.
[10] L Bourgeois and E Lunéville. The linear sampling method in a waveguide: a modal formulation. Inverse problems 24(1), 015018, 2008.
[11] L Bourgeois and E Lunéville. On the use of sampling methods to identify cracks in acoustic waveguides. Inverse Problems 28(10), 105011, 2012.
[12] L Bourgeois and E Lunéville. On the use of the linear sampling method to identify cracks in elastic waveguides. Inverse Problems 29(2), 025017, 2013.
[13] F Cakoni and D Colton. Qualitative Approach to Inverse Scattering Theory, Springer, 2016.
[14] F Cakoni, D Colton, and H Haddar. Inverse Scattering Theory and Transmission Eigenvalues, SIAM, 2016.
[15] F Cakoni, D Colton, and P Monk. The Linear Sampling Method in Inverse Electromagnetic Scattering, SIAM, 2011.
[16] J Chen and G Huang. A direct imaging method for inverse electromagnetic scattering problem in rectangular waveguide. Communications in Computational Physics 23, 1415–1433, 2017.
[17] L Chesnel, SA Nazarov and V Pagneux. Invisibility and perfect reflectivity in waveguides with finite length branches. SIAM Journal on Applied Mathematics 78(4), 2176–2199, 2018.
[18] D Colton and A Kirsch. A simple method for solving inverse scattering problems in the resonance region. Inverse Problems 12, 383–393, 1996.
[19] D Colton and R Kress. Inverse Acoustic and Electromagnetic Scattering Theory, Springer Nature, New York, 2019.
[20] ASBB Dhia, L Chesnel, and SA Nazarov. Perfect transmission invisibility for waveguides with sound hard walls. Journal de Mathématiques Pures et Appliquées 111(9), 79–105, 2018.
[21] R Griesmaier and C Schmiedecke. A factorization method for multi-frequency inverse source problems with sparse far-field measurements, SIAM J. Imag. Sci. 10, 2119–2139, 2017.
[22] A Haack, J Schreyer, and G Jackel. State-of-the-art of non-destructive testing methods for determining the state of a tunnel lining. Tunnelling and Underground Space Technology incorporating Trenchless Technology Research 4(10), 413–431, 1995.
[23] A Kirsch. Characterization of the shape of a scattering obstacle using the spectral data of the far-field operator, Inverse Problems 14, 1489–1512, 1998.
[24] A Kirsch and N Grinberg. The Factorization Method for Inverse Problems, Oxford University Press, Oxford, 2008.
[25] X Liu and S Meng. A multi-frequency sampling method for the inverse source problems with sparse measurements, arXiv:2109.01434.
[26] S Meng. A sampling type method in an electromagnetic waveguide. Inverse Probl. Imaging 15(4), 745–762, 2021.
[27] P Monk and V Selgas. Sampling type methods for an inverse waveguide problem. Inverse Probl. Imaging 6(4), 709–747, 2012.
[28] P Monk and V Selgas. An inverse acoustic waveguide problem in the time domain. Inverse Problems 32(5), 055001, 2016.
[29] P Monk, V Selgas, and F Yang. Near-field linear sampling method for an inverse problem in an electromagnetic waveguide. Inverse Problems 35(6), 065001, 2019.
[30] F Natterer. The Mathematics of Computerized Tomography, SIAM, 2001.
[31] P Rizzo, A Marzani, J Bruck, et al. Ultrasonic guided waves for nondestructive evaluation/structural health monitoring of trusses. Measurement science and technology 21(4), 045701, 2010.
[32] J Schöberl. Netgen an advancing front 2d/3d-mesh generator based on abstract rules. Computing and visualization in science 1(1), 41–52, 1997.
[33] J Sun and C Zheng. Reconstruction of obstacles embedded in waveguides. Contemporary Mathematics 586, 341–350, 2013.
[34] C Tsogka, DA Mitsoudis, and S Papadimitropoulos. Imaging extended reflectors in a terminating waveguide. SIAM J. Imag. Sci. 11(2), 1680–1716, 2018.
[35] Y Xu, C Mawata, and W Lin. Generalized dual space indicator method for underwater imaging. Inverse Problems 16(6),1761, 2000.

			$\displaystyle(Fg)(\sigma;x^{})=\int_{k_{-}}^{k_{+}}{-i}\mu_{1}(\omega_{\sigma\gamma})\Big{[}H(\sigma-\gamma)e^{i\theta}u_{p}^{s}\left(x^{};\omega_{\sigma\gamma}\right)-H(\gamma-\sigma)e^{-i\theta}\overline{u_{p}^{s}\left(x^{*};\omega_{\sigma\gamma}\right)}\Big{]}g(\gamma)\,\mbox{d}\gamma$
		$\displaystyle=$	$\displaystyle 1/2\int_{k_{-}}^{k_{+}}\int_{D}\psi_{1}(y_{\perp})e^{i(\sigma-\gamma)\|x^{}_{1}-y_{1}\|}e^{i\theta}{\psi_{1}(x^{}_{\perp})}f(y)g(\gamma)\,\mbox{d}y\,\mbox{d}\gamma$
		$\displaystyle=$	$\displaystyle\int_{D}e^{i\sigma\|x^{}_{1}-y_{1}\|}\overline{\sqrt{\psi_{1}(y_{\perp})}}\left[\left(\int_{k_{-}}^{k_{+}}\sqrt{\psi_{1}(y_{\perp})}e^{-i\gamma\|x^{}_{1}-y_{1}\|}g(\gamma)\,\mbox{d}\gamma\right)\frac{e^{i\theta}{\psi_{1}(x^{*}_{\perp})}f(y)}{2}\right]\,\mbox{d}y,$