Random matrix theory for low-frequency sound propagation in the ocean: a spectral statistics test

TL;DR

This paper evaluates a novel random matrix theory approach for modeling low-frequency sound propagation in the deep ocean, demonstrating its effectiveness in reproducing spectral statistics and mode coupling over long distances.

Contribution

It introduces and tests a RMT-based propagator for ocean sound propagation, showing its ability to match numerical solutions and capture spectral features.

Findings

RMT-based propagator aligns well with direct numerical solutions for long-range propagation.

Mode coupling is effectively modeled by the RMT approach.

Agreement decreases at shorter distances due to cross-mode correlations.

Abstract

Problem of long-range sound propagation in the randomly-inhomogeneous deep ocean is considered. We examine a novel approach for modeling of wave propagation, developed by K.C.Hegewisch and S.Tomsovic. This approach relies on construction of a wavefield propagator using the random matrix theory (RMT). We study the ability of the RMT-based propagator to reproduce properties of the propagator corresponding to direct numerical solution of the parabolic equation. It is shown that mode coupling described by the RMT-based propagator is basically consistent with the direct Monte-Carlo simulation. The agreement is worsened only for relatively short distances, when long-lasting cross-mode correlations are significant. It is shown that the RMT-based propagator with properly chosen range step can reproduce some coherent features in spectral statistics.

Equations

$$c(z,r)=c_0+\Delta c(z)+\delta c(z,r),$$

$$|\delta c{|}_{\text{max}}\ll |\Delta c{|}_{\text{max}}\ll c_0.$$

$$\frac{i}{k_0}\frac{\partial \Psi }{\partial r}=-\frac{1}{2k_0^2}\frac{{\partial }^2\Psi }{\partial z^2}+\left[U(z)+V(r,z)\right]\Psi ,$$

$$k_0=\frac{2\pi f}{c_0}.$$

$$U(z)=\frac{\Delta c(z)}{c_0},V(r,z)=\frac{\delta c(r,z)}{c_0}.$$

$${\Psi |}_{z=0}=0,{\frac{d\Psi }{dz}|}_{z=h}=0.$$

$$\Psi (r,z,k_0)=\sum _m{a}_m(r,k_0){\psi }_m(z)$$

$$-\frac{1}{2k_0^2}\frac{{\partial }^2{\psi }_m(z)}{\partial z^2}+U(z){\psi }_m(z)=E_m{\psi }_m(z).$$

$$a_m(r,k_0)=a_m(r=0,k_0)e^{-i{k}_0{E}_{m}r},$$

$$\Psi (r_{\mathrm{f}},z)=\hat{G}(r_0,r_{\mathrm{f}})\Psi (r_0,z).$$

$$G_{mn}(r_0,r_{\mathrm{f}})=\int {\psi }_m^{\ast }\hat{G}(r_0,r_{\mathrm{f}}){\psi }_{n}dz,$$

$$\hat{G}(r_0,r_J)=\hat{G}(r_{J-1},r_J)\hat{G}(r_{J-2},r_{J-1}\text{)}\text{...}\hat{G}(r_1,r_2)\hat{G}(r_0,r_1),r_J=r_{\mathrm{f}}.$$

$$\hat{G}(0,r_{\mathrm{f}})=\hat{G}(J\Delta r)=\prod _{j=1}^J{\hat{G}}_j(\Delta r),$$

$$\mathbf{G}(\Delta r)=\mathbf{\Lambda }(\mathbf{I}-i\mathbf{A}),$$

$${\Lambda }_{mn}={\delta }_{mn}{e}^{-i{k}_0{E}_m\Delta r},$$

$$A_{mn}=k_0{e}^{i{k}_0(E_m-E_n)r^{\prime }}\underset{r^{\prime }=0}{\overset{\Delta r}{\int }}{V}_{mn}(r^{\prime })d{r}^{\prime },$$

$$V_{mn}(r)=\int {\psi }_m^{\ast }(z)V(r,z){\psi }_n(z)dz.$$

$$A_{mn}(\Delta r,k_0)={\sigma }_{mn}(\Delta r,k_0)z_{mn}(k_0),$$

$$\mathbf{G}(\Delta r)=\mathbf{\Lambda }[\mathbf{I}+i\mathbf{A}(\Delta r)/2{]}^{-1}[\mathbf{I}-i\mathbf{A}(\Delta r)/2]$$

$$\vec{a}(r)=\mathbf{G}(r)\vec{a}(0).$$

$$c_{\mathrm{b}}(z)=c_0+\Delta c(z)=c_0\left[1+\frac{b^2}{2}(e^{-az}-\kappa {)}^2\right],$$

$$c_{\text{b}}(z)=c_0\left[1+\gamma (e^{-\eta }-1+\eta )\right],\eta =\frac{2(z-z_{\text{a}})}{B},$$

$$\delta c(r,z)=c_0\left(\frac{24.5}{g}\right)N^2\zeta ,$$

$$\zeta (r,z)=\frac{2B}{\pi }\sqrt{\frac{E\Delta k_l}{M}}{e}^{-3z/2B}\sum _{j=1}^{j_{\max }}\sum _{l=1}^{l_{\max }}\sin (j\pi \xi (z))\sqrt{\frac{I_{j,k_l}}{j^2+j_{\ast }^2}}\cos (k_{l}r+{\varphi }_{jl}),$$

$$I(j,k_l)=\frac{k_j}{k_l^2+k_j^2}+\frac{1}{2}\frac{k_l^2}{(k_l^2+k_j^2{)}^{3/2}}\text{ln}\frac{\sqrt{k_l^2+k_j^2}+k_j}{\sqrt{k_l^2+k_j^2}-k_j},$$

$$k_j=\frac{\pi j{f}_{\text{i}}}{N_{0}B}.$$

$$\delta c_{\mathrm{iw}}(z,r)=<\delta c_{\mathrm{iw}}(z)>+\sum _k{q}_k(r){\theta }_k(z).$$

$$K_{ij}=\frac{1}{l_{\text{max}}}\sum _{l=1}^{l_{\text{max}}}[\delta c_l(z_i)-<\delta c(z_i)>]\ast [\delta c_l(z_j)-<\delta c(z_j)>],$$

$${\sigma }_{mn}=\sum _{k=1}^{N_f}{Q}_{mn}^{(k)}{\sigma }_{mn}^{(k)},$$

$$Q_{mn}^{(k)}=\frac{1}{c_0}\int {\psi }_m^{\ast }{\theta }_k(z){\psi }_{n}dz,$$

Full text

{history}

1 Introduction

Development of efficient and computationally cheap methods for modeling of long-range sound propagation in randomly inhomogeneous waveguides is one of the most important problems of ocean acoustics. Indeed, even weak spatial variability of the sound speed, caused by oceanic internal waves, becomes a substantial factor on the ranges of few hundred kilometers [1, 2, 3, 4]. Therefore, calculations with range-independent or adiabatically varying models of medium couldn’t reproduce actual properties of sound propagation. On the other hand, acoustical experiments are commonly conducted with pulse broadband signals. In this way, the most simple approach based on the synthesis of an acoustic pulse from monochromatic components, together with Monte-Carlo sampling, can be excessively cumbersome. To avoid the above difficulties, one can utilize the kinetic approach based on construction of the master equation for mode amplitudes of a wavefield [4, 5, 6]. This approach can be used for studying spatial and temporal coherences of a wavefield [7, 8].

Another promising approach is modeling of sound propagation by means of the random matrix theory (or shortly RMT) [9, 10]. Its important advantage is the ability to model interference patterns and acoustic timefronts without solving the wave equation directly. In the RMT approach, complex-valued amplitudes of inter-mode transitions are considered as statistically independent random quantities, with variances determined by the random inhomogeneity of a waveguide. In this way, one basically ignores correlations between different transitions. These correlations, however, might be responsible for various coherent phenomena that can be observed in experiments. For example, they can give rise to regularly propagating beams formed by coherent ray clusters [11]. As it was shown in Ref. [12], coherent features are reflected in spectral properties of the propagator governing wave evolution. So, one can ask: whether the RMT approach is able to reproduce actual properties of sound propagation? This is namely the issue we address in the present work. We compare results obtained with the random matrix theory with those obtained by solving the parabolic equation. We restrict our attention by long-range sound propagation in the deep ocean. In particular, we use the canonical Munk waveguide and the waveguide with the biexponential sound-speed profile [2, 11]. Both models include random inhomogeneity caused by internal waves.

The manuscript is organized as follows. A wavefield propagator is introduced in the next section. Section 3 describes construction of the propagator by means of the random matrix theory. Section 4 is devoted to models of an underwater sound channel used in numerical simulation. In Section 5 we examine the ability of the RMT approach to describe dynamics of modal amplitudes. Some fundamental features of spectral analysis of the propagator are described in Section 6. Results of numerical simulation of spectral statistics are presented in Sections 7. In Summary, we give the main conclusions and outline perspectives of the future research.

2 Unitary propagator for sound propagation in the ocean and its spectral properties

Consider an underwater sound channel where spatial variability of sound speed is represented by the sum

[TABLE]

where $z$ is depth, $r$ is the range coordinate, $c_{0}$ is a reference sound speed, $\Delta c(z)$ corresponds to a background range-independent sound-speed profile, and $\delta c(r,z)$ is random sound-speed inhomogeneity caused by oceanic internal waves. Terms entering into Eq. (1) satisfy the inequality

[TABLE]

Owing to this inequality, we can invoke the small-angle approximation, when an acoustic wavefield is governed by the standard parabolic equation

[TABLE]

where wave function $\Psi$ is related to acoustic pressure $u$ by means of the formula $u=\Psi\exp(ik_{0}r)/\sqrt{r}$. Quantity $k_{0}$ is a reference wavenumber related to sound frequency $f$ as

[TABLE]

Functions $U(z)$ and $V(r,z)$ are determined by spatial sound-speed variations. In the the small-angle approximation they can be expressed as

[TABLE]

In the present paper we use idealistic perfectly-reflecting boundary conditions of the form

[TABLE]

Solution of the parabolic equation (3) can be represented as a sum over normal modes

[TABLE]

The normal modes and the corresponding eigenvalues satisfy the Sturm-Liouville problem

[TABLE]

In a range-independent waveguide ($\delta c=0$) we have

[TABLE]

i. e. modes are uncoupled and there is no transfer between different modes. Horizontal inhomogeneity of a waveguide leads to mode coupling.

As long as the boundary conditions (6) provide preservation of the wavefunction norm, general solution of the parabolic equation at range $r=r_{\mathrm{f}}$ can be formally written in terms of an unitary propagator $\hat{G}$ [13] acting as

[TABLE]

In fact, unitarity (i. e. norm preservation) of the propagator is provided by somewhat artificial assumptions. Indeed, the norm preservation means that we neglect sound absorption within the ocean bottom and the water volume. In the case of the deep ocean, interaction with the bottom doesn’t affect a large group of the lowest modes which are responsible for long-range propagation. As about the water volume attenuation, it is fairly weak for frequencies below 100 Hz and doesn’t influence process of sound scattering. Therefore, the approximation of the unitary evolution is quite reasonable for our purposes.

The propagator $\hat{G}$ can be represented in the matrix form using the basis of normal modes (8). The corresponding matrix elements are

[TABLE]

where $\hat{G}(r_{0},r_{\mathrm{f}})\psi_{n}$ is a solution of the parabolic equation at range $r=r_{\mathrm{f}}$ for the initial condition $\Psi(r=r_{0})=\psi_{n}$. Since $G_{mn}=G_{nm}^{*}$ for $m\neq n$, the matrix of the propagator is Hermitian.

3 Construction of the propagator by means of the random matrix theory

In the present section we give brief description of the random matrix approach proposed in Refs. [9, 10]. As a starting point, we express the propagator $G(r_{0},r_{\mathrm{f}})$ as a product of propagators for intermediate segments of a waveguide:

[TABLE]

Hereafter we set $r_{0}=0$ and consider the equally spaced partition of the interval $[0:r_{\mathrm{f}}]$. The spacing $\Delta r=r_{j}-r_{j-1}$ can be chosen large enough to provide statistical independence of propagators for neighboring segments. It is reasonable to assume that $\delta c(r)$ is a stationary stochastic process. Then we can replace $\hat{G}(r_{j-1},r_{j})$ by $\hat{G}(\Delta r)$ and rewrite Eq. (12) as

[TABLE]

where the index $j$ enumerates propagators corresponding to different realizations of $\delta c$.

To construct a realization of a single-segment propagator, it is reasonable to utilize the matrix representation (11). Then we can replace the propagators $\hat{G}$ in Eqs. (12) and (13) by the corresponding matrices $\mathbf{G}$. The first-order perturbation theory yields

[TABLE]

where $\mathbf{I}$ is the identity matrix, and $\mathbf{\Lambda}$ is a diagonal matrix with elements

[TABLE]

where $\delta_{mn}$ is the Kronecker symbol. In Eq. 14, $\mathbf{A}$ is an inhomogeneity-induced perturbation matrix whose elements are calculated as

[TABLE]

[TABLE]

The key idea of the random matrix approach is to treat matrix elements of the perturbation $\mathbf{A}$ as random quantities

[TABLE]

where $\sigma_{mn}$ is calculated from spectral properties of the random inhomogeneity, and $z_{mn}$ is a complex-valued Gaussian random variable with the unit variance.

The first-order perturbation theory basically doesn’t maintain unitarity. To restore unitarity without altering significantly physical properties, authors of Refs. [9, 10] offer to use the Cayley transform

[TABLE]

In this case the matrix $\mathbf{G}$ is guaranteed to be unitary.

Mode amplitudes of a wavefield can be combined into the vector $\vec{a}$, $\vec{a}\equiv(a_{1},a_{2},\text{...},a_{M})^{T}$. In accordance with Eq. (10), range evolution of this vector is governed by the equation

[TABLE]

So, a wavefield can be calculated by means of sequential multiplication of the vector of mode amplitudes by the propagator matrix. This algorithm is extremely fast if the matrix size (i. e. number of trapped modes) is not very large. Indeed, complexity of multiplication of a matrix $\mathbf{G}$ by the vector $\vec{a}$ scales as $M^{2}$, where $M$ is number of trapped modes.

4 Models of an underwater sound channel

In this work we consider two models of an underwater sound channel. The first one is the so-called biexponential (or shortly BEP) model originally introduced in Ref. [21] and further developed in Refs. [2, 11]. In the BEP model, the background sound-speed profile is determined by equation

[TABLE]

Here $c_{0}=1490$ m/s, $a=0.5$ km*-1*, $b=0.557$, and $\kappa=0.6065$. Ocean bottom is located at the depth of 4 km. An important advantage of the BEP model is the presence of analytical expressions for some important waveguide characteristics, like ray cycle length, ray travel times, e. t. c. [2, 11]

The second considered model is the celebrated canonical Munk waveguide with the background sound-speed profile

[TABLE]

where $c_{0}=1490$ m/s, $B=1$ km, $\gamma=0.0057$, and the ocean bottom is located at $z=5$ km. Both sound-speed profiles are shown in Fig. 1.

Sound-speed perturbation induced by internal waves is constructed according to the scheme presented in Refs. [10, 22]. In that scheme, one assumes linear dependence of sound-speed variations $\delta c$ on vertical displacement of a fluid parcel $\zeta$

[TABLE]

where $g=9.8$ ms*-2* is the gravitational acceleration, and $N$ is Väisäla-Brunt frequency. Wavefield of $\zeta$ is given by the double sum

[TABLE]

where $\xi(z)=e^{-z/B}-e^{-h/B}$, and $\phi_{jl}$ are random phases uniformly distributed over the interval $[0:2\pi]$. Spectral weights are described by formula

[TABLE]

where

[TABLE]

The following values of parameters are taken: $N_{0}=2\pi/10$ min, $f_{\text{i}}=1$ cycle per day, the Garrett-Munk energy $E=6.3*10^{-5}$, mode scaling number $M=(\pi j_{*}-1)/2j_{*}^{2}$, and the principle mode number $j_{*}=3$. We take 1000 values of the horizontal internal wave number $k_{l}$, which are equally spaced within the interval from $k_{l}=2\pi/100$ to $2\pi$ radians per km, with spacing $\Delta k_{l}$. A single realization of $\zeta(r,z)$ is plotted in Fig. 2.

Construction of a proper random matrix ensemble for the propagator $\hat{G}$ can be further facilitated by expanding the inhomogeneity over empirical orthogonal functions [23, 24]

[TABLE]

This expansion is also known as the Karhunen-Loeve expansion. Empirical orthogonal functions $\theta_{n}(z)$ are the eigenvectors of the covariance matrix $\mathbf{K}$ whose elements are given by

[TABLE]

where the summation goes over $l_{\text{max}}$ statistically independent realizations of и $\delta c(z)$, $\{z_{i}\}$ is a vector of depth values, and angular brackets denote statistical averaging. Taking into account that $\delta c$ is an oscillating function, we can set $<\delta c_{\mathrm{iw}}>=0$. Eigenvalues of $\mathbf{K}$ quantify contributions of the corresponding eigenvectors in the expansion (27). We retain 100 empirical eigenfunctions with the largest eigenvalues that provides almost complete reproduction of an internal wave field. Then variances in Eq. (18) are calculated as

[TABLE]

where

[TABLE]

and $\sigma_{mn}^{(k)}$ is variance of the integral

[TABLE]

Utilizing the Karhunen-Loeve expansion, we separate integration over $r$ and $z$ from each other, and thereby reduce computational efforts to calculate matrix elements $A_{mn}$.

5 RMT modeling vs direct calculation. I. Mode amplitudes

Now let’s consider how efficiently the RMT-based propagator reproduces evolution of the modal spectrum of a wavefield. We compare results obtained via the RMT with those obtained by means of direct Monte-Carlo sampling. The direct calculation was performed using the Crank-Nicolson scheme for numerical solution of the parabolic equation. The RMT modeling was carried out with two values of the propagator step $\Delta r$, 10 and 50 km. The case of $\Delta r=10$ km corresponds to smaller magnitudes of matrix elements $A_{mn}$, thus providing better accuracy of the first-order perturbation theory the RMT approach relies upon. However, it can inaccurately approximate of the actual propagator in the presence of cross-mode correlations whose length exceeds 10 km. From this viewpoint, the case of $\Delta r=50$ km looks as a more reasonable choice.

We consider two values of signal frequency $f$, 25 and 75 Hz. The case of $f=25$ Hz anticipates nearly-adiabatic propagation of modes, with relatively weak inter-mode coupling, while in the case of $f=75$ Hz ray-like effects are expected to be relevant. Ensemble of 1000 realizations was used for statistical averaging.

Figure 3 demonstrates range dependences of the ensemble-averaged first-mode energy for various waveguide models and frequencies. The calculations were performed using Eq. (20) for the initial state corresponding to the point source located at the waveguide axis, i. e. at the depth of 1000 meters. For all cases, the first-mode energy decays with range that reveals tendency towards the equipartition regime [5]. In the case of 25 Hz the rate of decay is lower indicating weaker scattering, than in the case of 75 Hz. Notably, direct Monte-Carlo simulations expose a little bit faster decay as compared with modelling via the random matrix theory. Thus, one can conclude that the random matrix theory slightly underestimates scattering. However, distributions of mode energies over modes show that the differences are mainly concerned with the lowest modes, while higher ones have nearly the same weights for all methods of simulation (see Fig. 4).

6 Propagator spectrum

Given a realization of random inhomogeneity, properties of sound propagation are reflected in the spectrum of the propagator $\hat{G}$. The corresponding eigenvalues and eigenfunctions obey the equation

[TABLE]

Taking into account the mathematical equivalence between the parabolic equation (3) and the Schrödinger equation, we can study spectral properties of the propagator by means of the methods developed in quantum mechanics. In particular, we can utilize the spectral theory of quantum chaos that is mainly based on the theory of random matrices.

The eigenvalues can be found as eigenvalues of the propagator matrix $\mathbf{G}$. Owing to the unitarity, they can be recast as

[TABLE]

Since eigenvalues of the propagator belong to the unit circle in the complex plane, its matrix corresponds to the circular ensemble of random matrices [14].

It is important to note the equivalence between spectrum of $\hat{G}(0,r_{\mathrm{f}})$ and spectrum of the Floquet operator governing wave evolution in a fictitious range-periodic waveguide with sound-speed inhomogeneity described as

[TABLE]

where $\delta c(r,z)$ corresponds to some individual realization of the actual (i. e. not range-periodic) waveguide. This property allows for interpretation of propagator spectrum in terms of the well-developed Floquet theory and thereby significantly facilitates the analysis of coherence in wave dynamics. Floquet theory was utilized for studying wave chaos in sound propagation in the ocean in Refs. [15, 16, 17]. It has to be kept in mind that this equivalence holds only if we restrict ourselves by the range interval $r\in[0:r_{\mathrm{f}}]$. However, the final range $r_{\mathrm{f}}$ of the propagator can be chosen arbitrarily. So, if we want to study long-range wave evolution, then we have to set a corresponding value of $r_{\mathrm{f}}$.

As it was argued in Refs. [3, 12], wave chaos associated with scattering on internal waves reveals itself in the statistics of level spacings. A level spacing for the circular ensemble is defined as

[TABLE]

where the sequence of eigenphases $\varphi_{m}$ is rearranged in the ascending order, $M$ is the total number of eigenvalues for a single realization of the propagator, equal to the number of trapped modes. Statistical distribution of level spacings is connected to all $m$-order correlation functions [14]

[TABLE]

where $P(\varphi_{1},\text{...},\varphi_{m})$ is the joint eigenphase probability distribution. Hence, level spacing statistics serves as a good indicator of differences between the spectrum of the propagator constructed via the RMT and the actual propagator obtained by solving the parabolic equation. Furthermore, level spacing distribution qualitatively depends on the strength of mode coupling. For instance, if scattering on inhomogeneity is weak and modes don’t significantly interact, the corresponding eigenphases of the propagator are statistically independent from each other, and level spacing distribution obeys the Poisson law

[TABLE]

In the opposite case of strong scattering and global cross-mode coupling, the neighboring eigenphases “repulse” that leads to level spacing statistics described by the Wigner surmise

[TABLE]

where constants $\alpha$ and $C$ depend on symmetries of the propagator. $\rho(s)$ is normalized as

[TABLE]

As the unitarity is the only constraint on the propagator, the propagator corresponds to the circular unitary ensemble (CUE), and we can set $\alpha=2$ and $C=4/\pi$ [18].

If the matrix of the propagator can be divided into two blocks, one corresponding to modes experiencing strong scattering, and near-diagonal one corresponding to modes maintaining coherence, then the resulting level spacing distribution has some intermediate form between the Poisson and Wigner statistics. In the high-frequency case one can use the Berry-Robnik distribution [19]

[TABLE]

where $v_{\mathrm{c}}$ and $v_{\mathrm{r}}$ denote fractions of strongly and weakly scattered modes, respectively, $v_{\mathrm{r}}+v_{\mathrm{c}}=1$. If there is good agreement between ray and wave descriptions, $v_{\mathrm{r}}$ is equal to relative volume of regular dynamics in phase space of ray equations, while $v_{\mathrm{c}}$ is equal to relative phase space volume attributed to chaotic ray dynamics. Berry-Robnik distribution undergoes the smooth transition from the Poisson to the Wigner law as $v_{\mathrm{r}}$ decreases from 1 to 0. Strictly speaking, the case of $v_{\mathrm{r}}=0$ corresponds to the Wigner distribution with $\alpha=1$, that is not identical but roughly similar to the case of $\alpha=2$. Thus, we can track scattering-induced decoherence of a wavefield by fitting level spacing distribution using the Berry-Robnik formula (39). It is worth reminding that long-range sound propagation is commonly realized with low-frequency sound, when one shouldn’t expect excellent correspondence between wave motion and the underlying ray dynamics. Therefore, the formula (39) allows us to find only approximate fraction of weakly scattered modes.

Another way to track the transition from the Poisson to Wigner statistics is the approximation of level spacing distribution using the Brody formula

[TABLE]

where

[TABLE]

$\Gamma(\mathrm{...})$ is the Euler gamma function. Eq. (40) undergoes the transition from the Poisson to the Wigner law with increasing of the parameter $\beta$ from 0 to 1. Unfortunately, the Brody distribution is empirical, and the parameter $\beta$ doesn’t have certain physical interpretation. However, formula (40) often provides more accurate fit of level spacing statistics than the Berry-Robnik and other formulae [14, 20]. So, if we are rather interested in comparison of level spacing statistics obtained using different methods, the Brody distribution is a proper choice for the fitting function.

Analysis of propagator eigenfunctions can be carried out without any auxiliary assumptions and approximations, in contrast to analysis of eigenvalues. Therefore, it can be considered as a more robust tool. Eigenfunctions of the propagator can be expanded over normal modes of the background waveguide

[TABLE]

where $C_{mn}$ is $m$-th element of the $n$-th eigenvector of the propagator matrix. Introduce a quantity

[TABLE]

that can serve as an identifier of a propagator eigenfunction [15]. In an unperturbed waveguide, only one normal mode contributes to each eigenfunction, and $\mu$ coincides with the number of this mode.

According to the definition (32), any eigenfunction of $\hat{G}$ remain invariant after propagation from $r=0$ to $r=r_{\mathrm{f}}$. Its form qualitatively depends on scattering on inhomogeneity. As scattering determines rate of mode coupling, an eigenfunction corresponding to propagation with strong impact of scattering is compound of many normal modes of the background waveguide. Hence, we can quantify strength of scattering by means of the participation ratio in the expansion (42). Participation ratio of the $n$-th eigenfunction is calculated as

[TABLE]

$\nu$ is equal to 1 in a range-independent waveguide, and increases as scattering intensifies.

Statistical analysis of the propagator’s spectrum for a model of the underwater sound channel in the Sea of Japan was performed in Ref. [12]. It was shown that spectral statistics qualitatively changes with increasing $r_{\mathrm{f}}$ indicating onset of global incoherent mode coupling, excepting for a small group of low-number modes associated with a weakly divergent beam.

7 RMT modeling vs direct calculation. II. Spectral statistics

The present section represents comparative analysis of the propagator spectral characteristics obtained via the direct calculation and those obtained by means of the random matrix theory. Level spacing distributions were averaged in the following way:

[TABLE]

where $\rho_{n}(s)$ is a level spacing distribution corresponding to the $n$-th realization of the propagator, $N=1000$ is number of realizations.

7.1 Signal frequency 25 Hz

Figure 5 demonstrates level spacing distributions for various ranges. Markedly, there is discrepancy between the results of direct simulation and RMT modeling for short ranges. In particular, direct simulation reveals peaks for relatively large spacing values (approximately 1.7 and 3.0) that are absent in the curves obtained with RMT. The discrepancy ceases with increasing range, and finally, for $r_{\mathrm{f}}=1000$ km, the distributions almost completely coincide. The distributions for $r_{\mathrm{f}}=1000$ km reveal substantial eigenvalue repulsion, as in the case of the Wigner surmise (37). However, fitting of the distributions by means of the Brody formula (40) results in strongly fluctuating range dependence of the Brody parameter $\beta$ (see Fig. 6). Almost the same behavior is also observed for the fitting using the Berry-Robnik formula (not shown). Presence of strong fluctuations can indicate inconsistence of level spacing statistics with the aforementioned analytical expressions.

Eigenfunction distributions in the plane of parameters $\mu$ and $\nu$ allow for visualization and physical interpretation of coherent phenomena [12, 25]. They are depicted in Fig. 7 for the BEP model and Fig. 8 for the Munk model. Each plot in Figs. 7 and 8 is superposition of eigenfunction distributions corresponding to individual realizations. Coherent features are seen as some structures with much higher density of plotting points as compared to fuzzy background. For instance, one can see that resulting distributions include vertical stripes, well-resolved or slightly smeared. Each such stripe corresponds to a set of eigenfunctions having nearly the same value of the identifier $\mu$. There are two kinds of the stripes, differing by interval of $\nu$ they belong to. The stripes adjoining to the minimally accessible value $\nu=1$ are formed by eigenfunctions that almost coincide with normal modes of the background waveguide. They are especially pronounced in the plots corresponding to $r_{\mathrm{f}}=50$ km.

There are also long vertical slightly smeared stripes that don’t adjoin $\nu=1$. They belong the interval $2\leq\nu\leq 6$ for the BEC model, and the interval $2\leq\nu\leq 8$ for the Munk waveguide. These stripes correspond to the phenomenon that was firstly described in Ref. [26] and referred to as mode-medium resonance. Mode-medium resonance means coherent dynamics of some compact group of modes. Onset of such group is associated with the principle of ray-mode duality, i. e. projection of underlying ray motion onto wave dynamics in the modal space. Particularly, mode-medium resonance for the propagator $\hat{G}(r_{\mathrm{f}})$ is a wave manifestation of the ray-medium resonance that arises in the one-step Poincaré map [12, 27] being the ray counterpart of the propagator $\hat{G}$. Ray-medium resonance corresponds to the condition

[TABLE]

where $l_{1}$ and $l_{2}$ are integers, $D$ is ray cycle length (ray double loop), $I$ is ray action determined by

[TABLE]

and $\chi(z)$ is ray grazing angle. The integration in Eq. (47) goes over one cycle of a ray trajectory in a background waveguide with $\delta c=0$. As long as the final range $r_{\mathrm{f}}$ can be chosen arbitrarily, we can call resonance (46) as quasi-resonance.

In the case of the BEP model, we can write down exact expressions for $D$ and $I$:

[TABLE]

[TABLE]

where

[TABLE]

is invariant in the range-independent waveguide,

[TABLE]

Ray action (47) allows one to establish ray-mode correspondence by means of the Einstein-Brillouin-Keller (EBK) quantization rule that is

[TABLE]

for purely-water rays, and

[TABLE]

for rays reflecting from the ocean surface. Replacing $m$ by $\mu$ in Eqs. (52) and (53), and using Eqs. (48)-(51), one can identify quasi-resonances responsible for vertical stripes for the BEP model. In the case of the Munk waveguide, dependence $D(I)$ can be found numerically [28]. So, one can find out that the plots corresponding to $r_{\mathrm{f}}=50$ km depict mode-medium quasi-resonance with $l_{1}=l_{2}=1$. Modes corresponding to this quasi-resonance obey to Eq. (52) in the case of the Munk model, and to Eq. (53) in the case of the BEP model.

Rays belonging to neighbourhood of the resonant action $I_{\text{res}}$ corresponding to some pair ($l_{1}$,$l_{2}$) form a cluster that maintains coherence until this quasi-resonance is well-isolated in the action space from quasi-resonances with other values of $l_{1}$, $l_{2}$ and $I_{\text{res}}$. Eqs. (52) and (53) allow one to project the cluster onto normal modes. Thus, the relating normal modes obey inequality

[TABLE]

where $m_{\text{res}}\simeq k_{0}I_{\text{res}}$ is mode number corresponding to the exact ray-medium quasi-resonance (46), and $\Delta m$ is related to half-width of ray-medium quasi-resonance $\Delta I$ as $\Delta m=k_{0}\Delta I$. Rays compounding a resonance-induced cluster have close travel times, providing spontaneous wavefield focusing in the time domain [11, 27]. This circumstance makes them observable in experiments. As it was shown in Ref. [27], density of ray-medium quasi-resonances linearly increases with increasing $r_{\mathrm{f}}$, therefore, all such clusters, as well as related modes, eventually overlap and loss their coherence.

It is worthwhile to mention that spectrum of the RMT-based propagator with $\Delta r=10$ km for the BEP model doesn’t have resonance-induced vertical stripe, in contrast to the case of $\Delta r=50$ km. In the Fig. 8(b) corresponding to the Munk waveguide, the stripe is visible but much shorter than in Figs. 8(a) and 8(c). Indeed, rays corresponding to the quasi-resonance have cycle length $D$ of 50 km. It is much larger than $\Delta r=10$ km. So, it turns out that the propagator with $\Delta r=10$ km ignores cross-mode correlations giving rise to mode-medium quasi-resonance. It leads to the conclusion that the correct reproduction of coherent resonance-induced phenomena requires $\Delta r$ to exceed the maximal value of cycle length for rays propagating without reaching the bottom. This rule, however, doesn’t apply if mode-medium quasi-resonance losses its stability already for $\Delta r<D$. It may take place for relatively high frequencies, in the presence of scattering on fine-scale inhomogeneities [3].

Another noticeable feature of eigenfunction distributions shown in Figs. 7 and 8 is the presence of well-resolved curved lines having form of bridges connecting points $(\mu=\mu^{\prime},\nu=1)$ and $(\mu=\mu^{\prime\prime},\nu=1)$. As it was shown in Ref. [12], these “bridges” correspond to resonant inter-mode transitions obeying

[TABLE]

Complete well-drawn “bridges” appear if resonance (54) is well-isolated, and the corresponding eigenfunction of the propagator can be fairly represented as sum

[TABLE]

Then normalization condition for eigenfunctions

[TABLE]

yields

[TABLE]

Rewriting $C_{mj}$ and $C_{nj}$ as

[TABLE]

we obtain

[TABLE]

Basically $\vartheta$ is a random quantity and depends on realization of the sound-speed inhomogeneity. Superposition of large number of propagator realizations uniformly fills full range of $\vartheta$ values, from $-\pi$ to $\pi$. It results in drawing the parametric curve $\mu(\vartheta),\nu(\vartheta)$ as a bridge-like pattern. If resonance (54) is poorly isolated, the bridge-like pattern becomes distorted or even completely dissolved. The latter eventually happens with increasing of $r_{\mathrm{f}}$ due to multiplication of triplets $(l,m,n)$ satisfying Eq. (54).

In general, appearance of resonance-assisted features is linked to the aforementioned spectral equivalence of the propagator $\hat{G}$ and the Floquet operator for the fictitious waveguide with sound speed inhomogeneity given by Eq. (34). As long as positions of quasi-resonances in the modal space depend on $r_{\mathrm{f}}$, we conclude that the resonance-induced structures don’t mean enhanced transfer between certain modes. Their actual meaning is indication of ability of related modal groups to maintain coherence for $r\leq r_{\mathrm{f}}$. Presence of such structures is an unambiguous signature of long-living cross-mode correlations. It is worth emphasizing that such structures are also observed in the plots obtained using the RMT approach. So, it turns out RMT is able to reproduce quite subtle coherent effects. This allows one to classify the RMT approach as a really powerful method of modeling.

As $r_{\mathrm{f}}$ increases number of quasi-resonances (46) and (54) grows. Since that isolation of individual quasi-resonances is violated. Overlapping of quasi-resonances results in their destruction [27]. So, the quasi-resonances are no more able to produce coherent structures in the $\mu$-$\nu$ plane. The resulting pattern acquires the form of a fuzzy spread boomerang, like in Figs. 7(g-i) and 8(g-i).

7.2 Signal frequency 75 Hz

Increasing of signal frequency substantially imrpoves correspondence between the RMT-based propagator and its directly-calculated counterpart. In the case of $f=75$ Hz the corresponding level spacing distributions almost fully coincide, as it is shown in Fig. 9. They all expose transition to Wigner statistics as $r_{\mathrm{f}}$ grows. It is well demonstrated by range dependence of the Brody parameter corresponding to the best fit (see Fig. 10). It also indicates good agreement between the RMT and Crank-Nicolson data. Notably, the curve corresponding to RMT data with $\Delta r=50$ km is closer to the results of direct solving that one corresponding to $\Delta r=10$ km. Fitting level spacing distribution by means of the Berry-Robnik formula (39), we can estimate fraction of modes which don’t experience strong scattering and propagate coherently. Figure 11 shows that this fraction rapidly decreases with range and completely vanishes for $r_{\mathrm{f}}$ of nearly 400 km.

Eigenfunction distributions in the $\mu$-$\nu$ plane exhibit signatures of mode-medium quasi-resonance only for relatively short ranges. RMT modeling is able to reproduce it, but with much less resolution than in the case of the direct solution. For longer ranges, the distributions rapidly achieve the boomerang-like pattern representing a kind of statistical equillibrium.

So, we can see that coherent structures are remarkably suppressed for $f=75$ Hz, as compared to the case of $f=25$ Hz. These drastic differences in spectral statistics indicate qualitatively different mechanismes of scattering, namely transition to ray-like scattering with impact of ray chaos.

8 Summary

In the present paper we examine the approach based on construction of a wavefield propagator using the random matrix theory (RMT). In particular, we study its ability to reproduce properties of the propagator constructed by means of direct numerical solution of the parabolic equation. It is found out that noticeable differences between the RMT-based propagator occur only for relatively short ranges, below 100 km, where cross-mode correlations play an essential role and give various coherent impacts to spectral statistics. The coherent impacts are caused by quasi-resonances having a deterministic origin and especially pronounced for frequencies of few tens Hz. To suppress influence of cross-mode correlations, the range step $\Delta r$ of a RMT-based propagator should exceed the maximal cycle length of rays propagating without reflections off the bottom. For longer ranges, the random matrix theory provides excellent agreement with direct solutions of the parabolic equation. Calculations of mode energies show that the RMT-based method slightly underestimates mode coupling, but the differences with the direct solutions are observed only for the lowest modes.

The main physical result of the present work is the evident influence of ray-based features onto propagator spectra for low signal frequencies. Somewhat unexpectedly, they are much more revealed in the propagator spectrum at the frequency of 25 Hz than at the frequency of 75 Hz. Mode-medium quasi-resonances, being modal counterparts of ray-medium quasi-resonances, are responsible for specific vertical stripes in distributions of propagator eigenfucntions in the plane with parameters $\mu$ and $\nu$ as coordinates. In the ray picture, the corresponding ray-medium quasi-resonances are mainly destroyed due to ray chaos. Appearance of resonance-assisted vertical stripes in the eigenfunction statistics for low frequencies indicates that stability of the corresponding mode-medium quasi-resonances is restored. It is closely related to the recovery of ray stability reported in Ref. [16, 29]. One of such mechanismes is ray scattering on fine-scale distortions of a sound-speed profile [30], which are also present in the models considered in the present paper. For low frequencies such distortions become irrelevant for wave dynamics. Therefore, refraction of a wavefield becomes meaningfully stabilized, as compared to the high-frequency regime. Stabilization of refraction is reflected in dynamics of modal amplitudes [31] and results in recovery of mode-medium resonance.

In the present paper we consider only relatively low frequencies. In the case of higher frequencies we can expect qualitatively different properties of sound scattering, for example, more significant influence of fine-scale structures in the sound-speed profile. It can be reflected in spectral statistics of a wavefield propagator. Another important issue is applicability of the RMT-based approach in the case of waveguides with a slowly-varying background sound-speed profile. Also, it is reasonable to extend applicability of the RMT approach onto shallow-sea propagation, when the unitarity of the propagator doesn’t hold. All these issues will be addressed in forthcoming works.

Acknowledgments

This work was supported by the Russian Foundation of Basic Research within the projects 16-35-60040 and 16-05-01074, and by the POI FEBRAS Program ’Mathematical simulation and analysis of dynamical processes in the ocean’ (№117030110034-7) Author is grateful to Steven Tomsovic for stimulating and fruitful discussions, and to Leonid Kon’kov for the assistance in preparation of figures.