Two Chebyshev Spectral Methods for Solving Normal Modes in Atmospheric Acoustics

The normal mode model is important in computational atmospheric acoustics. It is often used to compute the atmospheric acoustic field under a time-independent single-frequency sound source. Its solution consists of a set of discrete modes radiating into the upper atmosphere, usually related to the continuous spectrum. In this article, we present two spectral methods, the Chebyshev-Tau and Chebyshev-Collocation methods, to solve for the atmospheric acoustic normal modes, and corresponding programs are developed. The two spectral methods successfully transform the problem of searching for the modal wavenumbers in the complex plane into a simple dense matrix eigenvalue problem by projecting the governing equation onto a set of orthogonal bases, which can be easily solved through linear algebra methods. After the eigenvalues and eigenvectors are obtained, the horizontal wavenumbers and their corresponding modes can be obtained with simple processing. Numerical experiments were examined for both downwind and upwind conditions to verify the effectiveness of the methods. The running time data indicated that both spectral methods proposed in this article are faster than the Legendre-Galerkin spectral method proposed previously.


Introduction
The propagation of sound waves in the atmosphere is a basic subject of atmospheric acoustics [1]. Sound waves in the atmosphere undergo a series of complex processes, including ground reflection, atmospheric scattering, refraction, and absorption [2]. In fact, the propagation of sound waves in the atmosphere satisfies the wave equation, but it is difficult to strictly solve the wave equation. Thus, scientists make approximations to the wave equation for specific situations, thereby obtaining easy-to-solve equations, which can be solved numerically to obtain a solution of the sound field. Numerical sound fields have the advantages of intuitiveness and clarity, and they are widely used in acoustic research. Based on this idea of solving the numerical sound field, computational atmospheric acoustics, a sub-discipline of atmospheric acoustics, has been developed. Numerical models have many forms. Different models are suitable for different environments, and the results are not exactly the same. Mainstream numerical models include the parabolic Equation (PE) [3,4] model, the wavenumber integration method (the fast field program (FFP)) [5][6][7], and ray [8] and Gauss beam [2] approaches. The normal mode model is also a fundamental method for solving for the acoustic field in the atmosphere with a finite ground impedance and horizontally stratified sound speed [1,2,5]. A horizontally stratified atmosphere allows the wave equation to be solved by the separation-of-variables method. After using Hankel integral transforms, the sound field can be expressed in terms of the sum of normal modes [9,10]. When the ground impedance is complex or there is sound attenuation in the atmosphere, it is complicated to use the finite difference method to solve for the atmospheric normal modes, and the result is not very accurate [2].
In recent years, progress has been made on using spectral methods to solve underwater acoustic problems, and small-scale research has begun to link the spectral methods with the normal modes of underwater acoustics. Dzieciuch [11] developed MATLAB code for computing normal modes based on Chebyshev approximations. Although he only realized the calculation of the simple Munk waveguide, this was the first step in the application of the spectral methods to computational ocean acoustics. In 2016, Evans [12] used the Legendre-Galerkin spectral method to develop a sound propagation calculation program in a layered ocean environment. Subsequently, Tu et al. [13][14][15] used the Chebyshev-Tau spectral method to develop a program for calculating sound propagation in single-layer and layered ocean environments. They subsequently solved for the normal modes in underwater acoustics using the Legendre-Collocation method and proved that both of the spectral methods had high accuracy [16]. They also applied the spectral methods to the parabolic approximation of underwater acoustics [13,17,18]. The results of these studies indicated that it is feasible to apply spectral methods for the calculation of underwater acoustics, and in many cases, it has higher accuracy than the finite difference method. Monographs on spectral methods have also confirmed this [19][20][21][22]. Throughout the history of the development of atmospheric acoustics, many methods in underwater acoustics have been introduced [1,2]. In computational atmospheric acoustics, spectral methods are rarely used to calculate the numerical sound field. In 2017, Evans [23] successfully introduced the Legendre-Galerkin spectral method to construct atmospheric acoustic normal modes. He then further improved the method [24] and proved the convergence of the method [25]. In 2019, Sabatini [26] used the Chebyshev-Collocation method to discretize the governing equation of atmospheric acoustics into quadratic eigenvalue problems.
In this article, we propose two spectral methods for calculating atmospheric acoustic normal modes. The results are compared with Evans's code [24], the correctness of the two spectral methods proposed in this article was verified, and computational speeds of the two spectral methods were demonstrated. The text is organized as follows. Section 2 describes normal modes in the atmosphere mathematically. Section 3 provides brief descriptions of the Chebyshev-Tau and Chebyshev-Collocation spectral methods and introduces the discretization of atmospheric acoustic normal modes. In Section 4, two numerical experiments are shown to verify the correctness of the methods proposed in this article. Section 5 analyzes the running speed of the spectral methods, and Section 6 concludes this article.

Atmospheric Normal Modes
Acoustic theory reveals that the core of solving the acoustic field with a time-independent single-frequency sound source is the following wave Equation [5]: In the above homogeneous Helmholtz equation, p is the sound pressure in the frequency domain to be solved, and k represents the wavenumber, which is related to the frequency of the source and the spatial position. p and k are the functions of the spatial position, i.e., p(x, y, z) and k(x, y, z), respectively. We consider the medium of sound propagation to be the atmosphere depicted in Figure 1.
The ∆ operator in Equation (1) in cylindrical coordinates is taken to obtain the acoustic governing equation in the cylindrical coordinate system (r, z), where r is the range, and z is the depth. Considering the case in Figure 1 where the wavenumber k(z) is only related to depth (range-independent), Equation (1) where k(z) = ω/c(z), ω = 2π f is the angular frequency of the sound source, f is the frequency of the source, and c(z) is the sound speed profile. When considering the atten-uation of sound waves by the atmosphere, k(z) = [1 + iηα(z)]ω/c(z), where α(z) is the attenuating coefficient in units of dB per wavelength, and η = (40π log 10 e) −1 . Through separation of variables, the acoustic pressure p(r, z) can be decomposed as follows: where R(r) can be approximated by the analytical form of the Hankel function of the first kind H (1) 0 (k r m r) and ψ(z) satisfies the following modal equation: The modal equation is a Sturm-Liouville equation, and its characteristics are well known; that is, after adding appropriate boundary conditions, it has a series of modal solutions (k r , ψ), where k 2 r is a constant. When the considered medium has attenuation, k(z) is a complex function. The lower boundary of the atmosphere is the ground, and sound waves on the ground usually need to meet the following impedance boundary conditions: where Z is the normalized ground impedance. The upper boundary of the atmosphere can be regarded as a free boundary at infinity, or it can be called an acoustic half-space condition. In 2019, Sabatini [26] used the Chebyshev-Collocation method to discretize the governing equation with a half-space condition into quadratic eigenvalue problems, which were mathematically solved perfectly. However, in actual calculations, the complexity of this solution process is relatively high. In addition, the solution of the quadratic eigenvalue problem is inherently unstable [27][28][29].
Solving the standard Sturm-Liouville problem will yield multiple sets of solutions (k r m , ψ m ), m = 1, 2, . . . , where k r m is called the horizontal wavenumber and ψ m is called the eigenmode or mode. The modes of Equation (4) are arbitrary up to a nonzero scaling constant, so they should be normalized [5] as follows: Finally, the fundamental solution to the acoustic governing Equation (2) in the atmosphere can be approximated as follows [23]: where M is the number of modes used to synthesize the sound field and the Hankel function H 0 (k r m r) related to R(r) adopts its asymptotic form (see Equations (2.39), (5.13), and (5.14) of Reference [5]). The core of solving for the normal modes of atmospheric acoustics is the solution of the differential equations in Equations (4)- (6). Solving for the normal modes of the atmospheric acoustics requires the discretization of Equations (4)- (6). Traditionally, the domain of the problem solved by the spectral method is usually in the interval [−1, 1], so we first use 1]. Noting that dz/dx = H/2, we let the operators L, P, Q have the following forms: Equations (4)-(6) can be written in the following form: Next, we develop two spectral methods to solve this system.

Discretized Atmospheric Normal Modes by Two Spectral Methods
A spectral method is a kind of weighted residual method, and it can provide accurate solutions to differential equations [20,21]. In the spectral method, the unknown function to be solved ψ(x) is expanded by a set of linearly independent bases φ k (x). When the number of bases tends to infinity, an accurate representation of ψ(x) can be obtained. However, in actual calculations, it is usually necessary to truncate to the first N-order terms, thus obtaining an approximation of ψ(x) as follows: whereψ k represents the expansion coefficients. Obtaining the value ofψ k is equivalent to obtaining the approximate solution ψ N (x) of ψ(x). Inserting ψ N (x) from Equation (10) into Equation (9), Equation (9) is no longer strictly true, and there is a residual Res(x), defined as follows: To make ψ N (x) as close to ψ(x) as possible, we need to minimize the residual through a certain principle [22]. Setting the weighted integral of the residuals equal to zero is a widely used principle [20]: From Equation (10), the residual Res(x) can be minimized only by adjusting the value of the expansion coefficientsψ k . The choice of the weight function w(x) is also crucial.
In the two spectral methods developed in this article, the basis functions φ k (x) are both Chebyshev polynomials T k (x), and the difference is the selection of weight functions. The Chebyshev polynomial basis functions are provided in [14,[20][21][22].

Discretized Atmospheric Normal Modes by Chebyshev-Tau Spectral Method
In the Chebyshev-Tau spectral method, in addition to the basis functions being Chebyshev polynomials (φ k (x) = T k (x)), the weight functions are also Chebyshev polynomials (w(x) = T k (x)). Inserting Equations (10) and (11) into Equation (12), we obtain the new form of Equation (12) for the Chebyshev-Tau spectral method: where 1 √ 1−x 2 is the orthogonal weighting factor of the Chebyshev polynomial basis function space. This equation is also known as the weak form of Equation (9). It will form (N − 1) algebraic equations (excluding the boundaries of x), the two boundary conditions will produce two algebraic equations, and the unknowns to be solved for areψ 0 toψ N . The integral formulas listed in the above equations can be computed by the Gauss-Chebyshev-Lobatto quadrature [21] to obtain accurate results. To include the two end points of the domain x, the Chebyshev-Gauss-Lobatto nodes on x ∈ [−1, 1] are taken [21]: There are two forms of Gauss-Chebyshev-Lobatto quadrature [21]: where f (x) is the function to be integrated. In fact, using Chebyshev-Gauss-Lobatto nodes has additional advantages. As shown in Figure 2, Chebyshev-Gauss-Lobatto nodes are dense at both ends and sparse in the middle, so the discrete points are more densely distributed near the ground [0, h]. Owing to the interaction of temperature, humidity, turbulence, and other factors in the atmosphere near the ground, acoustic profiles are usually more complex and changeable. Figure 3 shows the entropy of the 2018 annual sound speed profile in Guangzhou, China (22.5 • N, 112.9 • E), calculated using the ERA-40 reanalysis data provided by the European Centre for Medium-Range Weather Forecasts (ECMWF). It can be seen from the figure that the entropy of sound speed near the ground is much higher than that at the upper air. Therefore, the more densely distributed discrete points of Chebyshev-Gauss-Lobatto nodes near the ground are conducive to accurately simulating the sound field.
Chebyshev-Gauss-Lobatto nodes We convert the original solution of the unknown function ψ(x) into solving for its expansion coefficientsψ k under the Chebyshev polynomial basis. The only difficulty is the discretization of the operator L. The conclusions used in the following text are directly given here. For detailed derivations, readers can refer to References [19][20][21][22].
The derivative d dx is included in the L operator, and the expanded coefficientsψ k of ψ (x) satisfy the following relationship withψ k : where " " denotes derivative, and the second formula is the vector form of the first formula, where column vectorsψ = [ψ 0 , · · · ,ψ N ] T andψ = [ψ 0 , · · · ,ψ N ] T , respectively.D is a square matrix of order (N + 1). To distinguish it from the differential matrix D in the Chebyshev-Collocation spectral method, a hat symbol is added to the relationship matrix. The known function k 2 (x) is included in the L operator. Letting v(x) = k 2 (x), there will be a product term y(x) = v(x)ψ(x), and the expanded coefficientsŷ k of y(x) satisfy the following relationship withψ k : Similarly, the second formula is the equivalent vector form, andĈ v is also a square matrix of order (N + 1), and the subscript v indicates that the known function in the operator is v(x).
We show the discretization of the operator L. In the Chebyshev-Tau spectral method, applying Equations (16) and (17) to Equation (9), we can obtain the discrete forms of the operator L and Equation (9) as follows: The boundary conditions produce algebraic equations about the expansion coefficients in the Chebyshev-Tau spectral method as follows. In the Tau method, the function ψ(x = ±1) in the boundary conditions is also expanded by Equation (10). The discretization of the boundary operators P and Q is similar to that of operator L, so the two boundary conditions generate two equations related toψ k . To facilitate the description of the processing of the boundary conditions, the following intermediate row vectors are defined as follows:t The matrix form of the discrete ground and air boundary conditions Equations (5) and (6) in the Chebyshev-Tau spectral method can be written as follows: The algebraic equations formed by these two boundary conditions and the (N − 1) algebraic equations obtained from the weak form are solved simultaneously, and we can then solve forψ k and obtain ψ(x). The row vectorsp andq are used to replace the last two rows of theL matrix in Equation (20), and the last two elements of the column vectorψ on the right-hand side of Equation (18) are replaced with 0, so that the boundary conditions are strictly met. We let the matrix composed of the first (N − 1) rows and (N − 1) columns ofL beL 11 . The matrix composed of the first (N − 1) rows and the last two columns ofL isL 12 . The row vectors composed of the first (N − 1) elements of the row vectorsp,q, andψ arep 1 ,q 1 , andψ 1 , respectively. The row vectors composed of the last two elements of the row vectorsp,q, andψ arep 2 ,q 2 , andψ 2 , respectively. Thus, Equation (18) can be changed to the following block form: According to the horizontal and vertical lines in the above formula, Equation (20) can be abbreviated as follows: Equation (21) can be solved as follows: Therefore, a set of (k 2 r ,ψ) can be solved for by the (N − 1)th-order matrix eigenvalue problem in Equation (22). For each set of eigenvalues/eigenvectors (k 2 r m ,ψ m ), an eigensolution (k r m , ψ m ) of Equation (4) can be obtained by Equation (10). In this process, each eigenmode should be normalized by Equation (7). Finally, the sound pressure field is obtained by applying Equation (8) to the chosen modes.

Discretized Atmospheric Normal Modes by Chebyshev-Collocation Spectral Method
The Collocation method uses the Dirac function δ(x) as the weight function in Equation (12). The characteristics of the δ(x) function are well known. In the Collocation method, Equation (12) becomes the following: The above formula shows that in the Collocation method, the weighted residual principle becomes that the residuals are all 0 at the selected discrete points x j . Its essence is to only make the original differential Equation (9) strictly hold on this set of discrete points, so as to solve for the function value ψ(x j ) of the modal function ψ(x) on this set of discrete points as an approximation. In the Collocation method, there is no need to expand the function to be sought as Equation (10). This is why the Collocation method is considered to be a special spectral method, sometimes called the pseudospectral method [21]. In the Chebyshev-Collocation method, we also take the discrete points of the Chebyshev-Gauss-Lobatto nodes in Equation (14). In this case, the only difficulty is the discretization of operator L. The conclusions used in the following text are directly given here as in the introduction of the Chebyshev-Tau spectral method. For a detailed derivation, readers can refer to References [19][20][21][22]. The derivative terms ψ (x) and ψ(x) have the following relationship: where ψ = [ψ (x 0 ), ψ (x 1 ), ψ (x 2 ), · · · , ψ (x N )] T represents the function value of the derivative term ψ (x). Similarly, ψ = [ψ(x 0 ), ψ(x 1 ), ψ(x 2 ), · · · , ψ(x N )] T . Matrix D is also called the Chebyshev-Collocation differential matrix. The product y(x) = v(x)ψ(x) can be processed as follows: where For the Collocation method, the boundary conditions are only related to the endpoints of the domain x, so the discrete points on the boundaries (x 0 and x N ) only need to satisfy the boundary conditions, not the differential equation. The discretized forms of the operators P and Q are similar to that of operator L. Similar to the Chebyshev-Tau spectral method, the operator L also needs to be discretized in the Chebyshev-Collocation method. With reference to Equations (24) and (25), in the Chebyshev-Collocation method, the L operator and Equation (9) have the following forms: To facilitate the description of the processing of the boundary conditions, the first and last rows of D are defined as row vectors d 1 and d 2 , respectively. The following intermediate (N + 1)-dimensional row vectors are defined as follows: The matrix form of the discrete ground and air impedance conditions Equations (5) and (6) in the Chebyshev-Collocation method can be written as follows: In the Collocation method, the row vectors p and q are used to replace the first row and the last row of the L matrix in Equation (26), so that the boundary conditions are satisfied. We let the block matrix formed by the second row to the N-th row of the matrix L be L 1 , and the column formed by the second to the N-th elements of ψ be ψ 1 . Equation (26) can then be written as follows: We only need to perform a simple row transformation and column transformation on Equation (28) to transform it to a form similar to Equation (21), and then we use the same method used for Equation (22) to find the eigenvalues/eigenvectors. A set of (k 2 r , ψ) can be solved by the (N − 1)th-order matrix eigenvalue problem in Equation (22). Each eigenmode should be normalized by Equation (7). Finally, the sound pressure field is obtained by applying Equation (8) to the chosen modes.

Numerical Experiment and Analysis
To verify the correctness of the spectral methods presented above in solving the normal modes of atmospheric acoustics, the authors developed the corresponding programs based on the above derivation. The programs based on Chebyshev-Tau spectral method and Chebyshev-Collocation method are called "AtmosCTSM" and "AtmosCCSM", respectively. The code was written in FORTRAN/MATLAB and is available at the author's GitHub homepage (https://github.com/tuhouwang/Atmospheric-normal-modes, accessed on 4 May 2021). For comparison, we considered the program "aaLG" based on the Legendre-Galerkin spectral method, which was developed by Richard B. Evans in FORTRAN and verified by comparison with PE and FFP [24].
Sound propagation is sensitive to the atmospheric state, particularly to the temperature and wind. Thus, to model sound propagation, one must know the state of the atmosphere at the time of propagation [30]. The two examples shown by Evans [23] can be used as benchmark examples. The source frequency of both cases was f = 100 Hz at a height of z s = 5 m above the ground. The normalized ground impedance Z (related to the constant G in Equation (6)) is the same as the value used by Gilbert [4], and the value is Z = 12.97 + 12.38i. In the following two experiments, the order of the spectral truncation in the three spectral methods was taken as N = 1500. Using the TL to express the acoustic field [5], the relationship between it and the sound pressure is TL = −20 log 10 (|p|/|p 0 |), where p 0 is the sound pressure at a distance of 1 m from the sound source.

Downwind Case
The first numerical experiment was a downwind case. The piecewise linear acoustic parameter profile used in this numerical experiment, which was presented by Evans [23], is shown in Table 1. The table clearly reveals that the thickness of the atmosphere is 700 m, and the artificial absorber layer is located between 700 and 2000 m. Table 1. Piecewise linear acoustic parameter profile used in experiment 1, cited from Evans [23].  Figure 4 shows the horizontal wavenumbers k r calculated by the Legendre-Galerkin spectral method and the two spectral methods developed in this article on the complex plane. The consistency of the eigenvalue distribution in the figure illustrates the correctness of the horizontal wavenumbers calculated by the three methods. Figure 5 shows the first four normal modes of experiment 1. It reveals that the modes obtained by the two spectral methods proposed in this article were highly consistent with those obtained by the Legendre-Galerkin spectral method. Figure 6 presents an overview of the acoustic fields obtained by the three methods. We used the first 552 modes with phase velocities between 341.7 and 391.2 m/s to synthesize the sound fields. The horizontal wavenumbers of these modes are shown in Figure 4. The acoustic fields calculated by the three methods were very similar. Figure 7 shows the TL curves versus the range for a receiver at a height of 1 m over the range interval 0-5 km. The results of the two spectral methods presented in this article were very similar to those of the Legendre-Galerkin spectral method, and there may have been small differences only in the acoustic shadow areas.

Upwind Case
The second numerical experiment is an upwind case. The piecewise linear acoustic parameter profile used in this numerical experiment was presented by Evans [23], and it is shown in Table 2. The table clearly illustrates that the thickness of the atmosphere was 900 m, and the artificial absorber layer was located between 900 and 2000 m. Table 2. Piecewise linear acoustic parameter profile used in experiment 2, cited from Evans [23].  Figure 9 presents the acoustic fields calculated by the three spectral methods, where 553 modes with phase velocities less than 393.2 m/s were used to synthesize the sound field. In the atmosphere layer, the sound fields calculated by the three methods were highly consistent. Figure 10 shows the TL curves versus the range for a receiver at a height of 1 m over the range interval 0-10 km from the Legendre-Galerkin spectral method; the figure shows that the results of several methods were very similar. The differences between the three methods are indistinguishable at this plotting accuracy, and there may have been small differences only in the acoustic shadow area.  From the numerical results displayed above, we can see that three methods with different theoretical foundations all yielded very similar acoustic fields and normal modes, regardless of whether the sound speed profile was downwind or upwind. The consistency of these three methods proved that the two spectral methods proposed in this article are feasible for solving the atmospheric normal modes.

Discussion of Computational Speed
To further compare the characteristics of the two spectral methods proposed in this article, we divided each method into four steps, and we discuss the running time and complexity of each part separately. The four steps are as follows: discretizing the equation, solving eigenvalue problems, obtaining normal modes, and synthesizing the sound field. Table 3 lists the time consumption of each step of the programs in the two experiments. The time listed in the table is the average of ten tests. In the tests, the three programs were run on a Dell XPS 8930 desktop computer equipped with an Intel i7-8700K CPU. The FORTRAN compiler used in the test was gfortran 7.5.0.
In terms of speed, the AtmosCCSM was slightly faster than the AtmosCTSM. This is because the Tau method requires forward and backward Chebyshev transformations, unlike the Collocation method. The aaLG program was much slower than the two newly developed programs. Solving the eigenvalue problem was the most time-consuming step for the three programs. Moreover, the aaLG program spent much more time than the other two programs on solving the eigenvalue problem. However, the aaLG uses a subroutine developed by Evans [24,31] to solve the matrix eigenvalue problem, while both the AtmosCTSM and AtmosCCSM solve the eigenvalue problem by calling the Lapack numerical library. In fact, matrix eigenvalue problems have the same computational complexity O(N 3 ). It is apparent from this table that the subroutine written by Evans is much slower than the Lapack numerical library, which is the main reason that the aaLG consumed much more time than the other two programs. We modified aaLG to also call the Lapack numerical library when solving the eigenvalue problem. The modified program was named "aaLG-M". The fourth column of Table 3 lists the running time of aaLG-M. The aaLG-M took roughly the same time to solve the matrix eigenvalue problem as the other two programs. However, the aaLG-M was still slower than the other two programs. The most significant difference in the running time between the three programs was in the first step (discretizing the equation). In the first step, each element of the matrix finally obtained by the Legendre-Galerkin spectral method must be numerically integrated for every piecewise linear acoustic profiles, which means the number of calculations is very large, and its computational complexity is O(N 2 ). In contrast, the two methods proposed in this article only need to perform simple interpolation of the acoustic profiles and matrix multiplication to obtain the discrete equations; the computational complexity of linear interpolation is O(N). In the mode-obtaining and sound-field-synthesizing steps, the two spectral methods devised in this article still required less time than the Legendre-Galerkin method; part of the reason for this is that the authors used certain skills in programming to optimize them. It is worth mentioning that the AtmosCTSM.m and AtmosCCSM.m programs (developed in MATLAB, which is better at matrix operations) could obtain the results of the above experiments in less than 4 s (run on the same platform in MATLAB 2019a), which is an attractive result.

Conclusions
In this article, we propose two spectral methods for solving for atmospheric acoustic normal modes. An artificial absorber layer was added above the atmosphere of interest to reduce the impact of the truncated half-space on the area of interest. Next, we designed two examples, performed a detailed analysis of the results of each example, and finally verified the correctness and reliability of the proposed methods. Tests on the running time of the programs developed based on the three spectral methods showed that, in terms of the running time, the methods proposed in this article had better speeds than the Legendre-Galerkin spectral method.