# Evaluation of Augmentation methods in classifying Autism Spectrum Disorders
from fMRI data with 3D Convolutional Neural Networks

###### Abstract

Classifying subjects as healthy or diseased using neuroimaging data has gained
a lot of attention during the last 10 years. Here we apply deep learning to
derivatives from resting state fMRI data, and investigate how different 3D
augmentation techniques affect the test accuracy. Specifically, we use resting
state derivatives from 1,112 subjects in ABIDE preprocessed to train a 3D
convolutional neural network (CNN) to perform the classification. Our results
show that augmentation only provide minor improvements to the test accuracy.

Index Terms—  Deep learning, autism, functional MRI, convolutional neural
network, augmentation

##  1 Introduction

Ever since the emergence of magnetic resonance imaging (MRI) in the 1980’s,
the absence of ionizing radiation and the flexibility of the acquisition
procedure have made this an increasingly important imaging modality in the
clinical sciences. The lack of contrast between different tissues in the brain
and the interference of the mineralized tissue around it when using x-ray
techniques, make MRI especially useful in neuroimaging.

While a wide variety of neurological conditions can be diagnosed with MRI,
psychiatric anomalies, however, have proven illusive to detect. Presumably
this is because these affect many systems distributed throughout the brain and
their manifestations are likely subtle as well as time variant. Functional MRI
(fMRI) is a technique that seems particularly suited to capture this
information and several large collaborative efforts have been made to collect
and share (resting state) fMRI data, for different psychiatric conditions [1].

ABIDE [2] is one such effort that make available data for 539 subjects with
Autism Spectrum Disorders (ASD) as well as 573 typical controls. We are using
machine learning in an endeavour to classify fMRI data according to the
presence or absence of ASD. The problem seems hard in that accuracies seldom
rise to more than 70% when the model classifies unseen data [3, 4, 5]. While
1,112 subjects is a very large fMRI dataset, it is still small from a deep
learning perspective (for example, the popular ImageNet database contains
several million images). To further increase the size of the training dataset,
and to make CNNs robust to transformations such as rotation, data augmentation
is often used [6]. In our recent work we investigated what the best 3D
augmentation is for brain tumor segmentation [7]. In this paper we instead
want to see if 3D augmentation can help train a better classifier as well as
what kind of augmentation techniques work the best. Thomas et al. [4] recently
used deep learning for the same ABIDE dataset, but do not mention anything
about augmentation.

We have therefore implemented a 3D convolutional neural network (CNN) of
modest size and trained it with resting state fMRI deriviatives from ABIDE,
using no augmentation as well as using different augmentation strategies.
Briefly, our results show some improvements with carefully chosen augmentation
strategies but the effects are in the order of magnitude of a percentage
point. Out of the commonly used augmentation techniques that we evaluated,
many proved to be deleterious in this context.

##  2 Methods

###  2.1 Data

ABIDE preprocessed [8]111http://preprocessed-connectomes-project.org/abide/
shares preprocessed data from resting state fMRI in various forms. As all the
preprocessing has been done, we can focus on the machine learning part, and
other researchers can use the same preprocessed data to reproduce our
findings. We downloaded 3D data of derivatives where the time dimension had
been collapsed into different forms of statistics.

Specifically, we downloaded fMRI data (derivatives) processed with the
Connectome Computation System (CCS) pipeline [9] which performs slice timing
correction, motion realignment and global intensity normalisation. The data
was cleaned from confounders by regressing with the head movement parameters,
the time dependent mean intensity as well as regressors for linear and
quadratic drift. Each time series was also band pass filtered (0.01 - 0.1 Hz).
This preprocessing corresponds to the strategy called global_filt. The data
was furthermore registered to the MNI152 brain template using boundary based
rigid body registration [10] for functional to anatomical registration, and
FLIRT and FNIRT for anatomical to template registration [11].

After preliminary testing of the 10 calculated derivatives available in ABIDE
preprocessed, we chose to use Regional Homogeneity (ReHo) derivative for
comparing different augmentation strategies. Each derivative volume from the
resting state fMRI data has a size of 61 x 73 x 61 voxels, which is fed into
the 3D CNN described below. The 539 subjects with ASD and the 573 controls
were split 70/15/15 into training, validation and test sets.

###  2.2 Deep learning

The 3D CNN was implemented using Keras and consists of 3 convolutional layers
(ReLU activation), max-pooling layers, a dense layer with 16 nodes, and a
final one-node layer with sigmoid activation. The first and second
convolutional layers contain 8 filters each (size 3 x 3 x 3), and the last
convolutional layer uses 16 filters. The total number of trainable parameters
in the 3D CNN is approximately $450$k. The CNN was trained with the Adam
optimizer with a learning rate of $10^{-5}$. To prevent overfitting, early
stopping was used with a patience of 50 epochs. The training was run until
validation accuracy did not improve, and the model then restored to the state
when the last improvement was seen. As an alternative, the models were also
trained for 150 epochs with no conditional stopping. To obtain more robust
estimates of the test accuracy, 10-fold cross validation was used and the mean
test accuracy was calculated.

The average training time for a single fold were between five minutes and 2.5
hours - depending on the type of on-the-fly augmentation employed - with one
Nvidia Tesla V100 graphics card for the early stopping models. For the
training with a fixed number of epochs, the average single fold training time
was at least 10 minutes but otherwise in the previously mentioned span. In the
longer training runs it is unlikely that the computation speed was bounded by
the speed of the graphics card, as the on-the-fly augmentations were done on
the CPU. In total we trained 300 3D CNNs in order to compare all settings.

###  2.3 Augmentation

There are many types of augmentation that can be useful in 3D. Rotation,
flipping and scaling are common for training 2D CNNs, and can easily also be
applied in 3D. Elastic deformations are common when training segmentation
networks, but perhaps not as common for classification. Brightness
augmentation can for example help if the data have been collected at several
different MR scanners, as they normally generate data with different
brightness [7]. While 2D augmentation is included in many deep learning
frameworks, the support for 3D augmentation is normally lacking, which forced
us to make our own implementation. The 3D augmentation techniques tested in
this study are:

  * •

Flipping: flipping of the x-axis or not, which seems to work better compared
to flipping all axes.

  * •

Rotation: rotation applied to each axis with angles randomly chosen from a
uniform distribution with range between -7.5 and 7.5 degrees, -15 and 15
degrees, -30 and 30 degrees or -45 and 45 degrees.

  * •

Scale: scaling applied to each axis by a factor randomly chosen from a uniform
distribution with range $\pm 10$% or $\pm 20$%.

  * •

Brightness: power-law $\gamma$ intensity transformation with its parameters
gain ($g$) and $\gamma$ chosen randomly between 0.8 - 1.2 from a uniform
distribution. The intensity ($I$) is randomly changed according to the
formula: $I_{new}=g\cdot I^{\gamma}$.

  * •

Elastic deformation: elastic deformation with square deformation grid with
displacements sampled from from a normal distribution with standard deviation
$\sigma=2$, 4, 6, or 8 voxels [12], where the smoothing is done by a spline
filter with order 3 in each dimension.

To investigate the effect of combining different types of augmentation, we
also trained the CNN with the two best performing augmentation approaches.

##  3 Results

The results from all the different augmentation techniques, as well as
baseline results obtained without augmentation, are presented in Figures 1
(early stopping) and 2 (fixed number of training epochs). In general the
augmentation does not have a large effect on the test accuracy. For early
stopping, random scaling seems to be the best single augmentation approach,
but the mean improvement over 10 cross-validation folds is only about 0.5
percentage units. Small elastic deformations also have a small positive
effect, while large deformations give the worst results. For a fixed number of
training epochs, elastic deformations seem to work best, with an improvement
of 2.2 percentage units.

##  4 Discussion

Compared to our previous work on augmentation for brain tumor segmentation
[7], where several augmentation techniques were shown to significantly improve
the segmentation accuracy, we only find minor improvements of the test
accuracy in this study (even though the training accuracy is well above 90%,
indicating over fitting). Volume classification is in general a harder problem
than volume segmentation, as each subject only represents a single training
example, and may partly explain our results.

Fig. 1: Test accuracy for classifying subjects as healthy or diseased for the
ABIDE dataset, for different data augmentation approaches. The error bar
represents the standard deviation over the 10 cross-validation folds. Note
that half of the augmentation approaches result in a test accuracy that is
lower compared to the baseline model trained without augmentation, but overall
the differences are small. These results were obtained when using early
stopping. Compared to no augmentation, the best augmentation approach
increases the test accuracy by 0.6 percentage units. Fig. 2: Test accuracy for
classifying subjects as healthy or diseased for the ABIDE dataset, for
different data augmentation approaches. The error bar represents the standard
deviation over the 10 cross-validation folds. These results were obtained when
using a fix number of epochs for each training. Compared to no augmentation,
the best augmentation approach increases the test accuracy by 2.2 percentage
units.

In this study brightness augmentation does not help at all, while it provided
a major improvement for brain tumor segmentation for MR images collected at
some 20 different sites [7]. A possible explanation is that the data in this
study are not raw MR images, since many preprocessing steps have been used to
normalize the intensities to a certain range, and to calculate different
derivatives. On the contrary, as the values in the derivative volumes have a
certain meaning, brightness augmentation can impair the performance.

Since all the subjects have been registered to MNI space, we hypothesized that
the results may be different if random transformations are applied to the test
volumes, but test time augmentation did not change our findings (results not
shown).

Our results are for a single preprocessing strategy (global signal regression
and bandpass filtering), and the preprocessing choice can at least in theory
affect how much the augmentation helps.

Our conclusion is that augmentation only provides minor improvements when
training 3D CNNs for classification of ASD versus controls, but the results
may be different for an easier task where the baseline test accuracy is for
example 80%. The results may also differ for other derivatives in ABIDE
preprocessed, and when using several derivatives at the same time using a
multi-channel 3D CNN.

##  5 Compliance with ethical standards

This research study was conducted retrospectively using human subject data
made available in open access by ABIDE preprocessed. Ethical approval was not
required as confirmed by the ethics committee of Linköping.

##  6 Acknowledgments

This work was supported by the Swedish research council grant 2017-04889, and
by the ITEA / VINNOVA funded project Automation, Surgery Support and Intuitive
3D visualization to optimize workflow in IGT SysTems (ASSIST). Anders Eklund
has previously received hardware from Nvidia, otherwise the authors have no
conflicts of interest to declare.

## References

  * [1] Russell A Poldrack and Krzysztof J Gorgolewski,  “Making big data open: data sharing in neuroimaging,”  Nature neuroscience, vol. 17, no. 11, pp. 1510–1517, 2014.
  * [2] Adriana Di Martino, Chao-Gan Yan, Qingyang Li, Erin Denio, Francisco X Castellanos, Kaat Alaerts, Jeffrey S Anderson, Michal Assaf, Susan Y Bookheimer, Mirella Dapretto, et al.,  “The autism brain imaging data exchange: towards a large-scale evaluation of the intrinsic brain architecture in autism,”  Molecular psychiatry, vol. 19, no. 6, pp. 659, 2014.
  * [3] Mohammad R Arbabshirani, Sergey Plis, Jing Sui, and Vince D Calhoun,  “Single subject prediction of brain disorders in neuroimaging: Promises and pitfalls,”  Neuroimage, vol. 145, pp. 137–165, 2017.
  * [4] Rajat Mani Thomas, Selene Gallo, Leonardo Cerliani, Paul Zhutovsky, Ahmed El-Gazzar, and Guido van Wingen,  “Classifying autism spectrum disorder using the temporal statistics of resting-state functional MRI data with 3D convolutional neural networks,”  Frontiers in psychiatry, vol. 11, pp. 440, 2020.
  * [5] Anibal Sólon Heinsfeld, Alexandre Rosa Franco, R Cameron Craddock, Augusto Buchweitz, and Felipe Meneguzzi,  “Identification of autism spectrum disorder using deep learning and the ABIDE dataset,”  NeuroImage: Clinical, vol. 17, pp. 16–23, 2018.
  * [6] Connor Shorten and Taghi M Khoshgoftaar,  “A survey on image data augmentation for deep learning,”  Journal of Big Data, vol. 6, no. 1, pp. 60, 2019.
  * [7] Marco Domenico Cirillo, David Abramian, and Anders Eklund,  “What is the best data augmentation approach for brain tumor segmentation using 3D U-Net?,”  in IEEE International Conference on Image Processing (ICIP), 2021\.
  * [8] Cameron Craddock, Yassine Benhajali, Carlton Chu, Francois Chouinard, Alan Evans, András Jakab, Budhachandra Singh Khundrakpam, John David Lewis, Qingyang Li, Michael Milham, Chaogan Yan, and Pierre Bellec,  “The neuro bureau preprocessing initiative: open sharing of preprocessed neuroimaging data and derivatives,”  Neuroinformatics, 2013.
  * [9] Ting Xu, Zhi Yang, Lili Jiang, Xiu-Xia Xing, and Xi-Nian Zuo,  “A connectome computation system for discovery science of brain,”  Science Bulletin, vol. 60, no. 1, pp. 86–95, 2015.
  * [10] Douglas N Greve and Bruce Fischl,  “Accurate and robust brain image alignment using boundary-based registration,”  Neuroimage, vol. 48, no. 1, pp. 63–72, 2009.
  * [11] Mark Jenkinson, Christian F Beckmann, Timothy EJ Behrens, Mark W Woolrich, and Stephen M Smith,  “FSL,”  Neuroimage, vol. 62, no. 2, pp. 782–790, 2012.
  * [12] Olaf Ronneberger, Philipp Fischer, and Thomas Brox,  “U-net: Convolutional networks for biomedical image segmentation,”  in International Conference on Medical image computing and computer-assisted intervention. Springer, 2015, pp. 234–241.


# Study of wave chaos in a randomly-inhomogeneous oceanic acoustic waveguide:
spectral analysis of the finite-range evolution operator

D.V. Makarov<EMAIL_ADDRESS>L.E. Kon’kov  M.Yu. Uleysky  V.I.Il’ichev
Pacific Oceanological Institute
of the Far-Eastern Branch of the Russian Academy of Sciences,
43 Baltiyskaya St., 690041, Vladivostok, Russia  P.S. Petrov  V.I.Il’ichev
Pacific Oceanological Institute
of the Far-Eastern Branch of the Russian Academy of Sciences,
43 Baltiyskaya St., 690041, Vladivostok, Russia  Far-East Federal University,
8 Sukhanova St., 690950, Vladivostok, Russia

###### Abstract

The proplem of sound propagation in an oceanic waveguide is considered.
Scattering on random inhomogeneity of the waveguide leads to wave chaos. Chaos
reveals itself in spectral properties of the finite-range evolution operator
(FREO). FREO describes transformation of a wavefield in course of propagation
along a finite segment of a waveguide. We study transition to chaos by
tracking variations in spectral statistics with increasing length of the
segment. Analysis of the FREO is accompanied with ray calculations using the
one-step Poincaré map which is the classical counterpart of the FREO.
Underwater sound channel in the Sea of Japan is taken for an example. Several
methods of spectral analysis are utilized. In particular, we approximate level
spacing statistics by means of the Berry-Robnik and Brody distributions,
explore the spectrum using the procedure elaborated by A. Relano with
coworkers (Relano et al, Phys. Rev. Lett., 2002; Relano, Phys. Rev. Lett.,
2008), and analyze modal expansions of the eigenfunctions. We show that the
analysis of FREO eigenfunctions is more informative than the analysis of
eigenvalue statistics. It is found that near-axial sound propagation in the
Sea of Japan preserves stability even over distances of hundreds kilometers.
This phenomenon is associated with the presence of a shearless torus in the
classical phase space. Increasing of acoustic wavelength degrades scattering,
resulting in recovery of localization near periodic orbits of the one-step
Poincaré map. Relying upon the formal analogy between wave and quantum chaos,
we suggest that the concept of FREO, supported by classical calculations via
the one-step Poincaré map, can be efficiently applied for studying chaos-
induced decoherence in quantum systems.

###### pacs:

05.45.Mt, 43.30.Cq, 03.65.Yz, 05.45.Ac, 43.30.Ft

##  I Introduction

Sound speed in the deep ocean typically has a minimum at some depth. This
results in formation of a refractive waveguide, the so-called underwater sound
channel, which prevents sound waves from contact with the absorbing bottom. As
sound absorption within water column is fairly weak, an underwater sound
channel enables sound propagation over distances of thousands kilometers. The
largest distance had been achieved using explosive charges in the seminal
experiment on sound transmission from Perth to Bermuda in 1960 Shockley  _et
al._ (1982); Munk  _et al._ (1988).

It is realized that small sound-speed variations induced by oceanic internal
waves lead to Lyapunov instability and chaos of sound rays. In the
mathematical sense, ray chaos is an analogue of classical chaos in Hamiltonian
systems. Indeed, ray motion in a waveguide is equivalent to motion of a point
particle in a potential well, and sound-speed variations along a waveguide
play the role of a nonstationary perturbation. Reciprocal Lyapunov exponent
for chaotic rays typically is of about several tens kilometers Makarov  _et
al._ (2009), therefore, the problem of ray chaos is mainly important for long-
range sound propagation. During the last two decades ray chaos in ocean
acoustics was an object of intense research, both theoretical and experimental
Smith  _et al._ (1992); Smirnov  _et al._ (2001); Brown  _et al._ (2003);
Makarov  _et al._ (2004); Beron-Vera and Brown (2004). Considerable attention
was paid to wave chaos Makarov  _et al._ (2009); Virovlyansky and Zaslavsky
(1999); Smirnov  _et al._ (2004, 2005); Kon’kov  _et al._ (2007); Makarov  _et
al._ (2008); Virovlyansky  _et al._ (2012). The term wave chaos relates to
wavefield manifestations of ray chaos. It was found that interference makes
wave refraction more regular than it is anticipated from ray modeling, albeit
influence of ray chaos persist even for very low sound frequencies Hegewisch
_et al._ (2005). This problem becomes especially important due to the growing
interest to hydroacoustical tomography, i. e. monitoring of environment using
sound signals. The classical scheme of tomography developed by Munk and Wunsch
Munk and Wunsch (1979) is based on computation of eigenrays connecting the
source and the receiver. It was shown in Tappert and Tang (1996) that ray
chaos leads to exponential proliferation of eigenrays with increasing
distance. As a result, the inverse problem becomes ill-posed, impeding
environment reconstruction. However, wave-based corrections “stabilize” wave
refraction, i. e. the standard semiclassical approximation typically
overestimates ray chaos. Thus, one needs either an improved version of the
semiclassical approximation for proper computation of eigenrays, or some
approach for making implications about eigenray stability relying upon wave
modeling, whereby a priori taking into account the wave-based suppression of
ray chaos.

Theory of ray and wave chaos extensively exploits methods borrowed from the
theory of Hamiltonian dynamical systems Makarov  _et al._ (2009); Zaslavsky
(2007). In particular, ray motion is often studied by means of the phase space
representation. This provides the clear geometric interpretation of ray
dynamics for toy models of deterministic range-periodic waveguides. For
instance, one can easily separate domains of initial conditions corresponding
to stable and chaotic rays from each other. However, realistic underwater
sound channels are not range-periodic, and their range inhomogeneity should be
rather described as a stochastic process. On the other hand, statistical
methods Dozier and Tappert (1978a); Colosi and Morozov (2009) are implicitly
based on the assumption of ergodic chaos and ignore the existence of phase
space domains of finite-range stability Wolfson and Tomsovic (2001); Makarov
_et al._ (2006). Therefore, it is reasonable to elaborate some theoretical
approach combining advantages of statistical and deterministic methods.
Recently it had been shown that a bridge between the deterministic and
statistical descriptions of wave propagation in random media can be built up
by introducing the finite-range evolution operator (hereafter FREO) Makarov
_et al._ (2010); Virovlyansky  _et al._ (2012). FREO acts as a propagator
determining wavefield transformation between two vertical sections of a
waveguide. Spectral analysis of the FREO in terms of the random matrix theory
allows one to estimate contribution of ray chaos to wave dynamics and track
transition to chaos with increasing distance between the sections. A somewhat
different approach had been offered in Hegewisch and Tomsovic (2012). These
approaches allow one to generalize the well-developed spectral theory of
quantum chaos (see, for instance, Stöckmann (2007)) on one-way wave
propagation in random media. On the other hand, concept of the FREO can serve
as a promising method for studying chaos-induced decoherence in randomly-
driven quantum systems Kolovsky (1997).

In the present paper we study spectral statistics of the FREO for the Sea of
Japan. Attention is concentrated on the track connecting the Gamov peninsula
and Kita-Yamato bank. The length of the track is about 350 km. Our interest to
this waveguide is motivated by results of the experiment conducted there in
2006 Bezotvetnykh  _et al._ (2009), indicating on high stability of near-axial
propagation. Similar behavior was observed in a earlier experiment with a
slightly different propagation track Spindel  _et al._ (2003). It should be
noted that stability of near-axial propagation is atypical for the deep ocean
Makarov  _et al._ (2009); Virovlyansky  _et al._ (2012).

The paper is organized as follows. The next section represents basic equations
describing long-range sound propagation in the ocean. In Section III, we
describe the waveguide used in the paper. Section IV contains description of
the FREO and methods of its analysis. Section V is devoted to the classical
counterpart of the FREO, namely the one-step Poincaré map. In Section VI, we
perform statistical analysis of the FREO for the underwater sound channel in
the Sea of Japan. In Conclusion we summarize and discuss the results obtained.

##  II General equations

Ocean is a layered media, and its horizontal variability is much weaker than
vertical one. This allows one to reduce the initial three-dimensional problem
of wave propagation in an underwater sound channel to a two-dimensional one by
assuming cylindrical symmetry and neglecting azimuthal coupling. Sound
refraction is governed by spatial variability of sound speed

$c(z,\,r)=c_{0}+\Delta c(z)+\delta c(z,\,r),$ (1)

where $z$ is depth, $r$ is range coordinate, $c_{0}$ is a reference sound
speed. Sound-speed variations obey the double inequality

$\lvert\delta c\rvert_{\text{max}}\ll\lvert\Delta c\rvert_{\text{max}}\ll
c_{0}.$ (2)

Left inequality implies that the range-dependent term can be treated as a weak
perturbation of the background sound-speed profile. This term is mainly
contributed from oceanic internal waves. Right inequality means that
variations of the refractive index are weak, and only those waves which
propagate with small angles with respect to the horizontal plane can avoid
contact with the absorbing ocean bottom. Thus, one can invoke the small-angle
approximation, in which an acoustic wavefield is governed by the standard
parabolic equation

$\frac{i}{k_{0}}\frac{\partial\Phi}{\partial
r}=-\frac{1}{2k_{0}^{2}}\frac{\partial^{2}\Phi}{\partial
z^{2}}+\left[U(z)+\varepsilon V(z,\,r)\right]\Phi,$ (3)

where wave function $\Phi$ is linked to acoustic pressure $u$ by means of the
formula $u=\Phi\exp(ik_{0}r)/\sqrt{r}$. Here the denominator $\sqrt{r}$
responds for the cylindrical spreading of sound. Quantity $k_{0}$ is the
reference wavenumber related to sound frequency $f$ as $k_{0}=2\pi f/c_{0}$.
Functions $U(z)$ and $V(z,r)$ are determined by spatial sound-speed
variations. In the the small-angle approximation they can be expressed as

$U(z)=\frac{\Delta c(z)}{c_{0}},\quad V(z,\,r)=\frac{\delta c(z,\,r)}{c_{0}}.$
(4)

According to (2), $\lvert V\rvert_{\text{max}}\ll\lvert U\rvert_{\text{max}}$,
that is, function $V(z,r)$ can be treated as a small perturbation. One can
easily see that the replacement

$k_{0}^{-1}\to\hbar,\quad r\to t$ (5)

transforms the parabolic equation (3) into the Schrödinger equation for a
particle with unit mass. This circumstance enables study of wave propagation
using the approaches developed in quantum mechanics. In this relationship,
function $U(z)$ serves as an unperturbed potential. As $r$ is a timelike
variable, $V(z,\,r)$ plays the role of a nonstationary perturbation.

In the short-wavelength limit $k_{0}\to\infty$, solution of the parabolic
equation (3) can be expressed as a sum of rays whose trajectories are governed
by the Hamiltonian

$H=\frac{p^{2}}{2}+U(z)+V(z,\,r),$ (6)

where $p=\tan\chi$, $\chi$ is ray grazing angle (i. e. angle with respect to
the horizontal plane). The respective Hamiltonian equations read

$\frac{dz}{dr}=\frac{\partial H}{\partial
p}=p,\quad\frac{dp}{dr}=-\frac{\partial H}{\partial z}=-\frac{\partial
U}{\partial z}-\frac{\partial V}{\partial z}.$ (7)

Owing to the analogy with classical mechanics, $p$ is referred to as ray
momentum.

##  III Model of a waveguide

###  III.1 Background sound-speed profile

Figure 1:  (a) Unperturbed sound-speed profile, (b) the first empirical
orthogonal function of the sound-speed perturbation (solid) and its smoothed
approximation (dotted).

Model of the underwater sound channel in the Sea of Japan was elaborated using
the hydrological data from the database Dat . Function $U(z)$ corresponding to
the background sound-speed profile was approximated by the expression

$U(z)=\left\\{\begin{aligned} &U_{1}(z),&z\leqslant z_{0},\\\
&U_{2}(z),&z>z_{0}.\end{aligned}\right.$ (8)

where

$\displaystyle U_{1}(z)$ $\displaystyle=\frac{c_{1}}{c_{0}}e^{-z/z_{1}},$ (9)
$\displaystyle U_{2}(z)$
$\displaystyle=\frac{c_{1}}{c_{0}}e^{-z_{0}/z_{1}}+\frac{g}{c_{0}}(z-z_{0}),$

$c_{0}=1455$ m/s, $c_{1}=70$ m/s, $z_{0}=250$ m is the depth of the channel
axis, i. e. the depth with the minimal sound speed, $z_{1}=30$ m, $g=0.017$
s-1. (see Fig. 1). The ocean bottom is assumed to be flat and placed at the
depth $h=3$ km. We consider only the deep-water propagation, albeit the source
in the aforementioned experiments Bezotvetnykh  _et al._ (2009); Spindel  _et
al._ (2003) was mounted into the bottom in the coastal zone near the Gamov
peninsula. The shallow-water part of the waveguide was relatively short, less
than 30 km, and didn’t have significant bottom features which could remarkably
alter ray stability.

Figure 2:  Ray cycle length vs initial ray momentum for the source located at
the channel axis.

Expressions (8) and (9) permit analytical derivation of the basic model
characteritics in the absence of horizontal inhomogeneity. One can introduce
the ray action

$I=\frac{1}{\pi}\int\limits_{z_{\text{min}}}^{z_{\text{max}}}\sqrt{2[E-U(z)]}\,dz,$
(10)

where $z_{\text{min}}$ and $z_{\text{max}}$ are the upper and lower ray
turning points, respectively, and $E=0.5p^{2}+U(z)$. Ray action describes
steepness of a ray trajectory and enters into the Einstein-Brillouin-Keller
quantization rule.

$k_{0}I_{m}=m-1/2,\quad m=1,2,\dots$ (11)

establishing the link between normal modes of the unperturbed waveguide and
modal rays. Here $I_{m}$ is the action of a modal ray, and both ray turning
points are assumed to be inside the water column, that is, rays undergo total
internal reflection due to smooth vertical gradient of the refractive index
$n=c_{0}/c(z,\,r)$. Normal modes are the solutions of the Sturm-Liouville
equation

$-\frac{1}{2k_{0}^{2}}\frac{\partial^{2}\phi_{n}(z)}{\partial
z^{2}}+U(z)\phi_{n}(z)=E_{n}\phi_{n}(z),$ (12)

with appropriate boundary conditions at the ocean surface and bottom. The main
contribution to the mode $m$ is given from the rays with $I\simeq I_{m}$
Virovlyansky and Zaslavsky (1999). Integration of (10) yields

$I=2z_{1}\sqrt{2E_{\text{min}}}\left[\sqrt{\varepsilon}\ln(\sqrt{\varepsilon}+\sqrt{\varepsilon-1})-\sqrt{\varepsilon-1}\vphantom{\sqrt{\mathstrut}}\right]+\\\
\frac{2\sqrt{2}c_{0}}{3g}(E-E_{\text{min}})^{3/2},$ (13)

where $\varepsilon=E/E_{\text{min}}$,
$E_{\text{min}}=(c_{1}/c_{0})\exp(-z_{0}/z_{1})$. Ray cycle length, i. e.
horizontal distance between two successive upper (or lower) ray turning
points, can be determined as

$D(E)=2\pi\frac{dI}{dE}=\\\
2z_{1}\sqrt{\frac{2}{E}}\ln\left(\sqrt{\varepsilon}+\sqrt{\varepsilon-1}\right)+\frac{2c_{0}\sqrt{2E_{\text{min}}}}{g}\sqrt{\varepsilon-1}.$
(14)

Fig. 2 represents dependence of ray cycle length on initial ray momentum for
the source located at the channel axis $z=250$ m. Initial momentum $p_{0}$
depends on $E$ as $p_{0}=\sqrt{2(E-E_{\text{min}})}$. Function $D(p_{0})$ is
nonmonotonic and has two extrema, the sharp maximum and the smooth minimum.
The latter one can give rise to a so-called weakly-divergent beam Brekhovskikh
_et al._ (1990); Smirnov  _et al._ (1999); Morozov and Colosi (2005); Petukhov
(2009). Its low divergence is associated with approximate equality of cycle
length values for rays forming the beam. It will be demonstrated in Section V
that the local minimum of $D(p_{0})$ plays a significant role in ray stabilty.

###  III.2 Sound-speed perturbation

Model of the internal-wave-induced sound-speed perturbation was built up in
several steps. Firstly, we computed the range-averaged profile of buoyancy
frequency, using the hydrological data Dat . Then, we calculated realizations
of the sound-speed perturbation using the method proposed in Colosi and Brown
(1998). In order to facilitate numerical simulation, the perturbation was
expanded over empirical orthogonal functions LeBlanc and Middleton (1980)

$\delta c(z,\,r)=\left<\delta c(z)\right>+\sum\limits_{n}b_{n}(r)Y_{n}(z).$
(15)

Empirical orthogonal functions $Y_{n}(z)$ are the eigenvectors of the
covariance matrix $\hat{K}$ with elements

$K_{ij}=\frac{1}{L}\sum\limits_{l=1}^{L}[\delta c_{l}(z_{i})-\left<{\delta
c}(z_{i})\right>][\delta c_{l}(z_{j})-\left<{\delta c}(z_{j})\right>],$ (16)

where index $l$ numbers $L$ statistically independent realizations of $\delta
c(z)$, $\\{z_{i}\\}$ is a vector of depth values, and angular brackets mean
ensemble average. As $\delta c$ is caused by internal waves, one can set
$\left<\delta c\right>=0$. Eigenvalues of the matrix $\hat{K}$ quantify
contributions from the corresponding eigenvectors in the expansion (15). It
was found that the contribution of the first orthogonal function prevails, and
one can fairly represent the sound-speed perturbation as the product

$V(z,\,r)=b_{1}(r)Y_{1}(z)$ (17)

where $b_{1}(r)$ is a random function. For simplicity, it is assumed that
$b_{1}(r)$ is a stochastic process with the exponentially-decaying
autocorrelation function

$\left<b_{1}(r)b_{1}(r^{\prime})\right>=\exp(-\lvert
r-r^{\prime}\rvert/\bar{r}),$ (18)

where the correlation length $\bar{r}$ is taken of 10 km, that is typical for
the deep ocean Dozier and Tappert (1978b). Then the realizations of $b_{1}(r)$
can be computed via the formula $b_{1}=\sqrt{2\bar{r}}\eta$, where $\eta$ is a
solution of the Ornstein-Uhlenbeck stochastic differential equation

$\frac{d\eta}{dr}=-\frac{1}{\bar{r}}\eta(r)+\frac{1}{\bar{r}}\xi(r).$ (19)

The procedure of solving this equation is described, for instance, in Mallick
and Marcq (2002). Here $\xi$ is a Gaussian white noise satisfying

$\left<\xi(r)\right>=0,\quad\left<\xi(r)\xi(r^{\prime})\right>=\delta(r-r^{\prime}).$
(20)

The resulting function $b_{1}(r)$ satisfies $\left<b_{1}^{2}(r)\right>\simeq
1$.

The calculated function $Y_{1}(z)$ is depicted in Fig. 1(b). It involves step-
like changes in the depth interval from 100 to 300 meters. These changes are
caused by the depth discretization of the hydrological data and physically
irrelevant. Wave modeling is insensitive to them, but ray-based calculations
can be significantly affected. Therefore, in ray calculations we use a smooth
analytical approximation

$Y_{1}=Ay_{\mathrm{a}}\exp(-y_{\mathrm{a}}^{n})+B\exp(-y_{\mathrm{b}}^{2}),$
(21)

where $A=0.0027$, $y_{\mathrm{a}}=z/z_{\mathrm{a}}$, $n=1.1$,
$z_{\mathrm{a}}=24$ m, $B=2\cdot 10^{-5}$,
$y_{\mathrm{b}}=(z-z_{\mathrm{c}})/z_{\mathrm{b}}$, $z_{\mathrm{b}}=50$ m,
$z_{\mathrm{c}}=200$ m.

##  IV Finite-range evolution operator

Finite-range evolution operator (FREO) had been introduced for studying wave
propagation in a randomly-inhomogeneous waveguide in Virovlyansky  _et al._
(2012); Makarov  _et al._ (2010). Its quantum-mechanical analogue was earlier
considered in Kolovsky (1997). Basically, a finite-range evolution operator
(FREO) is an element of one-parameter group generated by the operator in the
right-hand side of the equation (3). Consider a solution $\Phi(z,\,r)$ of the
parabolic equation (3) complemented with the standard boundary conditions

$\left.\Phi\right|_{z=0}=0,\quad\left.\frac{d\Phi}{dz}\right|_{z=h}=0$ (22)

and the initial condition $\Phi(z,\,r=0)=\bar{\Phi}(z)$, where $\bar{\Phi}(z)$
belongs to $L^{2}[0,h]$ and satisfies (22). Then FREO $\hat{G}(\tau)$ is
defined on the subspace of $L^{2}[0,h]$ (restricted by (22)) as

$\hat{G}(\tau)\bar{\Phi}(z)\equiv\left.\Phi(z,\,r)\right|_{r=\tau}.$ (23)

By definition, the FREO describes transformation of a wavefield in course of
propagation along a finite waveguide segment of length $\tau$. Each
realization of inhomogeneity produces its own realization of the FREO. Our
interest is concerned with statistical properties of the FREO and their
connection to classical ray stability.

Note that the choice of the hard wall boundary condition at the bottom (22) is
typical for the deep ocean acoustics problems when the attention is restricted
to the trapped modes, whose propagation is not affected by the bottom
interaction. Under these conditions no energy is absorbed by the bottom.
Hence, if the attenuation in the sea water is negligible (this is true for the
sound frequencies of our interest) and refraction index in (3) has no
imaginary part, then the FREO is a unitary operator.

FREO can be represented as a matrix in the basis of normal modes $\phi_{j}(z)$
satisfying the Sturm-Liouville problem (12,22). Matrix elements of the FREO
are given by

$G_{mn}(\tau)=\int\limits_{0}^{h}\phi_{m}(z)\hat{G}(\tau)\phi_{n}(z)\,dz.$
(24)

Thus, the matrix elements $G_{mn}$ are complex-valued amplitudes of modal
transitions. For the range-independent waveguide, the matrix of FREO is
diagonal with $|G_{mm}|=1$.

Eigenvalues and eigenvectors of the FREO obey the equation

$\hat{G}\Psi_{m}(z,\,r)=g_{m}\Psi_{m}(z,\,r).$ (25)

Owing to unitarity, eigenvalues $g_{m}$ can be recast as

$g_{m}=e^{-ik_{0}\epsilon_{m}},\quad\epsilon_{m}\in\Re.$ (26)

Since eigenvalues of the FREO belong to the unit circle in the complex plane,
the FREO corresponds to the circular ensemble Stöckmann (2007). The FREO has
much in common with the Floquet operator governing wave propagation in a
range-periodic waveguide Smirnov  _et al._ (2005); Kon’kov  _et al._ (2007)
and quantum dynamics in time-periodic systems. For instance, quantity
$\epsilon_{m}$ is the analogue of quasienergy in quantum mechanics.

Note that eigenvalues of $\hat{G}(\tau)$ may be easily computed using its
matrix representation $G_{mn}(\tau)$. To accomplish this one has to clip a
finite block of this (infinite) matrix corresponding to the trapped modes,
neglecting their interaction with high-order modes. This simplification is
reasonable and does not affect accuracy of the eigenvalues computation (and
the numerics confirms that) since the prevailing small-angle propagation
corresponds to the low-order modes.

###  IV.1 Level spacing statistics

Wave chaos reveals itself in the spectrum of the FREO, in particular, in the
statistics of level spacings. A level spacing is defined as

$\begin{gathered}s=\frac{k_{0}M(\epsilon_{m+1}-\epsilon_{m})}{2\pi},\quad
m=1,2,\dots,M,\\\
\epsilon_{M+1}=\epsilon_{1}+\frac{2\pi}{k_{0}}.\end{gathered}$ (27)

where $\epsilon_{m}$ increases with increasing $m$, $M$ is the total number of
eigenvalues for a single realization of the FREO, equal to the number of
trapped modes.

Level spacing statistics can be studied in terms of the random matrix theory
Stöckmann (2007). Regular dynamics implies that the matrix of the FREO
consists of separate independent blocks. Then level sequences contributed from
different blocks are statistically independent, therefore, the resulting level
spacing distribution obeys the Poisson law

$\rho(s)\sim\exp(-s).$ (28)

Under conditions of ergodic chaos, all normal modes are coupled to each other.
This results in repulsion of neighbouring levels, the phenomenon closely
related to spectral splittings induced by tunneling Landau and Lifshitz
(1977). In this case level spacing statistics is described by the Wigner-Dyson
distribution

$\rho(s)\sim s^{\zeta}\exp\left(-Cs^{2}\right),$ (29)

where constants $\zeta$ and $C$ depend on symmetries of the FREO. As the FREO
doesn’t possess the symmetry $r\to-r$, it corresponds to the circular unitary
ensemble (CUE) with $\zeta=2$ and $C=4/\pi$ Kolovsky (1997).

Figure 3:  Berry-Robnik distribution with $v_{\mathrm{r}}=0.9$,
$v_{\mathrm{r}}=0.5$ and $v_{\mathrm{r}}=0.1$.

The most interesting case is mixed phase space, with the coexistence of
regular and chaotic domains. Then level spacing statistics should be described
by some combination of Poisson and Wigner-Dyson laws. In the short-wavelength
limit one can use the Berry-Robnik distribution Berry and Robnik (1984)

$\rho(s)=\left[v_{\mathrm{r}}^{2}\operatorname{erfc}\left(\frac{\sqrt{\pi}}{2}v_{\mathrm{c}}s\right)+\right.\\\
\left.+\left(2v_{\mathrm{r}}v_{\mathrm{c}}+\frac{\pi}{2}v_{\mathrm{c}}^{3}s\right)\exp\left(-\frac{\pi}{4}v_{\mathrm{c}}^{2}s^{2}\right)\vphantom{\frac{\sqrt{\pi}}{2}}\right]\exp(-v_{\mathrm{r}}s),$
(30)

where $v_{\mathrm{r}}$ and $v_{\mathrm{c}}$ are relative phase space volumes
corresponding to regular and chaotic ray motion, respectively,
$v_{\mathrm{r}}+v_{\mathrm{c}}=1$. In (30), it is assumed that phase space
consists of only two distinct domains: one regular and one chaotic.

Taking into account the correspondence between rays and modes established by
the WKB representation (11), one can regard $v_{\mathrm{r}}$ as a fraction of
modes whose refraction is regular, and $v_{\mathrm{c}}$ as a fraction of modes
exhibiting wave chaos. In the limiting cases $v_{\mathrm{r}}=1$ and
$v_{\mathrm{c}}=1$, Berry-Robnik formula (30) reduces to the Poisson
distribution and the Wigner-Dyson distribution for the orthogonal ensemble
($\zeta=1$), respectively. As is shown in Fig. 3, Berry-Robnik distribution
undergoes smooth transition from the Poisson to the Wigner-Dyson law as
$v_{\mathrm{r}}$ decreases from 1 to 0. Wigner-Dyson distributions for
orthogonal ($\zeta=1$) and unitary ($\zeta=2$) ensembles are close to each
other, therefore, one can use (30) as an approximation for the circular
unitary ensemble. Berry-Robnik distribution is based on the assumption of the
total statistical independence of the matrix blocks corresponding to regular
and chaotic dynamics. This assumption is completely fulfilled only in the
semiclassical limit. As wave corrections grow, independence degrades due to
regular-to-chaotic tunneling Bäcker  _et al._ (2005). Hence, Berry-Robnik
formula cannot work perfectly for low-frequency sound propagation. In this
case, a better fit is expected with the Brody distribution Prosen and Robnik
(1994)

$\rho(s)=(\beta+1)A_{\beta}s^{\beta}\exp(-A_{\beta}s^{\beta+1}),$ (31)

where $A_{\beta}=[\Gamma(\frac{\beta+2}{\beta+1})]^{\beta+1}$, $\Gamma$ is the
Euler gamma function. $\beta=0$ corresponds to the Poisson distribution,
$\beta=1$ yields the Wigner-Dyson distribution with $\zeta=1$. Unfortunately,
the Brody distribution is semiempirical and doesn’t have an explicit physical
interpretation in the intermediate regime $0<\beta<1$.

Functions (30) and (31) should describe level spacing statistics for single
realizations of the FREO. Very unfortunately, this is problematic because one
encounters insufficiency of the statistical ensemble. Indeed, long-range sound
propagation is feasible only with low acoustic frequencies of tens or hundreds
Hz. In this frequency range, number of trapped modes doesn’t exceed several
hundreds. To resolve this problem, we consider ensemble-averaged level spacing
distribution

$\rho(s,\,\tau)=\lim\limits_{N\to\infty}\frac{1}{N}\sum\limits_{n=1}^{N}P_{n}(s,\,\tau),$
(32)

where $P_{n}(s,\,\tau)$ is a level spacing distribution corresponding to
$n$-th realization of FREO. Fitting the function $\rho(s,\,\tau)$ with the
Berry-Robnik distribution (30), one can estimate number regularly propagating
modes for various values of $\tau$ and, whereby, track the transtion to chaos
with increasing $\tau$. However, it should be noted that formula (32) enables
accurate estimate of $v_{\mathrm{r}}$ only if fluctuations of $v_{\mathrm{r}}$
are weak, otherwise one should take into account nonlinearity of $\rho$ as a
function of $v_{\mathrm{r}}$ in (30).

Likewise we can observe transformation of level spacing statistics by fitting
(32) with the Brody distribution (31). Then the value of $\beta$ corresponding
to the best fit should increase with increasing $\tau$ from 0 to 1, reflecting
wavefield stochastization. Such calculations were performed in Makarov  _et
al._ (2010); Virovlyansky  _et al._ (2012) for a waveguide with a perturbed
biexponential sound-speed profile.

###  IV.2 Eigenfunction statistics

Wave chaos is also reflected in eigenfunction statistics of the FREO. Each
eigenfunction can be expressed as a superposition of normal modes:

$\Phi_{m}(z)=\sum\limits_{n}c_{mn}\phi_{n}(z),$ (33)

where $c_{mn}$ is the $m$-th component of $n$-th eigenvector of the matrix
$\hat{G}$, $\phi_{n}(z)$ is the $n$-th normal mode.

There are many methods for identification of “chaotic” eigenfunctions. In the
present work we use only one of them. Ray chaos can lead to intense energy
exchange between normal modes Virovlyansky and Zaslavsky (1999); Makarov  _et
al._ (2008); Hegewisch and Tomsovic (2012), therefore, a “chaotic”
eigenfunction is a compound of many modes. The stronger chaos, the larger
number of contributing modes. Hence, we can characterize “chaoticity” by
estimating number of principal components Varga and Pipek (2003) in the
expansion (33). Number of principal components is calculated as

$\nu(n)=\left(\sum\limits_{m=1}^{M}\lvert c_{mn}\rvert^{4}\right)^{-1}.$ (34)

Number of principal components $\nu$ is equal to 1 in an unperturbed waveguide
and grows as scattering intesifies, tending asymptotically to $M$.

##  V One-step Poincaré map

In the geometrical acoustics approximation, dynamics of a wavepacket is
represented as motion of some ray bundle. The bundle remains compact in course
of propagation if ray dynamics is regular. Under conditions of ray chaos the
bundle rapidly diverges. Initial conditions for rays contributing to the
bundle can be found by means of the Husimi function

$W_{\mathrm{h}}(p,\,z)=\biggl{|}\frac{1}{\sqrt[4]{2\pi\Delta_{z}^{2}}}\int
dz^{\prime}\Phi_{0}^{*}(z^{\prime})\times\\\
\exp\left[ikp(z^{\prime}-z)-\frac{(z^{\prime}-z)}{4\Delta_{z}^{2}}\right]\biggr{|}^{2},$
(35)

projecting an initial wavepacket $\Phi_{0}$ onto phase space of ray equations
(7). In the limit $k_{0}\to\infty$, the transformation (35) turns
eigenfunctions $\Phi_{m}$ of the FREO into some phase space sets which are
invariant under the shift $r=0\to r=\tau$. The procedure proposed in Makarov
_et al._ (2006) allows one to find out these sets in a randomly-inhomogeneous
waveguide. This is the one-step Poincaré map (or the specific Poincaré map).
The ray analogue of the map looks as follows:

$p_{i+1}=p(r=\tau|\,p_{i},z_{i}),\quad z_{i+1}=q(r=\tau|\,p_{i},z_{i}),$ (36)

where $p(r=\tau|\,p_{i},z_{i})$ and $z(r=\tau|\,p_{i},z_{i})$ are the
solutions of ray equations (7) with initial conditions $p(r=0)=p_{i}$,
$z(r=0)=z_{i}$. Values of $p$ and $z$, calculated at the $i$-th step of
mapping, become the initial conditions for the $(i+1)$-th step. This procedure
is equivalent to the usual Poincaré map Zaslavsky (2007) for a range-periodic
waveguide with the ray Hamiltonian

$\bar{H}=\frac{p^{2}}{2}+U(z)+\tilde{V}(z,\,r).$ (37)

Here $\tilde{V}(z,\,r)$ is periodic in $r$ function

$\tilde{V}(z,\,r^{\prime}+n\tau)=V(z,\,r^{\prime}),\quad 0\leqslant
r^{\prime}\leqslant\tau,$ (38)

where $n$ is an integer. As it follows from (38), $\tilde{V}(z,\,r)$ is a
sequence of identical pieces of $V(z,\,r)$, each of them has the length
$\tau$. Thus we replace the original randomly-perturbed Hamiltonian system by
an equivalent periodically-perturbed one. This replacement is valid as long as
we restrict ourselves by considering dynamics within the range interval
$[0:\tau]$. One-step Poincaré map can be considered as the classical
counterpart of the FREO.

Owing to analogy with the usual Poincaré map, the main property of the one-
step Poincaré map can be formulated as follows: each point of a continuous
closed ray trajectory of the map (36) corresponds to a starting point of the
solution of (7) which remains stable by Lyapunov till the range $r=\tau$. The
inverse statement is not, in general, true. Hence, the one-step Poincaré map
provides a sufficient but not necessary criterion of stability.

Map (36) was studied in Makarov  _et al._ (2009, 2006); Makarov and Uleysky
(2006); Gan  _et al._ (2010); Gan and Lei (2011). Here we shall give its brief
description. It is reasonable to make canonical transformation of ray
variables from momentum-depth $(p-z)$ to the action–angle $(I,\vartheta)$.
This procedure provides more appropriate representation of ray equations. The
angle variable $\vartheta$ canonically conjugated to the action (10) is given
by

$\vartheta=\frac{\partial}{\partial I}\int\limits_{z_{0}}^{z}p\,dz.$ (39)

The transformed ray Hamiltonian is written as

$\bar{H}=H_{0}(I)+\tilde{V}(I,\vartheta,r).$ (40)

Ray equations in terms of the new variables:

$\begin{gathered}\frac{dI}{dr}=-\frac{\partial
H}{\partial\vartheta}=-\frac{\partial V}{\partial\vartheta},\\\
\frac{d\vartheta}{dr}=\frac{\partial H}{\partial I}=\omega(I)+\frac{\partial
V}{\partial I},\end{gathered}$ (41)

where $\omega=2\pi/D$ is spatial frequency of a ray trajectory in a waveguide.
Perturbation $\tilde{V}(I,\vartheta)$ can be expanded into a double Fourier
series

$\tilde{V}=\frac{1}{2}\sum\limits_{k,k^{\prime}=1}^{\infty}V_{k,k^{\prime}}e^{i(k\vartheta-k^{\prime}\Omega
r)}+\text{c.~{}c.},$ (42)

where $\Omega=2\pi/\tau$. Inserting (42) into (41), we obtain

$\displaystyle\frac{dI}{dr}$
$\displaystyle=-\frac{i}{2}\sum\limits_{k,k^{\prime}=1}^{\infty}kV_{k,k^{\prime}}e^{i(k\vartheta-k^{\prime}\Omega
r)}+\text{c.~{}c.},$ (43) $\displaystyle\frac{d\vartheta}{dr}$
$\displaystyle=\omega+\frac{1}{2}\sum\limits_{k,k^{\prime}=1}^{\infty}\frac{\partial
V_{k,k^{\prime}}}{\partial I}e^{i(k\vartheta-k^{\prime}\Omega
r)}+\text{c.~{}c.}$

If the condition

$k^{\prime}D(I)=k\tau$ (44)

is fulfilled, there occurs resonance. The pair of integers $k^{\prime}$ and
$k$ determines multiplicity of resonance $k:k^{\prime}$. Resonances occur at
certain values of the action, which correspond to the so-called resonant
torus. Ray dynamics in a small vicinity of a resonant torus with ray action
$I_{0}$ can be described using the so-called resonance approximation Zaslavsky
(2007), when one leaves only resonant terms in the r.h.s. of (43). It should
be mentioned that one and the same resonant torus corresponds to an infinite
number of resonances with multiplicities $(jk):(jk^{\prime})$, where $j$ is an
integer. However, resonance Fourier amplitudes $V_{k,k^{\prime}}$ rapidly
decrease with increasing $k$ and $k^{\prime}$, therefore, only few low-order
resonances influence significantly ray dynamics. Consequently, we can take
into account only some finite number of dominant resonances. For further
simplification, we make the following procedures:

  1. 1.

As $\tilde{V}$ is a smooth function of $z$ in the underwater sound channel
considered, the derivative $d\tilde{V}/dI$ is small compared with $\omega$,
and the sum in the second equation of (43) can be dropped out.

  2. 2.

Near resonance, spatial frequency $\omega$ can be expanded as
$\omega=\Omega+\omega_{I}^{\prime}\Delta I$, where
$\omega_{I}^{\prime}=d\omega/dI$.

Then introducing new variables

$\Delta I=I-I_{0},\quad\psi=k\vartheta-k^{\prime}\Omega r,$ (45)

and expressing $V_{k,k^{\prime}}$ as $\lvert
V_{k,k^{\prime}}\rvert\exp(i\zeta_{k,k^{\prime}})$, we can rewrite (43) as

$\begin{gathered}\frac{d(\Delta I)}{dr}=\sum\limits_{l=1}^{L}lk\lvert
V_{lk,lk^{\prime}}\rvert\sin(l\psi+\zeta_{lk,lk^{\prime}})=-\frac{\partial\tilde{H}}{\partial\psi},\\\
\frac{d\psi}{dr}=\omega_{I}\Delta I=\frac{\partial\tilde{H}}{\partial(\Delta
I)},\end{gathered}$ (46)

where $L$ is the number of dominant resonances, and

$\tilde{H}=\omega_{I}^{\prime}\frac{(\Delta
I)^{2}}{2}+\sum\limits_{l=1}^{L}lk\lvert
V_{lk,lk^{\prime}}\rvert\cos(l\psi+\zeta_{lk,lk^{\prime}}).$ (47)

If $L=1$, (47) turns into the universal Hamiltonian of nonlinear resonance
Makarov  _et al._ (2009); Zaslavsky (2007), and a phase space portrait of Eqs.
(46) contains the domain of finite motion enclosed by the separatrix and
corresponding to trapping into resonance. Then maximal value $\Delta I$ on the
separatrix is determined as

$\Delta I_{\text{max}}=4\sqrt{\frac{k\lvert
V_{k,k^{\prime}}\rvert}{\omega_{I}^{\prime}}}$ (48)

The terms with $l>1$ may deform the pendulum-like phase space portrait and,
moreover, result in the presence of additional separatrices inside the domain
of finite motion. The latter phenomenon can occur when the perturbation
oscillates with depth Makarov  _et al._ (2006).

Transition to global chaos in the one-step Poincaré map happens when
neighbouring dominant resonances overlap. The criterion of overlapping is the
well-known Chirikov criterion

$\frac{\Delta I_{\text{max}}}{\delta I}\geqslant 1.$ (49)

Here $\delta I$ is the distance between neighbouring dominant resonances in
the action space. Its variability with $\tau$ for $\tau>D$ is described by
equation

$\delta I=\frac{2\pi}{\omega_{I}^{\prime}\tau},$ (50)

that is, increasing of $\tau$ enhances resonance overlapping. Differences in
phase space patterns corresponding to different realizations of the
perturbation are associated with phase and amplitude fluctuations of Fourier-
amplitudes $V_{k,k^{\prime}}$. However, the contribution of these fluctuations
is limited, therefore, the ratio of the phase space volumes corresponding to
regular and chaotic motion is mainly controlled by $\tau$ and weakly varies
from one realization to another (see below for an illustration). This property
allows one to consider the one-step Poincaré map as an useful tool for
studying randomly-driven dynamical systems of various physical origins.

Resonance approximation presented above fails under violation of the
nondegeneracy condition $\omega_{I}^{\prime}\neq 0$. In this case one should
use more sophisticated approaches, one of them is presented in Rypina  _et
al._ (2007). We shall address this issue in the end of this Section.

Figure 4: Ray phase space portraits constructed via the one-step Poincaré map
(36). Each column corresponds a single realization of the sound-speed
perturbation. Value of $\tau$ is indicated in the left lower corner.

Fig. 4 illustrates phase space portraits calculated using the one-step
Poincaré map with three different realizations of the sound-speed
perturbation. The jump of the derivative $dU/dz$ at $z=z_{0}$ leads to fast
growth of numerical error in ray calculations. Therefore, we replaced the
expression (8) for $U(z)$ by the smoothed function

$U(z)=U_{1}(z)+\frac{1}{2}\left[1+\tanh\frac{z-z_{0}}{\Delta}\right]U_{2}(z),$
(51)

where $\Delta=1$ m. Each of the phase space portraits represents a mixed phase
space structure consisted of regular and chaotic domains. Phase portraits with
the same $\tau$ mainly differ only in angular locations of regular islands,
whereas their overall structure is very similar. The main regular domain is
placed near the point $z=z_{0}$, $p=0$ and corresponds to flat rays
intersecting the horizontal plane with the smallest angles. This circumstance
deserves especial attention because stability of flat rays is not typical for
sound propagation in the deep ocean. Numerous experiments on long-range sound
propagation in the North-Eastern Pacific Ocean (see, for instance, Spiesberger
and Tappert (1996); Worcester  _et al._ (1994, 1999); Wage  _et al._ (2005);
Grigorieva  _et al._ (2009)) indicate on strong irregularity of flat near-
axial rays, associated with ray chaos Simmen  _et al._ (1997); Beron-Vera  _et
al._ (2003). The “deterministic” mechanism of near-axial chaos is ray
scattering on vertical resonances caused by small-scale depth oscillations of
the sound-speed perturbation Hegewisch  _et al._ (2005); Kon’kov  _et al._
(2007); Makarov  _et al._ (2008). These oscillations are contributed from
high-number modes of an internal-wave field. In the Sea of Japan, the effect
of the high-number modes is weak, therefore, the sound-speed perturbation can
be fairly described by equation (17), where depth dependence is given by a
smooth function $Y_{1}(z)$. Weakness of high-number internal-wave modes is a
peculiar feature of the Sea of Japan. It is caused by the specific form of the
buoyancy frequency profile Spindel  _et al._ (2003). In particular, the
waveguide for internal waves, determined by the buoyancy frequency profile, is
too narrow. Consequently, it doesn’t efficiently focus high-number modes of
low-frequency internal waves. As low-frequency internal waves give the
dominant contribution into a total internal-wave field, depth oscillations of
the perturbation are suppressed.

Resonance overlapping is enhancing as $\tau$ grows, and stable islands
eventually submerge into the chaotic sea. However, there is a small region of
stability that survives for distances of hundreds kilometers, transforming
into a chain of islands around the region near $z=z_{0}$, $p=0$. This chain
corresponds to the smooth minimum of the function $D(p_{0})$ depicted in Fig.
2. In the absence of inhomogeneity, this minimum gives rise to the so-called
shearless torus in phase space, producing a weakly-divergent beam. It is
recognized that shearless tori can possess extraordinary persistence to chaos
Rypina  _et al._ (2007); del Castillo-Negrete  _et al._ (1996); Budyansky  _et
al._ (2009); Uleysky  _et al._ (2010a, b), therefore, a weakly-divergent beam
can survive in the presence of a sound-speed perturbation. Hence, formation of
a weakly-divergent beam can be considered as a possible mechanism responsible
for unusual stability of near-axial rays, observed in experiments Bezotvetnykh
_et al._ (2009); Spindel  _et al._ (2003).

##  VI Spectral statistics of FREO for the underwater sound channel in the Sea
of Japan

###  VI.1 General remarks

Figure 5:  Fraction of the phase space volume corresponding to regular motion
vs distance $\tau$ for various frequencies.

This section is devoted to numerical modeling of the FREO in the Sea of Japan.
Each realization of the FREO was represented as a matrix in the basis of
normal modes being solutions of the Sturm-Liouville problem (12) with the
boundary conditions (22). Only purely-water modes which propagate without
reaching the bottom were taken into account. They were selected using the
criterion $E_{m}<U(z=h)$, where $E_{m}$ is the $m$-th eigenvalue of the Sturm-
Liouville problem (12). This criterion ensues from the WKB approximation for
normal modes Virovlyansky and Zaslavsky (1999). Number of trapped modes $M$
depends on sound frequency. It is equal to 72 for $f=100$ Hz, 179 for $f=250$
Hz, 259 for $f=360$ Hz, and 361 for $f=500$ Hz. Statistical ensembles of the
FREOs, corresponding to the frequencies of 250, 360 and 500 Hz, were
calculated with 100 realizations of the perturbation. The ensemble
corresponding to 100 Hz was calculated with 500 realizations. For each
realization, we constructed a family of the FREOs $\hat{G}(\tau)$, where
$\tau=5,10,15,\dots,350$ km.

Figure 6: The Brody parameter $\beta$ vs distance $\tau$.

###  VI.2 Eigenvalue statistics

We calculated the ensemble-averaged level spacing distribution $\rho(s,\tau)$
using the formula (32) and fitted it, for each value of $\tau$, by means of
the Berry-Robnik (30) and Brody (31) distributions. Thus, we obtained
dependencies of the regular phase space volume $v_{\mathrm{r}}$ and Brody
parameter $\beta$ on distance $\tau$. As is shown in Fig. 5, $v_{\mathrm{r}}$
rapidly decreases in the first 100–150 km. Then, it becomes almost constant.
This may indicate on the influence of the long-living stable islands in the
vicinity of the weakly-divergent beam. Notably, the curves corresponding to
250, 360 and 500 Hz are very close to each other, whereas the curve
corresponding to 100 Hz lies above them and undergoes strong fluctuations
which persist even with increasing number of realizations. It should be
mentioned that the above estimates of $v_{\mathrm{r}}$ have limited accuracy
because the assumptions underlying the formula (30) are satisfied only
approximately. Therefore, the results obtained using the Berry-Robnik
distribution are rather qualitative than quantitative, especially for low
sound frequencies.

Figure 7:  Relano parameter $\alpha$ vs distance $\tau$.

Approximation of level spacing statistics by means of the Brody distribution
leads to qualitatively similar results. They are presented in Fig. 6. For
frequencies 250, 360 and 500 Hz, the Brody parameter $\beta$ gradually grows
from 0 to 1, reflecting transformation from Poissonian to Wigner-like
statistics. In the case of $f=100$ Hz the growth is remarkably slower and
affected by strong fluctuations. The fluctuations of $v_{\mathrm{r}}$ and
$\beta$ may be induced by regular-to-chaotic tunneling and dynamical
localization, the phenomena whose influence on spectral statistics is still
not well understood.

In addition, we used the method of spectral analysis, developed by A. Relano
with coworkers in Relaño  _et al._ (2002). In this method, one firstly
constructs a series

$\delta_{n}=\sum\limits_{i=1}^{n}(s_{i}-\left<s\right>),$ (52)

where $n=1,2,\dots,N-1$, $N$ is the total number of eigenvalues. Then, making
discrete Fourier transform

$\bar{\delta}_{k}=\frac{1}{\sqrt{N}}\sum\limits_{n}\delta_{n}\exp\left(\frac{2\pi
ikn}{N}\right),$ (53)

one finds the power spectrum

$S(k)=\lvert\bar{\delta}_{k}\rvert^{2}.$ (54)

Generally, ensemble-averaged spectrum obeys a power low

$\left<S(k)\right>\sim k^{-\alpha}.$ (55)

Relano with coworkers found that regular dynamics corresponds to $\alpha=2$,
and global chaos results in $\alpha=1$. In the mixed regime, $\alpha$ takes on
an intermediate value between 1 and 2 Relaño (2008). Fig. 7 demonstrates that
$\alpha$ decreases with increasing $\tau$ for the frequencies of 360 and 500
Hz, indicating gradual transition to chaos. However, $\alpha$ varies
relatively slowly and remains near the middle value 1.5 for all distances
considered, despite of the marked changes in the classical phase space
portrait (see Fig. 4). Analogous dependencies for the frequencies of 100 and
250 Hz exhibit strong fluctuations and, therefore, are not presented in the
figure. This implies that the method developed in Relaño  _et al._ (2002);
Relaño (2008) provides good agreement only for relatively short wavelengths.

###  VI.3 Eigenfunction analysis

Analysis of eigenfunctions possesses some advantages as compared with analysis
of eigenvalues. The main advantage is the possibility to associate each
eigenfunction with some set of normal modes and, whereby, associate it with a
certain geometry of propagation. We can properly classify eigenfunctions,
taking into account the interplay with normal modes. Such classification can
be used for finding wavepacket configurations whose dynamics is expected to be
less or more regular. Proper classification can be obtained using the
parameter $\mu$ Smirnov  _et al._ (2005). It is defined as

$\mu=\sum\limits_{m=1}^{M}\lvert c_{m}\rvert^{2}m.$ (56)

In an unperturbed waveguide, only one normal mode contributes to each
eigenfunction, and $\mu$ coincides with the number of this mode. Taking into
account the quantization rule (11), we obtain the formula

$\left<I\right>=\frac{\mu}{k_{0}}+\frac{1}{2k_{0}}$ (57)

that gives mean action corresponding to an eigenfunction. According to (57),
the parameter $\mu$ determines phase space location of the eigenfunction and
can serve as its identificator.

Figure 8: Distribution of eigenfunctions in the $\mu$–$\nu$ plane, where the
parameter $\mu$ is given by (56), and $\nu$ is a number of principal
components (34). Distance values: (a) $\tau=10$ km, (b) $\tau=35$ km, (c)
$\tau=100$ km, (d) $\tau=350$ km. The sound frequency is 500 Hz.

Figures 8 and 9 illustrate eigenfunction distributions in the $\mu$–$\nu$
plane, where $\nu$ is number of principal components (34), for the frequencies
of 500 and 100 Hz, respectively.

Let us firstly consider Fig. 8 corresponding to the frequency of 500 Hz. Mode
coupling is relatively weak for small values of $\tau$, therefore, the
distribution is mainly concentrated near $\nu=1$. We want to focus attention
on the dense spots elongated in the $\nu$-direction. They correspond to
eigenfunctions reflecting mode-medium resonance Virovlyansky and Zaslavsky
(1999) being the wave counterpart of ray-medium resonance (44). Indeed, an
isolated ray-medium resonance produces oscillations of ray action inside the
interval $I_{0}-\Delta I_{\text{max}}\leqslant I\leqslant I_{0}+\Delta
I_{\text{max}}$, where $I_{0}$ is a resonance action, $\Delta I_{\text{max}}$
is resonance width in the action space. According to the principle of ray-mode
duality, these oscillations of action correspond to coherent transitions
between the normal modes whose numbers satisfy the inequality

$m_{0}-\Delta m\leqslant m\leqslant m_{0}+\Delta m,$ (58)

where $m_{0}=k_{0}I_{0}+1/2$, $\Delta m=k_{0}\Delta I_{\text{max}}$.
Resonance-induced modal transtions give rise to eigenfunctions with $\mu\simeq
m_{0}$ and $\nu$ varying from 1 to $\Delta m$. As long as resonance values of
the action, being determined by $\tau$, are the same for all realizations of
the perturbation, these eigenfunctions form vertically-elongated
concentrations of points in the $\mu$–$\nu$ plane. Location of mode-medium
resonances along the $\mu$-axis can be found using the formula

$k^{\prime}D(I=\left<I\right>)=k\tau,$ (59)

where $\left<I\right>$ is linked to $\mu$ by (57). Concentrations induced by
mode-medium resonance disappear with increasing $\tau$ due to overlapping of
mode-medium resonances and delocalization Berman and Kolovskiĭ (1992).
Delocalization leads to abrupt growth of number of principal components. It
eventually subjects all eigenfunctions in the interval between $\mu\simeq 100$
and $\mu\simeq 300$, resulting in the “boomerang” pattern in the $\mu$–$\nu$
plane, as illustrated in Fig. 8(d). Left and right ends of the “boomerang” are
formed by weakly scattered eigenfunctions. The left end corresponds to the
almost horizontal near-axial propagation, that is, its regularity can be
associated with long-living stable islands in the one-step Poincaré map.

Figure 9: The same as in Fig. 8, but for the frequency of 100 Hz. (a)
$\tau=10$ km, (b) $\tau=35$ km, (c) $\tau=100$ km, (d) $\tau=350$ km.

Eigenfunction distribution in the $\mu$–$\nu$ plane for the frequency of 100
Hz possesses a more complicated structure. It is exceptionally regular for
$\tau=10$ km and $\tau=35$ km, as shown in Figs. 9(a) and 9(b). Mode-medium
resonances reveal themselves as “stalagmites”. Each stalagmite is drawn by a
family of distinct weakly biased lines. Contours of the most pronounced
stalagmite for $\tau=35$ km are somewhat disordered and smeared. Some traces
of stalagmite-like patterns survive even for distances of hundreds kilometers,
despite of global overlapping of ray-medium resonances. Persistence of
stalagmites for large $\tau$ indicates on the presence of eigenstates
localized near periodic orbits of the one-step Poincaré map. As long as these
eigenstates correspond to nonspeading wavepackets, we can associate such
localization with suppresion of ray chaos and recovery of regular refraction
in the vicinities of the periodic orbits. An analogous phenomenon had been
reported in Kon’kov  _et al._ (2007) for range-periodic waveguides.

Besides of stalagmites, Figs. 9(a) and 9(b) illustrate the pattens in the form
of “bridges”. For instance, a pronounced “bridge” in the left part of Fig.
9(a) connects the points $\mu=5$, $\nu=1$ and $\mu=30$, $\nu=1$. The
eigenfunctions producing the “bridges” are consisted of normal modes
satisfying the condition

$k_{0}(E_{m}-E_{n})=\frac{2\pi l}{\tau},\quad m>n.$ (60)

The aforementioned “bridge” in the Fig. 9(a) satisfies (60) with $m=30$,
$n=5$, and $l=9$. Condition (60) is equivalent to quantum resonance between
two energy levels. In the ray limit

$k_{0}(E_{m}-E_{n})\to(m-n)\frac{dE}{dI}\equiv\frac{2\pi(m-n)}{D},$ (61)

and condition (60) reduces to (44). Each realization of the FREO can yield
eigenfunctions manifesting resonance (60). If resonance (60) corresponding to
some numbers ($l$, $m$, $n$) is localized, i. e. the modes $m$ and $n$ are not
affected by other resonances, then the respective eigenfunction is a
superposition of normal modes $m$ and $n$,

$\Phi_{\text{res}}(z)\simeq c_{m}\phi_{m}+c_{n}\phi_{n},\quad\lvert
c_{m}\rvert^{2}+\lvert c_{n}\rvert^{2}\simeq 1.$ (62)

As this takes place, the ratio of amplitudes $\lvert c_{m}\rvert/\lvert
c_{n}\rvert$ is determined by the phase of the resonance harmonics of the
perturbation. For the perturbation (17), amplitude of the $l$-th resonance
harmonics reads

$B_{l}=\frac{1}{2\pi}\int\limits_{0}^{2\pi/\tau}b_{1}(r)\exp\left(-i\frac{2\pi
lr}{\tau}\right)dr.$ (63)

Phase of the resonance harmonics is a random quantity with uniform probability
density in the range $[0:2\pi]$. Each value of the phase uniquely determines
the values of $\mu$ and $\nu$ via the formulae

$\mu=\lvert c_{m}\rvert^{2}m+\lvert c_{n}\rvert^{2}n,$ (64) $\nu=\left(\lvert
c_{m}\rvert^{4}+\lvert c_{n}\rvert^{4}\right)^{-1}.$ (65)

It turnes out that quantities $\mu$ and $\nu$ are correlated for each
eigenfunction corresponding to localized resonance (60). This results in
formation of “bridges” in the $\mu$–$\nu$ plane. Localization becomes violated
with increasing $\tau$, therefore, correlation between $\nu$ and $\mu$ ceases,
and ordered “bridges” transform into disordered clouds of points.

Figure 10: Mean number of principal components as function of distance.

The above analysis shows that there are qualitative differences in sound
scattering for different frequencies. To estimate these differences
quantitatively, we need a suitable parameter characterizing scattering
strength. For example, chaos-induced phase space delocalization can be
measured by the ensemble-averaged number of principal components, divided by
the total number of trapped modes. As it follows from Fig. 10, the rate of
delocalization increases with increasing sound frequency. This implies that
chaotic diffusion associated with ray chaos degrades with increasing sound
wavelength. It can be thought of as a manifestation of dynamical localization
Stöckmann (2007), an analogue of Anderson localization, when destructive
interference supresses wavepacket spreading.

Figure 11: Fraction of strongly-localized eigenfunctions as function of
distance. The criterion of strong localization is the inequality $\nu\leqslant
2$.

In practice, it is useful to know fraction of the eigenfunction ensemble
corresponding to regular propagation. This quantity can be regarded as an
analogue of the parameter $v_{\mathrm{r}}$ in the Berry-Robnik distribution.
It can be estimated by means of the cumulative distribution function

$F(\nu)=\int\limits_{1}^{\nu}\rho(\nu^{\prime})\,d\nu^{\prime},$ (66)

where $\rho(\nu^{\prime})$ is the probability density function of $\nu$. We
can conditionally distinguish two regimes of localization: strong localization
and moderate localization. Strong localization implies that eigenfunction of
FREO is close to one of normal modes of the unperturbed waveguide. To select
strongly-localized eigenfunctions, we can use the inequality $\nu\leqslant 2$.
Dependence of $F(2)$ on the horizontal distance $\tau$ is depicted in Fig. 11.
Evidently, fraction of strongly-localized eigenfunctions is much larger for
$f=100$ Hz than for higher frequencies. The curves corresponding to the higher
frequencies are close to each other. They drop down in the range of small
$\tau$, then decreasing of $F(2)$ becomes very slow. Slow decreasing can be
linked to the presence of the long-living islands of stability in the
neibourhood of the weakly-divergent beam.

Figure 12: Fraction of moderately localized eigenfunctions vs distance. The
criterion of moderate localization is the inequality $\nu\leqslant 0.1M$,
where $M$ is the number of trapped modes.

In the regime of moderate localization, mode coupling can be sufficiently
strong, but an eigenfunction occupies relatively small phase space volume,
that is, number of principal components is limited. We use the inequality
$\nu\leqslant 0.1M$ as the criterion of moderate localization. As is
demonstrated in Fig. 12, almost all eigenfunctions corresponding to the
frequency of 100 Hz are moderately-localized. This is not the case of higher
frequencies, when fraction of moderately-localized eigenfunctions
significantly decreases with $\tau$. We can conclude that decreasing of sound
frequency can result in remarkable suppression of wavepacket diffusion in
phase space.

##  VII Conclusion

Spectral analysis of the finite-range evolution operator prompts a way to
explore wave dynamics in a randomly-inhomogeneous waveguide by means of the
quasideterministic approach, involving resonances, periodic orbits, phase
space portraits, e.t.c. This approach was originally proposed in Virovlyansky
_et al._ (2012); Makarov  _et al._ (2010). The present work is devoted to its
further development. We demonstrate various methods of spectral analysis. For
instance, fitting of level spacing statistics by means of the Berry-Robnik
distribution yields approximate estimate for fraction of regularly propagating
normal modes of a waveguide. Also we use the method developed by A. Relano
with coworkers, and consider distribution of eigenfunctions in the $\mu$–$\nu$
space. In our opinion, the latter approach is the most robust, albeit its
analytical description is lacking. The most important advantage of the
eigenfunction analysis is the possibility to study separately scattering of
different modes of a waveguide and find out the modes corresponding to regular
propagation. A detailed view of the eigenfunction distribution in the
$\mu$–$\nu$ space suggests that the mechanism of the chaos onset with
increasing distance can be associated with overlapping and delocalization of
mode-medium resonances. The approach based on the statistical analysis of
level spacings by means of the Berry-Robnik formula gives qualitative
description of the transition but doesn’t ensure quantitative agreement. The
Relano method and the method based on the Brody distribution also give the
qualitative description but cannot make any quantitative estimates due to
their semiempirical nature.

We consider the underwater sound channel in the Sea of Japan as an example.
Our analysis shows that almost horizontal near-axial sound propagation
preserves regularity over distances of hundreds kilometers. There are two
factors responsible for near-axial stability. The first factor is the peculiar
hydrological structure in the region considered, resulting in the absence of
vertical oscillations of the sound-speed perturbation. This circumstance leads
to a qualitatively different scenario of ray chaos, as compared with the well-
known acoustic experiments in the North-Eastern Pacific Ocean. The second
factor is the formation of the weakly-divergent beam supported by long-living
stable islands in classical phase space. Weakly-divergent beam occurs in the
vicinity of the shearless torus. Non-dispersive motion of quantum wavepackets
near shearless tori was earlier observed in Kudo and Monteiro (2008). Also, we
should emphasize that sound propagation with the frequency of 100 Hz reveals
exceptionally high degree of stability, indicating on different physics of
scattering for low sound frequencies. This phenomenon can be associated with
strong dynamical localization.

We suppose that applicability of the approach presented in this paper is not
limited to the problems of wave propagation in random media. An analogue of
the FREO can be readily used for studying of noise-induced quantum transport
and related phenomena. We suggest that combination of deterministic and
statistical approaches should provide an insightful view on details of
dynamics, especially in the presence of intermittency or synchronization.

This work was supported by the Russian Foundation of Basic Research (projects
09-05-98608 and 09-02-01258-a), the Federal Program “World Ocean” and the
“Dynasty” Foundation. Authors are grateful to A.R. Kolovsky, S. Tomsovic, O.A.
Godin, K.V. Koshel and V.V. Novotryasov for helpful comments concerning the
subject of the research.

## References

  * Shockley  _et al._ (1982) R. C. Shockley, J. Northrop, P. G. Hansen, and C. Hartdegen, J. Acoust. Soc. Am. 71, 51 (1982).
  * Munk  _et al._ (1988) W. H. Munk, W. C. O’Reilly, and J. L. Reid, J. Phys. Oceanography 18, 1876 (1988).
  * Makarov  _et al._ (2009) D. Makarov, S. Prants, A. Virovlyansky, and G. Zaslavsky,  _Ray and wave chaos in ocean acoustics: chaos in waveguides_ , Series on complexity, nonlinearity and chaos (World Scientific, Singapore, 2009) p. 388.
  * Smith  _et al._ (1992) K. B. Smith, M. G. Brown, and F. D. Tappert, J. Acoust. Soc. Am. 91, 1939 (1992).
  * Smirnov  _et al._ (2001) I. P. Smirnov, A. L. Virovlyansky, and G. M. Zaslavsky, Phys. Rev. E 64, 036221 (2001).
  * Brown  _et al._ (2003) M. G. Brown, J. A. Colosi, S. Tomsovic, A. L. Virovlyansky, M. A. Wolfson, and G. M. Zaslavsky, J. Acoust. Soc. Am. 113, 2533 (2003).
  * Makarov  _et al._ (2004) D. V. Makarov, M. Y. Uleysky, and S. V. Prants, Chaos 14, 79 (2004).
  * Beron-Vera and Brown (2004) F. J. Beron-Vera and M. G. Brown, J. Acoust. Soc. Am. 115, 1068 (2004).
  * Virovlyansky and Zaslavsky (1999) A. L. Virovlyansky and G. M. Zaslavsky, Phys. Rev. E 59, 1656 (1999).
  * Smirnov  _et al._ (2004) I. P. Smirnov, A. L. Virovlyansky, and G. M. Zaslavsky, Chaos 14, 317 (2004).
  * Smirnov  _et al._ (2005) I. P. Smirnov, A. L. Virovlyansky, M. Edelman, and G. M. Zaslavsky, Phys. Rev. E 72, 026206 (2005).
  * Kon’kov  _et al._ (2007) L. E. Kon’kov, D. V. Makarov, E. V. Sosedko, and M. Y. Uleysky, Phys. Rev. E 76, 056212 (2007).
  * Makarov  _et al._ (2008) D. V. Makarov, L. E. Kon’kov, and M. Y. Uleysky, Acoust. Phys. 54, 382 (2008).
  * Virovlyansky  _et al._ (2012) A. L. Virovlyansky, D. V. Makarov, and S. V. Prants, Phys. Usp. 55 (2012).
  * Hegewisch  _et al._ (2005) K. C. Hegewisch, N. R. Cerruti, and S. Tomsovic, J. Acoust. Soc. Am. 117, 1582 (2005).
  * Munk and Wunsch (1979) W. Munk and C. Wunsch, Deep Sea Res. Part A. 26, 123 (1979).
  * Tappert and Tang (1996) F. D. Tappert and X. Tang, J. Acoust. Soc. Am. 99, 185 (1996).
  * Zaslavsky (2007) G. M. Zaslavsky,  _The physics of chaos in Hamiltonian systems_ (Imperial College Press, London, 2007) p. 328.
  * Dozier and Tappert (1978a) L. B. Dozier and F. D. Tappert, J. Acoust. Soc. Am. 63, 353 (1978a).
  * Colosi and Morozov (2009) J. A. Colosi and A. K. Morozov, J. Acoust. Soc. Am. 126, 1026 (2009).
  * Wolfson and Tomsovic (2001) M. A. Wolfson and S. Tomsovic, J. Acoust. Soc. Am. 109, 2693 (2001).
  * Makarov  _et al._ (2006) D. V. Makarov, M. Y. Uleysky, M. V. Budyansky, and S. V. Prants, Phys. Rev. E 73, 066210 (2006).
  * Makarov  _et al._ (2010) D. V. Makarov, L. E. Kon’kov, and M. Y. Uleysky, ArXiv e-prints (2010), arXiv:1008.3037 [nlin.CD] .
  * Hegewisch and Tomsovic (2012) K. C. Hegewisch and S. Tomsovic, Europhys. Lett. 97, 34002 (2012).
  * Stöckmann (2007) H. J. Stöckmann,  _Quantum Chaos: An Introduction_ (Cambridge University Press, Cambridge, 2007) p. 384.
  * Kolovsky (1997) A. R. Kolovsky, Phys. Rev. E 56, 2261 (1997).
  * Bezotvetnykh  _et al._ (2009) V. Bezotvetnykh, A. Burenin, Y. Morgunov, and Y. Polovinka, Acoust. Phys. 55, 376 (2009).
  * Spindel  _et al._ (2003) R. C. Spindel, J. Na, P. H. Dahl, S. Oh, C. Eggen, Y. G. Kim, V. A. Akulichev, and Y. N. Morgunov, IEEE J. Ocean. Engin. 28, 297 (2003).
  * (29) “Oceanography and marine environment of the Far Eastern Region of Russia (proj. leader Rostov I. D.),” http://www.pacificinfo.ru/en.
  * Brekhovskikh  _et al._ (1990) L. M. Brekhovskikh, V. V. Goncharov, S. A. Dremuchev, V. M. Kurtepov, V. G. Selivanov, and Y. A. Chepurin, Sov. Phys. Acoust. 36, 461 (1990).
  * Smirnov  _et al._ (1999) I. Smirnov, J. Caruthers, and A. Khil’ko, Radiophys. Quantum Electron. 42, 864 (1999).
  * Morozov and Colosi (2005) A. K. Morozov and J. A. Colosi, J. Acoust. Soc. Am. 117, 1611 (2005).
  * Petukhov (2009) Y. Petukhov, Acoust. Phys. 55, 785 (2009).
  * Colosi and Brown (1998) J. A. Colosi and M. G. Brown, J. Acoust. Soc. Am. 103, 2232 (1998).
  * LeBlanc and Middleton (1980) L. R. LeBlanc and F. H. Middleton, J. Acoust. Soc. Am. 67, 2055 (1980).
  * Dozier and Tappert (1978b) L. B. Dozier and F. D. Tappert, J. Acoust. Soc. Am. 64, 533 (1978b).
  * Mallick and Marcq (2002) K. Mallick and P. Marcq, Phys. Rev. E 66, 041113 (2002).
  * Landau and Lifshitz (1977) L. D. Landau and E. M. Lifshitz,  _Course of theoretical physics._ , Vol. 3: Quantum mechanics. Nonrelativistic theory (Pergamon Press, Oxford, 1977).
  * Berry and Robnik (1984) M. V. Berry and M. Robnik, J. Phys. A: Math. Gen. 17, 2413 (1984).
  * Bäcker  _et al._ (2005) A. Bäcker, R. Ketzmerick, and A. G. Monastra, Phys. Rev. Lett. 94, 054102 (2005).
  * Prosen and Robnik (1994) T. Prosen and M. Robnik, J. Phys. A: Math. Gen. 27, 8059 (1994).
  * Varga and Pipek (2003) I. Varga and J. Pipek, Phys. Rev. E 68, 026202 (2003).
  * Makarov and Uleysky (2006) D. Makarov and M. Uleysky, J. Phys. A: Math. Gen. 39, 489 (2006).
  * Gan  _et al._ (2010) C. Gan, Q. Wang, and M. Perc, J. Physics A: Math. Theor. 43, 125102 (2010).
  * Gan and Lei (2011) C. Gan and H. Lei, J. Sound Vibr. 330, 2174 (2011).
  * Rypina  _et al._ (2007) I. I. Rypina, M. G. Brown, F. J. Beron-Vera, H. Koçak, M. J. Olascoaga, and I. A. Udovydchenkov, Phys. Rev. Lett. 98, 104102 (2007).
  * Spiesberger and Tappert (1996) J. L. Spiesberger and F. D. Tappert, J. Acoust. Soc. Am. 99, 173 (1996).
  * Worcester  _et al._ (1994) P. F. Worcester, B. D. Cornuelle, J. A. Hildebrand, J. William S. Hodgkiss, T. F. Duda, J. Boyd, B. M. Howe, J. A. Mercer, and R. C. Spindel, J. Acoust. Soc. Am. 95, 3118 (1994).
  * Worcester  _et al._ (1999) P. F. Worcester, B. D. Cornuelle, M. A. Dzieciuch, W. H. Munk, B. M. Howe, J. A. Mercer, R. C. Spindel, J. A. Colosi, K. Metzger, T. G. Birdsall, and A. B. Baggeroer, J. Acoust. Soc. Am. 105, 3185 (1999).
  * Wage  _et al._ (2005) K. E. Wage, M. A. Dzieciuch, P. F. Worcester, B. M. Howe, and J. A. Mercer, J. Acoust. Soc. Am. 117, 1565 (2005).
  * Grigorieva  _et al._ (2009) N. S. Grigorieva, G. M. Fridman, J. A. Mercer, R. K. Andrew, M. A. Wolfson, B. M. Howe, and J. A. Colosi, J. Acoust. Soc. Am. 125, 1919 (2009).
  * Simmen  _et al._ (1997) J. Simmen, S. M. Flatte, and G.-Y. Wang, J. Acoust. Soc. Am. 102, 239 (1997).
  * Beron-Vera  _et al._ (2003) F. J. Beron-Vera, M. G. Brown, J. A. Colosi, S. Tomsovic, A. L. Virovlyansky, M. A. Wolfson, and G. M. Zaslavsky, J. Acoust. Soc. Am. 114, 1226 (2003).
  * del Castillo-Negrete  _et al._ (1996) D. del Castillo-Negrete, J. M. Greene, and P. J. Morrison, Physica D 91, 1 (1996).
  * Budyansky  _et al._ (2009) M. V. Budyansky, M. Y. Uleysky, and S. V. Prants, Phys. Rev. E 79, 056215 (2009).
  * Uleysky  _et al._ (2010a) M. Y. Uleysky, M. V. Budyansky, and S. V. Prants, Phys. Rev. E 81, 017202 (2010a).
  * Uleysky  _et al._ (2010b) M. Uleysky, M. Budyansky, and S. Prants, J. of Exp. and Theor. Phys. 111, 1039 (2010b).
  * Relaño  _et al._ (2002) A. Relaño, J. M. G. Gómez, R. A. Molina, J. Retamosa, and E. Faleiro, Phys. Rev. Lett. 89, 244102 (2002).
  * Relaño (2008) A. Relaño, Phys. Rev. Lett. 100, 224101 (2008).
  * Berman and Kolovskiĭ (1992) G. P. Berman and A. R. Kolovskiĭ, Sov. Phys. Usp. 35, 303 (1992).
  * Kudo and Monteiro (2008) K. Kudo and T. S. Monteiro, Phys. Rev. E 77, 055203 (2008).