# Coupled Mode Sound Propagation in Inhomogeneous Stratified Waveguides
## 2 Method
### Two-Dimensional Two-Way Coupled Mode Equations
Thus, one can obtain \(\mathrm{H}=C^{\mathrm{T}}\mathbf{A}^{-1}\mathbf{C}\). From Equations (18) and (20), we obtain the first-order coupled mode equations  
\[\frac{\mathrm{d}}{\mathrm{d}x} \begin{pmatrix}\mathbf{P}\\ \mathbf{S}\end{pmatrix} = \begin{pmatrix}-\mathbf{A}^{-1}\mathbf{C}&\mathbf{A}^{ -1}\\ \mathrm{L}-\mathrm{K}&C^{\mathrm{T}}\mathbf{A}^{-1}\end{pmatrix}\begin{pmatrix} \mathbf{P}\\ \mathbf{S}\end{pmatrix} . \tag{22}\]  
We note that, to calculate acoustic fields, only four of the nine matrices in Equation (8) remain to be computed. The first-order coupled system Equation (22) preserves the energy conservation property but does not contain any interface matrix. For details of energy conservation of Equation (22), see Appendix B. It is the proper application of \(\mathbf{S}\) that gives the coupled mode equations a simple form and substantially reduces computing power, especially in long-range application. The first-order differential equations account for source and radiation conditions. Directly integrating Equation (22) starting from the radiation condition may encounter the problem of numerical divergence owing to the exponential growth of evanescent modes. In order to avoid numerical problems, we introduce admittance matrix \(\mathbf{Y}\) and propagator matrix \(\mathbf{M}\) by \(\mathbf{S}(x)=\mathbf{Y}(x)\mathbf{P}(x)\) and \(\mathbf{P}(x)=\mathbf{M}(x)\mathbf{P}(0)\), respectively. \(\mathbf{Y}\) satisfies a Riccati equation and \(\mathbf{M}\) a first-order differential equation coupled to \(\mathbf{Y}\):  
\[\mathbf{Y}^{\prime}=-\mathbf{Y}\mathbf{A}^{-1}\mathbf{Y}+\mathbf{Y}\mathbf{A}^ {-1}\mathbf{C}+C^{\mathrm{T}}\mathbf{A}^{-1}\mathbf{Y}+\mathrm{L}-\mathrm{K}, \tag{23}\]  
\[\mathbf{M}^{\prime}=-\mathbf{A}^{-1}\mathbf{C}\mathbf{M}+\mathbf{A}^{-1} \mathbf{Y}\mathbf{M}. \tag{24}\]  
\(\mathbf{Y}\) can be computed by integrating from right to left with the initial value \(\mathbf{Y}(x_{R})=\sqrt{\mathbf{A}(x_{R})\Big{(}\mathrm{L}(x_{R})-\mathrm{K}( x_{R})\Big{)}}\). Eigenvalues of matrix \(\mathbf{A}^{-1}(x_{R})\mathbf{Y}(x_{R})\) are equal to \(\mathrm{i}k_{x}\), where \(k_{x}\) are horizontal wavenumbers in region \(x>x_{R}\) (for more details, see Appendix C). \(\mathbf{M}(x)\) can be solved from left to right with the initial value \(\mathbf{M}(0)=\mathrm{l}\), where \(\mathrm{l}\) is the identity matrix. In this paper, numerical integration of Equations (23) and (24) is performed using the Magnus method [25; 26]. Further, the reflection and transmission matrices that characterize the scattering properties of the scattering region can be computed. Reflection matrix \(\mathrm{R}\) is defined by \(\mathbf{P}_{r}(0)=\mathrm{RP}_{i}(0)\), and transmission matrix \(\mathrm{T}\) is defined by \(\mathbf{P}_{i}(x_{R})=\mathrm{TP}_{i}(0)\), where \(\mathbf{P}_{i}\), \(\mathbf{P}_{r}\), and \(\mathbf{P}_{t}\) denote the modal amplitudes of the incident, reflected and transmitted waves, respectively. We obtain  
\[\mathrm{R}=\Big{[}\mathbf{Y}_{0}+\mathbf{Y}(0)\Big{]}^{-1}\Big{[}\mathbf{Y}_{0 }-\mathbf{Y}(0)\Big{]}, \tag{25}\]  
and  
\[\mathrm{T}=\mathbf{M}(x_{R})(\mathrm{l}+\mathrm{R}), \tag{26}\]where \(Y(0)\) is the admittance matrix at \(x=0\), and \(Y_{0}=\sqrt{A(0)\left(L(0)-K(0)\right)}\) is the local admittance matrix at \(x=0\). Eventually, one can write the acoustic pressure in the multilayered waveguide as  
\[p(x,z)= \begin{cases}\mathbf{\Psi}^{\mathrm{T}}(z;x)M(x)(1+R)\mathbf{\mathrm{P}}_{i}(0), &0\leq x\leq x_{R},\\ \mathbf{\Psi}^{\mathrm{T}}(z;x_{R})\mathrm{e}^{A^{-1}(x_{R})Y(x_{R})(x-x_{R})} \mathbf{\mathrm{T}}\mathbf{\mathrm{P}}_{i}(0),&x>x_{R},\end{cases}  \tag{27}\]  
where \(\mathbf{\mathrm{T}}\) is a vector with elements being \(\psi_{n}\). To summarize, once the geometry and medium parameters of the waveguide are identified, admittance matrix \(Y\) and propagator matrix \(M\) are available for all \(x\). Then, reflection matrix \(R\) and transmission matrix \(T\) can be deduced. Finally, the sound field generated by any incident wave can be obtained.