An Efficient Method for Local Base Transform in Pekeris Waveguide with Radiation Condition

Jianxin Zhu^*, Haibo Chen

Department of Mathematics, Jinan University, Guangzhou, China.

DOI: 10.4236/jamp.2020.85060                

Abstract

In this paper, for the fast computation of the coordinates under the basis of the eigenfunctions of Helmholtz operator, we derive the conjugate operator with the radiation boundary condition. Further, we prove the cross orthogonality between the linearly-independent eigenfunctions of the Helmholtz operator and the linearly-independent eigenfunctions of the conjugate operator. The numerical simulations demonstrate the effectiveness of our treatment.

Keywords

Helmholtz Operator, Radiation Boundary Condition, Conjugate Operator, Eigenfunction Expansion, Orthogonality, Pekeris Waveguide

1. Introduction

It is known that the computation of the wave propagation plays a particularly prominent role in ocean acoustics or optical waveguides [1] [2] [3]. Because of the inhomogeneity of the ocean medium, the sound velocity forms the underwater channel with the regular change of the ocean, and the remote acoustic propagation of hundreds of or even thousands of kilometers can be carried out by using the underwater channel. The study of acoustic wave propagation in the ocean is of great significance to such things as underwater positioning, communication and ranging [4].

The propagation of sound waves in the medium is mathematically formulated as a well-known Helmholtz equation [5]. When the propagation area or the solution region is bounded, the linearly-independent eigenfunctions of the Helmholtz operator are weighted orthogonal, so one can use the marching algorithm to solve the Helmholtz equation quickly [1], and a good numerical solution has been obtained. However, in reality, it is extremely difficult to deal with such a problem when considering the propagation region is unbounded in the transverse direction, such as the infinite depth of the ocean (including the submarine rock layer) or open light region in the upper or lower direction [6]. For the convenience of calculation, we need to truncate the unbounded waveguide area into a bounded waveguide area. It is common practice to add some artificial boundary conditions to the outside of the truncation area [7]. In order for the acoustic propagation of the truncated area to be consistent with that before truncation, special processing needs to be added to the artificial boundary so that the sound waves do not have obvious reflection when they leave the area, and the sound waves can be absorbed effectively.

For the addition of artificial boundary conditions, we use the radiation boundary condition (RBC) [8] to truncate the unbounded waveguide. However, the orthogonality of between the linearly-independent eigenfunctions of the Helmholtz operator is lost, so it is difficult to fast compute the coordinate coefficients under the local eigenfunctions’ bases when the Helmholtz equation with the RBC is solved by some numerical marching method, such as the Operator Marching Method and the One-way Method [9].

For this reason, the conjugate operator of the Helmholtz operator with the RBC is constructed in this paper. Further, it is proved that there is the crossly-orthogonal property between the linearly-independent eigenfunctions of the Helmholtz operator and their linearly-independent conjugate eigenfunctions. Furthermore, a simple and convenient formula is given to the coordinate transformation under two different eigenfunction bases. Finally, the numerical experiment results verify the correctness of this treatment.

2. Mathematical Treatment

2.1. Eigenfunctions of the Helmholtz Operator

For simplicity, we consider the wave propagation in the Pekeris waveguide. The mathematical model is as follows:

$(\begin{array}{l}\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial z^2}+{\kappa }_1^{2}u=0,0<z<G;\\ \frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial z^2}+{\kappa }_2^{2}u=0,z>G;\\ u\left(x,0\right)=0,\text{the}\text{boundary}\text{condition}\text{at}z=0;\\ \underset{z\to G^{-}}{\lim }u=\underset{z\to G^{+}}{\lim }u,\underset{z\to G^{-}}{\lim }\frac{1}{{\rho }_1}\cdot \frac{\partial u}{\partial z}=\underset{z\to G^{+}}{\lim }\frac{1}{{\rho }_2}\cdot \frac{\partial u}{\partial z},\text{the}\text{interface}\text{conditions}\text{at}z=G;\\ \underset{z\to +\infty }{\lim }u\left(x,z\right)=0.\end{array}$ (1)

Where z is the depth axis pointing downward with the ocean surface at $z=0$, and $x$ is the range variable in horizontal direction. One layer with density ${\rho }_1$ is located in $0<z<G$, the other with density ${\rho }_2$ is located in $z>G$ ; and their interface is flat at $z=G$. Here, ${\kappa }_1$ and ${\kappa }_2$ are represented as the wavenumbers in two different subregions, respectively.

We introduce the RBC for truncation at $z>G$, that is $\frac{\partial u}{\partial z}=i{\gamma }_2\cdot u$ at $z\ge H$, where $H>G$, $i=\sqrt{-1}$, ${\gamma }_2=\sqrt{{\kappa }_2^2-\lambda }$, and $\lambda$ is the eigenvalue of the Helmholtz operator $P$. Here $P=\frac{{\text{d}}^2}{\text{d}{z}^2}+{\kappa }^2\left(z\right)$, where $\kappa \left(z\right)$ is ${\kappa }_1$ as $0<z<G$ and is ${\kappa }_2$ as $z>G$. Then Equation (1) is approximately transformed into a complex partial differential equation in the bounded region as follows:

$(\begin{array}{l}\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial z^2}+{\kappa }_1^{2}u=0,0<z<G;\\ \frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial z^2}+{\kappa }_2^{2}u=0,G<z<H;\\ u\left(x,0\right)=0,\text{the}\text{boundary}\text{condition}\text{at}z=0;\\ \underset{z\to G^{-}}{\lim }u=\underset{z\to G^{+}}{\lim }u,\underset{z\to G^{-}}{\lim }\frac{1}{{\rho }_1}\cdot \frac{\partial u}{\partial z}=\underset{z\to G^{+}}{\lim }\frac{1}{{\rho }_2}\cdot \frac{\partial u}{\partial z},\text{the}\text{interface}\text{conditions}\text{at}z=G;\\ \frac{\partial u}{\partial z}=i{\gamma }_2\cdot u,z\ge H,\text{the}\text{boundary}\text{condition}\text{at}z\ge H.\end{array}$ (2)

For the coordinate transformation between two different bases, it is necessary to quickly calculate the coordinate coefficients of any given function under a basis. Here the basis is formed by the linearly-independent eigenfunctions of the Helmholtz operator $P$. If we assume the solution of Equation (2) to be $u=\varphi \left(z\right)\cdot exp\left(i\sqrt{\lambda }\cdot x\right)$, then the eigenvalue problem of the Equation (2) is derived as follows:

$(\begin{array}{l}\frac{{\text{d}}^2\varphi }{\text{d}{z}^2}+{\kappa }_1^2\varphi =\lambda \varphi ,0<z<G;\\ \frac{{\text{d}}^2\varphi }{\text{d}{z}^2}+{\kappa }_2^2\phi =\lambda \varphi ,G<z<H;\\ \varphi \left(0\right)=0,\text{the}\text{boundary}\text{condition}\text{at}z=0;\\ \underset{z\to G^{-}}{\lim }\varphi =\underset{z\to G^{+}}{\lim }\varphi ,\underset{z\to G^{-}}{\lim }\frac{1}{{\rho }_1}\cdot \frac{\text{d}\varphi }{\text{d}z}=\underset{z\to G^{+}}{\lim }\frac{1}{{\rho }_2}\cdot \frac{\text{d}\varphi }{\text{d}z},\text{the}\text{interface}\text{conditions}\text{at}z=G;\\ \frac{\text{d}\varphi }{\text{d}z}=i{\gamma }_2\varphi ,z\ge H,\text{the}\text{boundary}\text{condition}\text{at}z\ge H.\end{array}$ (3)

Where $\varphi \left(z\right)$ is the eigenfunction corresponding to the eigenvalue $\lambda$ for the eigenvalue problem of Equation (2).

Given $u\left(x,z\right)\approx {\displaystyle {\sum }_{j=1}^m}{c}_j\left(x\right){\varphi }_j\left(z\right)$ [10], if the solution functions ${\varphi }_j\left(z\right)$, $\left(j=\mathrm{1,2,}\cdots \mathrm{,}m\right)$ of Equation (3) are crossly orthogonal to their conjugate functions ${\phi }_i\left(z\right)$, $\left(i=\mathrm{1,2,}\cdots \mathrm{,}m\right)$ with the weight $\frac{1}{\rho \left(z\right)}$, that is, ${\displaystyle {\int }_0^H}\frac{1}{\rho \left(z\right)}\cdot {\varphi }_j\left(z\right)\overline{{\phi }_i\left(z\right)}\text{d}z=0$ as $i\ne j$, then we have $c_j\left(x\right)\approx \frac{{\displaystyle {\int }_0^H}\frac{1}{\rho \left(z\right)}\cdot u\left(x\mathrm{,}z\right)\overline{{\phi }_j\left(z\right)}\text{d}z}{{\displaystyle {\int }_0^H}\frac{1}{\rho \left(z\right)}{\varphi }_j\left(z\right)\overline{{\phi }_j\left(z\right)}\text{d}z}$, where each eigenfunction ${\varphi }_j\left(z\right)$ corresponds to each eigenvalue ${\lambda }_j$, respectively, $\left(j=\mathrm{1,2,}\cdots \mathrm{,}m\right)$ ; and the overline represents the conjugate operation.

The following processes focus on constructing the conjugate eigenfunctions ${\phi }_i\left(z\right)$, $\left(i=\mathrm{1,2,}\cdots \mathrm{,}m\right)$ corresponding to the eigenfunctions ${\varphi }_j\left(z\right)$, $\left(j=\mathrm{1,2,}\cdots \mathrm{,}m\right)$, and further prove their cross orthogonality.

For Equation (3), we have the general solutions as follows:

$\varphi \left(z\right)=(\begin{array}{l}{F}_1\cdot \sin \left(\sqrt{{\kappa }_1^2-\lambda }\cdot z\right),0\le z<G;\\ F_2\cdot \exp \left(i{\gamma }_2\cdot z\right)+F_3\cdot \exp \left(-i{\gamma }_2\cdot z\right),G<z<H;\end{array}$ (4)

where

$(\begin{array}{l}{F}_2=F_1/2\cdot exp\left(-i{\gamma }_2\cdot G\right)\cdot \left[sin\left({\gamma }_1\cdot G\right)-i\cdot \frac{{\rho }_2}{{\rho }_1}\cdot \frac{{\gamma }_1}{{\gamma }_2}\cdot cos\left({\gamma }_1\cdot G\right)\right]\mathrm{,}\\ F_2=F_1/2\cdot exp\left(i{\gamma }_2\cdot G\right)\cdot \left[sin\left({\gamma }_1\cdot G\right)+i\cdot \frac{{\rho }_2}{{\rho }_1}\cdot \frac{{\gamma }_1}{{\gamma }_2}\cdot cos\left({\gamma }_1\cdot G\right)\right]\mathrm{,}\end{array}$ (5)

and ${\gamma }_1=\sqrt{{\kappa }_1^2-\lambda }$.

By the RBC at $z=H$, namely ${\frac{\text{d}\varphi }{\text{d}z}|}_{z=H}=i{\gamma }_2\cdot \varphi \left(H\right)$, we can obtain the following equation of the eigenvalue $\lambda$ :

$\text{tan}\left({\gamma }_1\cdot G\right)=-i\cdot \frac{{\rho }_2}{{\rho }_1}\cdot \frac{{\gamma }_1}{{\gamma }_2}\mathrm{.}$ (6)

2.2. The Construction of the Conjugate Operator

Let the operator $P$ be with the RBC at $z\mathrm{<mo>\; =\; </mo>}H$ and the interface conditions at $z\mathrm{<mo>\; =\; </mo>}G$, we have

$P\varphi =\frac{\text{d}}{\text{d}z}\left[\frac{\text{d}\varphi }{\text{d}z}\right]+{\hat{\kappa }}^2\left(z\right)\varphi \mathrm{,}$ (7)

where

$\hat{\kappa }\left(z\right)=(\begin{array}{cc}{\kappa }_1,& 0<z<G;\\ {\kappa }_2,& G<z<H.\end{array}$ (8)

Let $\varphi$ and $\phi$ be the eigenfunction and the corresponding conjugate eigenfunction of the operator $P$, respectively, where the boundary conditions are

given as $\varphi \left(0\right)=\mathrm{0,}\frac{\text{d}\varphi }{\text{d}z}\left(H\right)=i{\gamma }_2\cdot \varphi \left(H\right)$.

Denote

$\rho \left(z\right)=(\begin{array}{cc}{\rho }_1,& 0<z<G;\\ {\rho }_2,& G<z\le H.\end{array}$ (9)

Consider

$\begin{array}{c}\langle P\varphi \mathrm{,}\phi \rangle ={\displaystyle {\int }_0^H}\frac{1}{\rho \left(z\right)}\bar{\phi }\cdot P\varphi \text{d}z={\displaystyle {\int }_0^H}\frac{1}{\rho \left(z\right)}\cdot \left[\frac{{\text{d}}^2\varphi }{\text{d}{z}^2}+{\hat{\kappa }}^2\left(z\right)\varphi \right]\bar{\phi }\text{d}z\\ ={\displaystyle {\int }_0^H}\frac{1}{\rho }\cdot \frac{{\text{d}}^2\varphi }{\text{d}{z}^2}\cdot \bar{\phi }\text{d}z+{\displaystyle {\int }_0^H}\frac{1}{\rho }\cdot {\hat{\kappa }}^2\left(z\right)\varphi \cdot \bar{\phi }\text{d}z\\ ={\displaystyle {\int }_0^G}\frac{1}{{\rho }_1}\cdot \frac{{\text{d}}^2\varphi }{\text{d}{z}^2}\cdot \bar{\phi }\text{d}z+{\displaystyle {\int }_G^H}\frac{1}{{\rho }_2}\cdot \frac{{\text{d}}^2\varphi }{\text{d}{z}^2}\cdot \bar{\phi }\text{d}z+{\displaystyle {\int }_0^H}\frac{1}{\rho }\cdot {\hat{\kappa }}^2\left(z\right)\varphi \cdot \bar{\phi }\text{d}z\\ \triangleq S_1+S_2+S_3\mathrm{,}\end{array}$ (10)

where

$S_1={\displaystyle {\int }_0^G}\frac{1}{{\rho }_1}\cdot \frac{{\text{d}}^2\varphi }{\text{d}{z}^2}\cdot \bar{\phi }\text{d}z\mathrm{,}{S}_2={\displaystyle {\int }_G^H}\frac{1}{{\rho }_2}\cdot \frac{{\text{d}}^2\varphi }{\text{d}{z}^2}\cdot \bar{\phi }\text{d}z\mathrm{,}{S}_3={\displaystyle {\int }_0^H}\frac{1}{\rho }{\hat{\kappa }}^2\left(z\right)\cdot \varphi \cdot \bar{\phi }\text{d}z\mathrm{.}$

Then, by the integration by parts and making use of the interface and the boundary conditions, we have

$\begin{array}{c}{S}_1={\displaystyle {\int }_0^G}\frac{1}{{\rho }_1}\cdot \frac{{\text{d}}^2\varphi }{\text{d}{z}^2}\cdot \bar{\phi }\text{d}z\\ ={\left(\frac{1}{{\rho }_1}\bar{\phi }\frac{\text{d}\varphi }{\text{d}z}\right)|}_0^{G^{-}}-{\displaystyle {\int }_0^G}\frac{1}{{\rho }_1}\frac{\text{d}\bar{\phi }}{\text{d}z}\cdot \left[\frac{\text{d}\varphi }{\text{d}z}\right]\text{d}z\\ ={\left(\frac{1}{{\rho }_1}\bar{\phi }\frac{\text{d}\varphi }{\text{d}z}\right)|}_{z=G^{-}}-{\displaystyle {\int }_0^G}\frac{1}{{\rho }_1}\frac{\text{d}\bar{\phi }}{\text{d}z}\cdot \left[\frac{\text{d}\varphi }{\text{d}z}\right]\text{d}z\\ ={\left(\frac{1}{{\rho }_1}\bar{\phi }\frac{\text{d}\varphi }{\text{d}z}\right)|}_{z=G^{-}}-{\left(\frac{1}{{\rho }_1}\frac{\text{d}\bar{\phi }}{\text{d}z}\cdot \varphi \right)|}_0^{G^{-}}+{\displaystyle {\int }_0^G}\varphi \cdot \frac{1}{{\rho }_1}\frac{{\text{d}}^2\bar{\phi }}{\text{d}{z}^2}\text{d}z\\ ={\left(\frac{1}{{\rho }_1}\bar{\phi }\frac{\text{d}\varphi }{\text{d}z}\right)|}_{z=G^{-}}-{\left(\frac{1}{{\rho }_1}\frac{\text{d}\bar{\phi }}{\text{d}z}\cdot \varphi \right)|}_{z=G^{-}}+{\displaystyle {\int }_0^G}\varphi \cdot \frac{1}{{\rho }_1}\frac{{\text{d}}^2\bar{\phi }}{\text{d}{z}^2}\text{d}z\mathrm{,}\end{array}$ (11)

where $\bar{\phi }\left(0\right)=0$ is required.

Also, we have

$\begin{array}{c}{S}_2={\displaystyle {\int }_G^H}\frac{1}{{\rho }_2}\cdot \left[\frac{\text{d}}{\text{d}z}\left[\frac{\text{d}\varphi }{\text{d}z}\right]\right]\cdot \bar{\phi }\text{d}z\\ ={\left(\frac{1}{{\rho }_2}\bar{\phi }\frac{\text{d}\varphi }{\text{d}z}\right)|}_{G^{+}}^H-{\displaystyle {\int }_G^H}\frac{1}{{\rho }_2}\frac{\text{d}\varphi }{\text{d}z}\frac{\text{d}\bar{\phi }}{\text{d}z}\text{d}z\\ ={\left(\frac{1}{{\rho }_2}\bar{\phi }\frac{\text{d}\varphi }{\text{d}z}\right)|}_{G^{+}}^H-{\left(\frac{1}{{\rho }_2}\frac{\text{d}\bar{\phi }}{\text{d}z}\cdot \varphi \right)|}_{G^{+}}^H+{\displaystyle {\int }_G^H}\varphi \frac{1}{{\rho }_2}\frac{{\text{d}}^2\bar{\phi }}{\text{d}{z}^2}\text{d}z\mathrm{.}\end{array}$ (12)

Thus, we have

$\begin{array}{c}\langle P\varphi \mathrm{,}\phi \rangle =S_1+S_2+S_3\\ ={\left(\frac{1}{{\rho }_1}\bar{\phi }\frac{\text{d}\varphi }{\text{d}z}\right)|}_{z=G^{-}}-{\left(\frac{1}{{\rho }_1}\frac{\text{d}\bar{\phi }}{\text{d}z}\cdot \varphi \right)|}_{z=G^{-}}+{\displaystyle {\int }_0^G}\varphi \cdot \frac{1}{{\rho }_1}\frac{{\text{d}}^2\bar{\phi }}{\text{d}{z}^2}\text{d}z\\ +{\left(\frac{1}{{\rho }_2}\bar{\phi }\frac{\text{d}\varphi }{\text{d}z}\right)|}_{G^{+}}^H-{\left(\frac{1}{{\rho }_2}\frac{\text{d}\bar{\phi }}{\text{d}z}\cdot \varphi \right)|}_{G^{+}}^H+{\displaystyle {\int }_G^H}\varphi \frac{1}{{\rho }_2}\frac{{\text{d}}^2\bar{\phi }}{\text{d}{z}^2}\text{d}z\\ +{\displaystyle {\int }_0^H}\frac{1}{\rho }{\hat{\kappa }}^2\left(z\right)\varphi \cdot \bar{\phi }\text{d}z\mathrm{.}\end{array}$ (13)

Applying the interface and the boundary conditions:

$\varphi \left(G^{-}\right)=\varphi \left(G^{+}\right)\mathrm{,}\frac{1}{{\rho }_1}{\varphi }_z\left(G^{-}\right)=\frac{1}{{\rho }_2}{\varphi }_z\left(G^{+}\right)\text{and}{\varphi }_z\left(H\right)=i{\gamma }_2\varphi \left(H\right)\mathrm{,}$ (14)

the above expression is simplified as follows:

$\begin{array}{c}\langle P\varphi \mathrm{,}\phi \rangle ={\left(\frac{1}{{\rho }_1}\bar{\phi }\frac{\text{d}\varphi }{\text{d}z}\right)|}_{z=G^{-}}+{\left(\frac{1}{{\rho }_2}\bar{\phi }\frac{\text{d}\varphi }{\text{d}z}\right)|}_{z=H}-{\left(\frac{1}{{\rho }_2}\bar{\phi }\frac{\text{d}\varphi }{\text{d}z}\right)|}_{z=G^{+}}\\ -{\left(\frac{1}{{\rho }_1}\frac{\text{d}\bar{\phi }}{\text{d}z}\cdot \varphi \right)|}_{z=G^{-}}-{\left(\frac{1}{{\rho }_2}\frac{\text{d}\bar{\phi }}{\text{d}z}\cdot \varphi \right)|}_{z=H}+{\left(\frac{1}{{\rho }_2}\frac{\text{d}\bar{\phi }}{\text{d}z}\cdot \varphi \right)|}_{z=G^{+}}\\ +{\displaystyle {\int }_0^G}\varphi \cdot \frac{1}{{\rho }_1}\frac{{\text{d}}^2\bar{\phi }}{\text{d}{z}^2}\text{d}z+{\displaystyle {\int }_G^H}\varphi \frac{1}{{\rho }_2}\frac{{\text{d}}^2\bar{\phi }}{\text{d}{z}^2}\text{d}z+{\displaystyle {\int }_0^H}\frac{1}{\rho }{\hat{\kappa }}^2\left(z\right)\varphi \cdot \bar{\phi }\text{d}z\end{array}$

$\begin{array}{l}=\left[{\left(\frac{1}{{\rho }_1}\bar{\phi }\frac{\text{d}\varphi }{\text{d}z}\right)|}_{z=G^{-}}-{\left(\frac{1}{{\rho }_2}\bar{\phi }\frac{\text{d}\varphi }{\text{d}z}\right)|}_{z=G^{+}}\right]\\ +\left[{\left(\frac{1}{{\rho }_2}\frac{\text{d}\bar{\phi }}{\text{d}z}\cdot \varphi \right)|}_{z=G^{+}}-{\left(\frac{1}{{\rho }_1}\frac{\text{d}\bar{\phi }}{\text{d}z}\cdot \varphi \right)|}_{z=G^{-}}\right]\\ +\left[{\left(\frac{1}{{\rho }_2}\bar{\phi }\frac{\text{d}\varphi }{\text{d}z}\right)|}_{z=H}-{\left(\frac{1}{{\rho }_2}\frac{\text{d}\bar{\phi }}{\text{d}z}\cdot \varphi \right)|}_{z=H}\right]\\ +{\displaystyle {\int }_0^G}\varphi \cdot \frac{1}{{\rho }_1}\frac{{\text{d}}^2\bar{\phi }}{\text{d}{z}^2}\text{d}z+{\displaystyle {\int }_G^H}\varphi \frac{1}{{\rho }_2}\frac{{\text{d}}^2\bar{\phi }}{\text{d}{z}^2}\text{d}z+{\displaystyle {\int }_0^H}\frac{1}{\rho }{\hat{\kappa }}^2\left(z\right)\varphi \bar{\phi }\text{d}z\mathrm{.}\end{array}$ (15)

For simplicity, the interface and the boundary conditions are given to the conjugate eigenfunction $\phi$ as follows:

let

${\frac{1}{{\rho }_1}\bar{\phi }\frac{\text{d}\varphi }{\text{d}z}|}_{z=G^{-}}-{\left(\frac{1}{{\rho }_2}\bar{\phi }\frac{\text{d}\varphi }{\text{d}z}\right)|}_{z=G^{+}}=\mathrm{0,}$ (16)

it is equivalent to

$\bar{\phi }\left(G^{-}\right)=\bar{\phi }\left(G^{+}\right)\text{or}\phi \left(G^{-}\right)=\phi \left(G^{+}\right)\mathrm{<mo>\; ;\; </mo>}$ (17)

let

${\frac{1}{{\rho }_2}\frac{\text{d}\bar{\phi }}{\text{d}z}\cdot \varphi |}_{z=G^{+}}-{\left(\frac{1}{{\rho }_1}\frac{\text{d}\bar{\phi }}{\text{d}z}\cdot \varphi \right)|}_{z=G^{-}}=\mathrm{0,}$ (18)

it is equivalent to

$\frac{1}{{\rho }_1}{\bar{\phi }}_z\left(G^{-}\right)=\frac{1}{{\rho }_2}{\bar{\phi }}_z\left(G^{+}\right)\text{or}\frac{1}{{\rho }_1}{\phi }_z\left(G^{-}\right)=\frac{1}{{\rho }_2}{\phi }_z\left(G^{+}\right)\mathrm{<mo>\; ;\; </mo>}$ (19)

and let

${\frac{1}{{\rho }_2}\bar{\phi }\frac{\text{d}\varphi }{\text{d}z}|}_{z=H}-{\left(\frac{1}{{\rho }_2}\frac{\text{d}\bar{\phi }}{\text{d}z}\cdot \varphi \right)|}_{z=H}=\mathrm{0,}$ (20)

it is equivalent to

${\bar{\phi }}_z\left(H\right)=i{\gamma }_2\bar{\phi }\left(H\right)\text{or}{\phi }_z\left(H\right)=-i\overline{{\gamma }_2}\phi \left(H\right)\mathrm{.}$ (21)

Thus, by the interface and the boundary conditions of the function $\phi$ as above, then we have

$\begin{array}{c}\langle P\varphi \mathrm{,}\phi \rangle ={\displaystyle {\int }_0^H}\frac{1}{\rho }\varphi \left[\frac{\text{d}}{\text{d}z}\left(\frac{\text{d}\bar{\phi }}{\text{d}z}\right)+{\hat{\kappa }}^2\left(z\right)\bar{\phi }\right]\text{d}z\\ \triangleq {\displaystyle {\int }_0^H}\frac{1}{\rho }\varphi \cdot \overline{Q\phi }\text{d}z\\ \mathrm{<mo>\; =\; </mo>}\langle \varphi \mathrm{,}Q\phi \rangle \mathrm{,}\end{array}$ (22)

where

$\overline{Q\phi }=\frac{{\text{d}}^2\bar{\phi }}{\text{d}{z}^2}+{\hat{\kappa }}^2\left(z\right)\bar{\phi }\mathrm{.}$ (23)

Finally, we have

$\langle P\varphi \mathrm{,}\phi \rangle =\langle \varphi \mathrm{,}Q\phi \rangle \mathrm{.}$ (24)

In fact, there is the following relation between the operators P and Q:

$\bar{Q}=\text{P}=\frac{{\text{d}}^2}{\text{d}{z}^2}+{\hat{\kappa }}^2\left(z\right)\mathrm{.}$ (25)

Namely, we have

${\displaystyle {\int }_0^H}\frac{1}{\rho \left(z\right)}\cdot \bar{\phi }\cdot P\varphi \text{d}z={\displaystyle {\int }_0^H}\frac{1}{\rho \left(z\right)}\cdot \overline{Q\phi }\cdot \varphi \text{d}z\mathrm{.}$ (26)

Obviously, $\phi$ is the conjugate function of $\varphi$, that is also called the conjugate eigenfunction of the eigenfunction $\varphi$. In this paper, we assume that $\rho \left(z\right)$ is a piecewise constant function. Under above assumptions, $\bar{\phi }$ should satisfy the following ordinary differential equation:

$(\begin{array}{l}{\bar{\phi }}_{zz}+{\kappa }_1^2\bar{\phi }=\lambda \bar{\phi },0<z<G;\\ {\bar{\phi }}_{zz}+{\kappa }_2^2\bar{\phi }=\lambda \bar{\phi },G<z<H;\\ \bar{\phi }\left(0\right)=0;\\ \underset{z\to G^{-}}{lim}\bar{\phi }=\underset{z\to G^{+}}{lim}\bar{\phi }\mathrm{,}\underset{z\to G^{-}}{lim}\frac{1}{{\rho }_1}\frac{\text{d}\bar{\phi }}{\text{d}z}=\underset{z\to G^{+}}{lim}\frac{1}{{\rho }_2}\frac{\text{d}\bar{\phi }}{\text{d}z}\mathrm{<mo>\; ;\; </mo>}\\ {\bar{\phi }}_z=i{\gamma }_2\bar{\phi }\mathrm{,}z\ge H\mathrm{.}\end{array}$ (27)

We have the similar solution form of Equation (27) as Equation (3):

$\bar{\phi }\left(z\right)=(\begin{array}{l}C\sin \left({\gamma }_1\cdot z\right),0\le z<G,\\ C_1{\text{e}}^{i{\gamma }_2\cdot z}+C_2{\text{e}}^{-i{\gamma }_2\cdot z}\mathrm{,}G<z\le H\mathrm{,}\end{array}$ (28)

where the parameters $C\mathrm{,}{C}_1$ and $C_2$ are to be determined.

In general, we choose $C=1$, then we can write down the solution as

$\bar{\phi }\left(z\right)=(\begin{array}{l}\sin \left({\gamma }_1\cdot z\right),0\le z<G,\\ C_1{\text{e}}^{i{\gamma }_2\cdot z}+C_2{\text{e}}^{-i{\gamma }_2\cdot z},G<z\le H.\end{array}$ (29)

Making ${\bar{\phi }}_z\left(H\right)=i{\gamma }_2\bar{\phi }\left(H\right)$, we can get the same nonlinear equation of the eigenvalue $\lambda$ as Equation (6). As a result, we give Theorem 1 as follows.

Theorem 1: For the above-mentioned Helmholtz operator P, then its conjugate operator is Q.

Remark 1: Although the two operators P and Q have the same representation, their RBCs are different.

2.3. Proof of the Cross Orthogonality between the Eigenfunctions and Their Conjugate Eigenfunctions

In this section, we prove the cross orthogonality of the linearly-independent eigenfunctions ${\varphi }_j\left(z\right)$ and their linearly-independent conjugate eigenfunctions ${\phi }_i\left(z\right)$, $\left(i,j=1,2,\cdots ,m;i\ne j\right)$.

Consider the eigenvalues and eigenfunctions of both ${\varphi }_j\left(z\right)$ and ${\phi }_i\left(z\right)$, which satisfy

$P{\varphi }_j={\lambda }_j{\varphi }_j\mathrm{,}\overline{Q{\phi }_i}={\lambda }_i\overline{{\phi }_i}\mathrm{,}\text{for}{\lambda }_i\ne {\lambda }_j\mathrm{,}i\ne j\mathrm{.}$

This give rises to

${\displaystyle {\int }_0^H}\frac{1}{\rho }\overline{{\phi }_i}\cdot {\varphi }_j\text{d}z=0.$ (30)

This is because

$\langle P{\varphi }_j\mathrm{,}{\phi }_i\rangle ={\displaystyle {\int }_0^H}\frac{1}{\rho }\cdot \overline{{\phi }_i}\cdot \text{P}{\varphi }_j\text{d}z={\displaystyle {\int }_0^H}\frac{1}{\rho }\cdot \overline{{\phi }_i}\cdot {\lambda }_j{\varphi }_j\text{d}z={\lambda }_j{\displaystyle {\int }_0^H}\frac{1}{\rho }\cdot \overline{{\phi }_i}\cdot {\varphi }_j\text{d}z\mathrm{,}$ (31)

and

$\langle {\varphi }_j\mathrm{,}Q{\phi }_i\rangle ={\displaystyle {\int }_0^H}\frac{1}{\rho }\cdot \overline{\text{Q}{\phi }_i}\cdot {\varphi }_j\text{d}z={\displaystyle {\int }_0^H}\frac{1}{\rho }\cdot {\lambda }_i\overline{{\phi }_i}\cdot {\varphi }_j\text{d}z={\lambda }_i{\displaystyle {\int }_0^H}\frac{1}{\rho }\cdot \overline{{\phi }_i}\cdot {\varphi }_j\text{d}z\mathrm{.}$ (32)

By

$\langle P{\varphi }_j\mathrm{,}{\phi }_i\rangle =\langle {\varphi }_j\mathrm{,}Q{\phi }_i\rangle \mathrm{,}$ (33)

we have

${\lambda }_j{\displaystyle {\int }_0^H}\frac{1}{\rho }\cdot \overline{{\phi }_i}\cdot {\varphi }_j\text{d}z={\lambda }_i{\displaystyle {\int }_0^H}\frac{1}{\rho }\cdot \overline{{\phi }_i}\cdot {\varphi }_j\text{d}z\mathrm{,}$ (34)

that is,

$\left({\lambda }_i-{\lambda }_j\right)\cdot {\displaystyle {\int }_0^H}\frac{1}{\rho }\cdot \overline{{\phi }_i}\cdot {\varphi }_j\text{d}z=0.$ (35)

Thus

${\displaystyle {\int }_0^H}\frac{1}{\rho }\cdot \overline{{\phi }_i}\cdot {\varphi }_j\text{d}z=0.$ (36)

At last, we also have the following conclusion.

Theorem 2: For the linearly-independent eigenfunctions ${\varphi }_j\left(z\right)$ and their linearly-independent conjugate eigenfunctions ${\phi }_i\left(z\right)$, $\left(i\mathrm{,}j=\mathrm{1,2,}\cdots \mathrm{,}m\mathrm{<mo>\; ;\; </mo>}i\ne j\right)$, as mentioned above, then there is the cross orthogonality, that is, ${\displaystyle {\int }_0^H}\frac{1}{\rho }\cdot \overline{{\phi }_i}\cdot {\varphi }_j\text{d}z=0$, $\left(i\mathrm{,}j=\mathrm{1,2,}\cdots \mathrm{,}m\mathrm{<mo>\; ;\; </mo>}i\ne j\right)$.

3. Numerical Examples

In this section, we will verify the method presented in the previous sections. Considering a Pekeris waveguide, we let ${\kappa }_1=16,{\kappa }_2=0.7\times 16,{\rho }_1=1,{\rho }_2=1.7,G=1,H=2.5$. These parameters are taken from reference [1]. By the nonlinear Equation (6), we define

$f\left(\lambda \right)=\tan \left({\gamma }_1\cdot G\right)+i\cdot \frac{{\rho }_2}{{\rho }_1}\cdot \frac{{\gamma }_1}{{\gamma }_2}\text{and}\text{AE}=\left|f\left(\lambda \right)\right|.$

We find out the roots of the equation $f\left(\lambda \right)=0$ by Newton’s iteration method and perform ten iterations for each root, where the initial values for the leaky modes are used by the asymptotic solutions, and initial values for the propagating modes are chosen by some appropriate values within the interval $\left[min\left({\kappa }_1^2\mathrm{,}{\kappa }_2^2\right)\mathrm{,}max\left({\kappa }_1^2\mathrm{,}{\kappa }_2^2\right)\right]$. The details of choosing initial values are listed in the reference [3]. As a result, we obtain the eigenvalues’ distributions of the operators $P$ and $\bar{Q}$, which are shown in Figure 1. And these distributions really fit the physical meaning.

To numerically verify the cross orthogonality, we list some eigenvalues of the operator P in detail as follows:

${\lambda }_1=\mathrm{248.4729239386205,}$

${\lambda }_2=\mathrm{225.4329511872122,}$

${\lambda }_3=-298.5970311054733+26.80248573108962i\mathrm{,}$

${\lambda }_4=-160.3961093395291+21.76064335348847i\mathrm{,}$

${\lambda }_5=\mathrm{132.0985417466585,}$

${\lambda }_6=\mathrm{186.1806848042874,}$

Figure 1. The original eigenvalues of the operator $P$ represented by “o” and the eigenvalues of the operator $\bar{Q}$ represented by “+”.

${\lambda }_7=-41.927161372191101+16.370663337533994i\mathrm{,}$

${\lambda }_8=56.807916692538242+10.204330264562326i\mathrm{,}$

where their error estimations about AE are listed in Table 1.

Then we choose the adaptive Gauss integral (AGI) method and numerical computation to verify the orthogonality of the eigenfunctions’ set

$\left\{{\varphi }_j\left(z\right)\mathrm{<mo>\; |\; </mo>}j=\mathrm{1,2,}\cdots \mathrm{,8}\right\}$ and its conjugate set $\left\{{\phi }_i\left(z\right)\mathrm{<mo>\; |\; </mo>}i=\mathrm{1,2,}\cdots \mathrm{,8}\right\}$, that is, we may find the appropriate integration method to let ${\displaystyle {\int }_0^H}\frac{1}{\rho }\cdot \overline{{\phi }_i}\cdot {\varphi }_j\text{d}z$ be close to 0

when $i\ne j$. And the numerical results shown in Table 2 verify the cross orthogonality. The real and the imaginary components of the original and the conjugate eigenfunctions from the above modes are shown in Figures 2-5, respectively.

From Figures 2-5, we can see that the eigenfunctions of the propagating modes change relatively stable than the ones of the leaky modes. Furthermore, for leaky modes, the larger the absolute values of the eigenvalues are, the more drastic the change of the eigenfunctions are. And these changes are in accordance with the physical significance. Thus, if there is an eigenfunction corresponding to a leaky mode, it leads to verification difficulty of the cross orthogonality between the eigenfunction $\varphi \left(z\right)$ and its conjugate one $\phi \left(z\right)$, since there is the numerical computation difficulty of the integration for the violent oscillation eigenfunctions.

4. Conclusion

In this paper, firstly, the conjugate operator Q of the Helmholtz operator P with the RBC is constructed. Secondly, it is theoretically proved that there is the cross orthogonality between the linearly-independent eigenfunctions of the operator P

Table 1. Approximate errors of the equation $f\left(\lambda \right)=0$.

Table 2. The orthogonality between different modes.

Figure 2. The real and the imaginary components of original eigenfunctions for propagating modes, where the left and the right figures are represented as the real and the imaginary components of the original eigenfunctions for propagating modes (corresponding to ${\lambda }_1\mathrm{,}{\lambda }_2\mathrm{,}{\lambda }_5$ and ${\lambda }_6$ ), respectively.

Figure 3. The real and the imaginary components of original eigenfunctions for leaky modes, where the left and the right figures are represented as the real and the imaginary components of the original eigenfunctions for leaky modes (corresponding to ${\lambda }_3\mathrm{,}{\lambda }_4\mathrm{,}{\lambda }_7$ and ${\lambda }_8$ ), respectively.

and the linearly-independent eigenfunctions of the operator Q. Finally, Numerical simulation results demonstrate that the cross orthogonality has been almost founded when the RBC is used to truncate the unbounded domain. Namely, it is possible that the high-precision computation of the coordinate transformation under two different local bases is helpful for obtaining more actual propagation behavior. However, we also see that the AGI method has great influence on the

Figure 4. The real and the imaginary components of conjugate eigenfunctions for propagating modes, where the left and the right figures are represented as the real and the imaginary components of the conjugate eigenfunctions for propagating modes (corresponding to ${\lambda }_1\mathrm{,}{\lambda }_2\mathrm{,}{\lambda }_5$ and ${\lambda }_6$ ), respectively.

Figure 5. The real and the imaginary components of conjugate eigenfunctions for leaky modes, where the left and the right figures are represented as the real and the imaginary components of the conjugate eigenfunctions for leaky modes (corresponding to ${\lambda }_3\mathrm{,}{\lambda }_4\mathrm{,}{\lambda }_7$ and ${\lambda }_8$ ), respectively.

accuracy of results, especially for the oscillation eigenfunctions corresponding to the leaky modes with the large absolute values of the eigenvalues. Furthermore, the accuracy of ${\displaystyle {\int }_0^H}\frac{1}{\rho }\cdot \bar{\phi }\cdot \varphi \text{d}z$ is largely dependent on the precision of eigenfunctions.

Therefore, the higher-precision integral method is needed to study in the future.

Acknowledgements

This work was supported by the Ministry of Science and Technology of the People’s Republic of China under Grant 2018YFB1107600.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Lu, Y. and Zhu, J. (2004) A Local Orthogonal Transform for Acoustic Waveguides with an Internal Interface. Journal of Computational Acoustics, 12, 37-53. https://doi.org/10.1142/S0218396X04002183
[2]	Jensen, F., Kuperman, W., Porter, M. and Schmidt, H. (2011) Computational Ocean Acoustics. 2nd Edition, Springer-Verlag, New York. https://doi.org/10.1007/978-1-4419-8678-8
[3]	Zhu, J. and Lu, Y. (2008) Aymptotic Solutions of the Leaky Modes and PML Modes in a Pekeris Waveguide. Wave Motion, 45, 207-216. https://doi.org/10.1016/j.wavemoti.2007.06.001
[4]	Etter, P.C. (2013) Underwater Acoustic Modeling and Simulation. CRC Press, Boca Raton. https://doi.org/10.1201/b13906
[5]	Zhu, J. and Wang, G. (2015) High-Precision Computation of Optical Propagation in Inhomogenneous Waveguides. Journal of the Optical Society of America, 32, 1653-1660. https://doi.org/10.1364/JOSAA.32.001653
[6]	Hu, J. and Menyuk, C.R. (2009) Understanding Leaky Modes: Slab Waveguide Revisited. Advances in Optics and Photonics, 1, 58-106. https://doi.org/10.1364/AOP.1.000058
[7]	Han, H. (2009) Artificial Boundary Method: Numerical Solutions of Partial Differential Equations in Unbounded Domains. Tsinghua University Press, Beijing.
[8]	Hagstrom, T. and Lau, S. (2007) Radiation Boundary Conditions for Maxwell’s Equations: A Review of Accurate Time-Domain Formulations. Journal of Computational Mathematics, 25, 305-336. https://doi.org/10.21236/ADA470448
[9]	Lu, Y.Y. and McLaughlin, J.R. (1996) The Riccati Method for the Helmholtz Equation. The Journal of the Acoustical Society of America, 100, 1432-1446. https://doi.org/10.1121/1.415990
[10]	Shevchenko, V. (1971) The Expansion of the Fields of Open Waveguides in Proper and Improper Modes. Radiophys Quantum Electronics, 14, 972-977. https://doi.org/10.1007/BF01029499