It is known that the computation of the wave propagation plays a particularly prominent role in ocean acoustics or optical waveguides   . 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 .
The propagation of sound waves in the medium is mathematically formulated as a well-known Helmholtz equation . 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 , 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 . 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 . 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)  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 .
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:
Where z is the depth axis pointing downward with the ocean surface at , and is the range variable in horizontal direction. One layer with density is located in , the other with density is located in ; and their interface is flat at . Here, and are represented as the wavenumbers in two different subregions, respectively.
We introduce the RBC for truncation at , that is at , where , , , and is the eigenvalue of the Helmholtz operator . Here , where is as and is as . Then Equation (1) is approximately transformed into a complex partial differential equation in the bounded region as follows:
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 . If we assume the solution of Equation (2) to be , then the eigenvalue problem of the Equation (2) is derived as follows:
Where is the eigenfunction corresponding to the eigenvalue for the eigenvalue problem of Equation (2).
Given , if the solution functions , of Equation (3) are crossly orthogonal to their conjugate functions , with the weight , that is, as , then we have , where each eigenfunction corresponds to each eigenvalue , respectively, ; and the overline represents the conjugate operation.
The following processes focus on constructing the conjugate eigenfunctions , corresponding to the eigenfunctions , , and further prove their cross orthogonality.
For Equation (3), we have the general solutions as follows:
By the RBC at , namely , we can obtain the following equation of the eigenvalue :
2.2. The Construction of the Conjugate Operator
Let the operator be with the RBC at and the interface conditions at , we have
Let and be the eigenfunction and the corresponding conjugate eigenfunction of the operator , respectively, where the boundary conditions are
given as .
Then, by the integration by parts and making use of the interface and the boundary conditions, we have
where is required.
Also, we have
Thus, we have
Applying the interface and the boundary conditions:
the above expression is simplified as follows:
For simplicity, the interface and the boundary conditions are given to the conjugate eigenfunction as follows:
it is equivalent to
it is equivalent to
it is equivalent to
Thus, by the interface and the boundary conditions of the function as above, then we have
Finally, we have
In fact, there is the following relation between the operators P and Q:
Namely, we have
Obviously, is the conjugate function of , that is also called the conjugate eigenfunction of the eigenfunction . In this paper, we assume that is a piecewise constant function. Under above assumptions, should satisfy the following ordinary differential equation:
We have the similar solution form of Equation (27) as Equation (3):
where the parameters and are to be determined.
In general, we choose , then we can write down the solution as
Making , we can get the same nonlinear equation of the eigenvalue 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 and their linearly-independent conjugate eigenfunctions , .
Consider the eigenvalues and eigenfunctions of both and , which satisfy
This give rises to
This is because
At last, we also have the following conclusion.
Theorem 2: For the linearly-independent eigenfunctions and their linearly-independent conjugate eigenfunctions , , as mentioned above, then there is the cross orthogonality, that is, , .
3. Numerical Examples
In this section, we will verify the method presented in the previous sections. Considering a Pekeris waveguide, we let . These parameters are taken from reference . By the nonlinear Equation (6), we define
We find out the roots of the equation 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 . The details of choosing initial values are listed in the reference . As a result, we obtain the eigenvalues’ distributions of the operators and , 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:
Figure 1. The original eigenvalues of the operator represented by “o” and the eigenvalues of the operator represented by “+”.
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
and its conjugate set , that is, we may find the appropriate integration method to let be close to 0
when . 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 and its conjugate one , since there is the numerical computation difficulty of the integration for the violent oscillation eigenfunctions.
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 .
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 and ), 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 and ), 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 and ), 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 and ), 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 is largely dependent on the precision of eigenfunctions.
Therefore, the higher-precision integral method is needed to study in the future.
This work was supported by the Ministry of Science and Technology of the People’s Republic of China under Grant 2018YFB1107600.
 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.
 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.