Paraxial coupling of propagating modes in threedimensional waveguides with random boundaries
Abstract
We analyze long range wave propagation in threedimensional random waveguides. The waves are trapped by top and bottom boundaries, but the medium is unbounded in the two remaining directions. We consider scalar waves, and motivated by applications in underwater acoustics, we take a pressure release boundary condition at the top surface and a rigid bottom boundary. The wave speed in the waveguide is known and smooth, but the top boundary has small random fluctuations that cause significant cumulative scattering of the waves over long distances of propagation. To quantify the scattering effects, we study the evolution of the random amplitudes of the waveguide modes. We obtain that in the long range limit they satisfy a system of paraxial equations driven by a Brownian field. We use this system to estimate three important modedependent scales: the scattering mean free path, the crossrange decoherence length and the decoherence frequency. Understanding these scales is important in imaging and communication problems, because they encode the cumulative scattering effects in the wave field measured by remote sensors. As an application of the theory, we analyze time reversal and coherent interferometric imaging in strong cumulative scattering regimes.
Key words. Waveguides, random media, asymptotic analysis.
AMS subject classifications. 76B15, 35Q99, 60F05.
1 Introduction
We study long range scalar (acoustic) wave propagation in a threedimensional waveguide. The setup is illustrated in Figure LABEL:fig:schematic, and it is motivated by problems in underwater acoustics. We denote by the range, the main direction of propagation of the waves. The medium is unbounded in the crossrange direction , but it is confined in depth by two boundaries which trap the waves, thus creating the waveguide effect.
The acoustic pressure field is denoted by , and it satisfies the wave equation
\hb@xt@.01(1.1) 
in a medium with wave speed . The excitation is due to a source located in the plane , emitting the pulse . The medium is quiescent before the source excitation,
\hb@xt@.01(1.2) 
The bottom of the waveguide is assumed rigid
\hb@xt@.01(1.3) 
and we take a pressure release boundary condition at the perturbed top boundary
\hb@xt@.01(1.4) 
Perturbed means that the boundary has small fluctuations around the mean depth ,
\hb@xt@.01(1.5) 
We choose this setup for simplicity. The results extend readily to other boundary conditions and to fluctuating bottoms. Such boundaries were considered recently in [1, 11], in twodimensional waveguides. Extensions to media with small dependent fluctuations of the wave speed can also be made using the techniques developed in [12, 6, 8, 9, 10].
The goal of our study is to quantify the effect of scattering at the surface. Because in applications it is not feasible to know the boundary fluctuations in detail, we model them with a random process. The solution of equations (LABEL:eq:form.1)(LABEL:eq:form.4) is therefore a random field, and we describe in detail its statistics at long ranges, where cumulative scattering is significant. We use the results for two applications: time reversal and sensor array imaging.
Our method of solution uses a change of coordinates to straighten the boundary. The transformed problem has a simple geometry but a randomly perturbed differential operator. Its solution is given by a superposition of propagating and evanescent waveguide modes, with random amplitudes. We show that in the long range limit these amplitudes satisfy a system of paraxial equations that are driven by a Brownian field. The detailed characterization of the statistics of follows from this system. It involves the calculation of the modedependent scattering mean free path, which is the distance over which the modes lose coherence; the modedependent decoherence length, which is the crossrange offset over which the mode amplitudes decorrelate; and the modedependent decoherence frequency, which is the frequency offset over which the mode amplitudes decorrelate. These scales are important in studies of time reversal and imaging, because they dictate the resolution of focusing and the robustness (statistical stability) of the results with respect to realizations of the random fluctuations of the boundary.
The paper is organized as follows: We begin in section LABEL:sect:HOMWG with the description of the reference pressure field in ideal waveguides with planar boundaries. The random field derived in section LABEL:sect:RANDWG may be viewed as a perturbation of , in the sense that it is decomposed in the same waveguide modes. However, the amplitudes of the modes are random and coupled. Because the fluctuations of the boundary are small, we consider in section LABEL:sect:LR a long range limit, so that we can observe significant cumulative scattering. The statistics of the wave field at such long ranges is described in section LABEL:sect:STAT. The results are summarized in section LABEL:sect:FM and are used in sections LABEL:sect:TR and LABEL:sect:IM for analyzing time reversal and imaging with sensor arrays. We end with a summary in section LABEL:sect:SUM.
2 Wave propagation in ideal waveguides
The pressure field in ideal waveguides, with planar boundaries, is given by
\hb@xt@.01(2.1) 
with Fourier coefficients satisfying a separable problem for the Helmholtz equation
\hb@xt@.01(2.2) 
with boundary conditions
\hb@xt@.01(2.3) 
and outgoing radiation conditions at . The Fourier transform of the source
\hb@xt@.01(2.4) 
is compactly supported in , for any and . Here is the central frequency and is the bandwidth.
2.1 Propagating and evanescent modes
The solution of the Helmholtz equation (LABEL:eq:HOM.2) is a superposition of propagating modes, and infinitely many evanescent ones,
\hb@xt@.01(2.5) 
The decomposition is in the orthonormal basis of the eigenfunctions of the selfadjoint differential operator in ,
with eigenvalues that are simple [16].
To simplify the analysis, we assume in this paper that the wave speed is homogeneous
\hb@xt@.01(2.6) 
The results for variable wave speeds are similar in all the essential aspects. The simplification brought by (LABEL:eq:HOM.6) amounts to having explicit expressions of the eigenfunctions, which are independent of the frequency
\hb@xt@.01(2.7) 
The eigenvalues are
\hb@xt@.01(2.8) 
where is the wavenumber, and only the first of them are nonnegative
\hb@xt@.01(2.9) 
The notation stands for the integer part. We suppose for simplicity that remains constant in the bandwidth , and write from now on .
The propagating components in (LABEL:eq:HOM.4) satisfy the twodimensional Helmholtz equation
\hb@xt@.01(2.10) 
with outgoing, radiation conditions at . The evanescent components solve
\hb@xt@.01(2.11) 
with decay condition at Here we introduced the coefficients of the source profile in the basis of the eigenfunctions
\hb@xt@.01(2.12) 
and the mode wavenumbers
\hb@xt@.01(2.13) 
We assume that none of the vanishes in the bandwidth, so that there are no standing waves. That is to say,
\hb@xt@.01(2.14) 
2.2 The paraxial regime
We now introduce the paraxial scaling for the ideal waveguide, with the source emitting a beam that propagates along the axis. As we show below, this happens when the crossrange profile of the source is larger than the wavelength. The source generates a quasiplane wave, with slowly varying envelope satisfying a Schrödingerlike equation.
Explicitly, we assume that the source is of the form
\hb@xt@.01(2.15) 
where is a small dimensionless parameter defined as the ratio of the central wavelength and the transverse width of the source. Standard diffraction theory gives that the Rayleigh length for a beam with initial width is of the order of
The Rayleigh length is the distance along the axis from the beam waist to the place where the beam area is doubled by diffraction. Therefore, we look at the wavefield at crossrange scales, similar to , and at range scale, similar to the Rayleigh length. We rename the field in this scaling as
\hb@xt@.01(2.16) 
The Fourier coefficients of (LABEL:eq:poeps) are given by the scaled version of (LABEL:eq:HOM.4)
\hb@xt@.01(2.17) 
with propagating mode amplitudes satisfying the scaled equation (LABEL:eq:HOM.10), with the source replaced by . They can be written as
in terms of the outgoing Green’s function
Here is the Hankel function of the first kind, and because , we can use its asymptotic form for a scaled range
The propagating components of the wave field become
with
\hb@xt@.01(2.18) 
for . The evanescent components are obtained similarly from (LABEL:eq:HOM.11)
and
\hb@xt@.01(2.19) 
for . These modes are exponentially damped and can be neglected.
In summary, the paraxial approximation of the wave field is given by
\hb@xt@.01(2.20) 
It is a superposition of forward going modes with complex valued amplitudes given by (LABEL:eq:HOM.15), and solving the paraxial equations
\hb@xt@.01(2.21) 
with initial conditions
\hb@xt@.01(2.22) 
3 Wave propagation in random waveguides
In this section we consider a waveguide with fluctuating boundary and analyze the wave field under the following scaling assumptions:

The transverse width of the source and the central wavelength satisfy
\hb@xt@.01(3.1) as in the previous section. The correlation length of the boundary fluctuations is similar to ,
\hb@xt@.01(3.2) so that there is a nontrivial interaction between the boundary fluctuations and the wavefield. Here is the scaled orderone correlation length defined below.

The scale of the propagation distance is much larger than . More precisely,
\hb@xt@.01(3.3) Recall that the Rayleigh length for a beam with initial width and central wavelength is of the order of in absence of random fluctuations. The highfrequency scaling assumption (LABEL:eq:sc.3) ensures that the propagation distance is similar to the Rayleigh length.

The amplitude of the boundary fluctuations is small, of the order of . As we will show, this scaling is precisely the one that gives a cumulative scattering effect of order one after the propagation distance .
We use the hyperbolicity of the problem to truncate mathematically the boundary fluctuations to the range interval . The bound is the maximum range of the fluctuations that can affect the waves up to the observation time of order . The lower bound in the range interval coincides with the location of the source. It is motivated by two facts: First, we observe the waves at positive ranges. Second, the backscattered field is negligible in the scaling regime defined above, as we show later in section LABEL:sect:LR.2.
The boundary fluctuations are modeled with a random process
\hb@xt@.01(3.4) 
The process is bounded, zeromean, stationary and mixing, meaning in particular that its covariance is integrable^{1}^{1}1More precisely, is a mixing process, with , as stated in [13, 4.6.2].. Because our method of solution flattens the boundary by changing coordinates, we require that is twice differentiable, with almost surely bounded derivatives. Its covariance function is given by
\hb@xt@.01(3.5) 
and we denote by its integral over ,
\hb@xt@.01(3.6) 
Our assumption on the differentiability of implies that is four times differentiable. Note that is the maximum of the integrated covariance , so we have
\hb@xt@.01(3.7) 
We define the scaled square amplitude and correlation length of the boundary fluctuations through the equations
\hb@xt@.01(3.8) 
3.1 Change of coordinates
We introduce the change of coordinates from to , with
\hb@xt@.01(3.9) 
It straightens the boundary to , for any and . The pressure field in the new coordinates is denoted by
\hb@xt@.01(3.10) 
It satisfies the simple boundary conditions
\hb@xt@.01(3.11) 
and the partial differential equation
\hb@xt@.01(3.12) 
derived from (LABEL:eq:form.1) and (LABEL:eq:RAND.6) using the chain rule. Here and are the gradient and Laplacian operators in and is the source of the form (LABEL:scaledsource).
When substituting the model (LABEL:eq:RAND.1) in (LABEL:eq:RAND.8), we obtain that satisfies a randomly perturbed problem
\hb@xt@.01(3.13) 
The higherorder terms denoted by the dots are
They come from the expansions in of the coefficients in (LABEL:eq:RAND.8), and are negligible in the limit considered in section LABEL:sect:LR.
3.2 Wave decomposition
Equation (LABEL:eq:RAND.9) is not separable, but we can still write its solution in the basis of the eigenfunctions (LABEL:eq:HOM.7). The expansion is similar to (LABEL:eq:HOM.4)
\hb@xt@.01(3.14) 
We define the forward and backward going wave mode amplitudes and by
\hb@xt@.01(3.15) 
so that the complex valued amplitudes of the propagating modes can be written as
Definition (LABEL:eq:def.ab) implies that
\hb@xt@.01(3.16) 
This equation is needed to specify uniquely the propagating mode amplitudes, because they each satisfy a single boundary condition in the range of the fluctuations. To derive these boundary conditions, let us observe that and must be constant in and in , because the boundary is flat outside . Moreover, the radiation conditions
imply that the mode amplitudes satisfy
\hb@xt@.01(3.17)  
\hb@xt@.01(3.18) 
The last equation is the boundary condition for . The boundary value follows from the jump conditions across the plane of the source in equation (LABEL:eq:RAND.9). We have
with defined by (LABEL:eq:HOM.14). This gives
and therefore
\hb@xt@.01(3.19) 
Substituting (LABEL:eq:RAND.10) in (LABEL:eq:RAND.9), and using the orthonormality of the eigenfunctions , we find that the wave mode amplitudes solve paraxial equations coupled by the random fluctuations in ,
\hb@xt@.01(3.20)  
\hb@xt@.01(3.21) 
We dropped the higherorder terms that do not play a role in the limit , and replaced the equality with the approximate sign. To simplify our notation, we omit henceforth all the arguments in the equations, except those of . The arguments will be spelled out only in definitions.
The coupling matrix in (LABEL:eq:RAND.14)(LABEL:eq:RAND.15) is given by
\hb@xt@.01(3.22) 
It takes this simple diagonal form because we assumed a homogeneous background speed . If we had a variable speed , the matrix would not be diagonal, and the modes with would be coupled. However, the results of the asymptotic analysis below would still hold, because the coupling would become negligible in the limit considered in section LABEL:sect:LR, due to rapid phases arising in the right hand sides of (LABEL:eq:RAND.14), (LABEL:eq:RAND.15).
The equations for the evanescent components are obtained similarly,
\hb@xt@.01(3.23) 
and they are augmented with the decay conditions as , for all .
4 The limit process
We characterize next the wave field in the asymptotic limit . We begin with the paraxial long range scaling that gives significant net scattering, and then take the limit. The scaling has already been described at the beginning of section LABEL:sect:RANDWG.
4.1 Asymptotic scaling
We obtain from (LABEL:eq:RAND.14)(LABEL:eq:RAND.17) that the propagating mode amplitudes satisfy the block diagonal system of partial differential equations
\hb@xt@.01(4.1) 
for . Again, the approximate sign means equal to leading order.
Because the right hand side in (LABEL:eq:LR.1) is small, of order , and has zero statistical expectation, it follows from [8, Chapter 6] that there is no net scattering effect until we reach ranges of order . Thus, we let
\hb@xt@.01(4.2) 
with scaled range independent of . The source directivity in the range direction suggests observing the wavefield on a crossrange scale that is smaller than that in range. We choose it as
\hb@xt@.01(4.3) 
with scaled crossrange independent of , to balance the two terms in the paraxial operators in (LABEL:eq:LR.1).
Our goal is to characterize the limit of the mode amplitudes in the paraxial long range scaling regime (LABEL:eq:LR.2)(LABEL:eq:LR.3). We denote them by
\hb@xt@.01(4.4) 
and obtain from (LABEL:eq:LR.1)(LABEL:eq:LR.3) that they satisfy the scaled system
\hb@xt@.01(4.5) 
for , with initial conditions
\hb@xt@.01(4.6) 
and end conditions
\hb@xt@.01(4.7) 
4.2 The random propagator
Let us rewrite (LABEL:eq:LR.5) in terms of the random propagator matrix , the solution of the initial value problem
\hb@xt@.01(4.8) 
Here is the identity matrix, is the Dirac delta distribution in , and and are matrices with entries given by partial differential operators in , with deterministic coefficients. We can define them from (LABEL:eq:LR.5) once we note that the solution
\hb@xt@.01(4.9) 
follows from
\hb@xt@.01(4.10) 
Here is the vector of backward going amplitudes at the beginning of the randomly perturbed section of the waveguide, and it can be eliminated using the boundary identity
\hb@xt@.01(4.11) 
The initial conditions are given in (LABEL:eq:LR.6).
We obtain from (LABEL:eq:LR.5) that and have the block form
\hb@xt@.01(4.12) 
where the bar denotes complex conjugation. The blocks are diagonal, with entries
\hb@xt@.01(4.13) 
and
\hb@xt@.01(4.14) 
The entries of the diagonal blocks depend only on the mode indices and the frequency, via . The entries of the offdiagonal blocks are rapidly oscillating, due to the large phases proportional to .
The symmetry relations satisfied by the blocks in and imply that the propagator has the form
\hb@xt@.01(4.15) 
with complex, diagonal blocks and .
4.3 The diffusion limit
The limit of as is a multidimensional Markov diffusion process, with entries satisfying a system of ItôSchrödinger equations. This follows from the diffusion approximation theorem [14, 15], see also [8, Chapter 6], applied to system (LABEL:eq:LR.8).
When computing the generator of the limit process, we obtain that due to the fast phases in the offdiagonal blocks of and , the forward and backward going amplitudes decouple as . This implies that there is no backscattered field in the limit, because the backward going amplitudes are set to zero at . Equation (LABEL:eq:LR.10) simplifies as
\hb@xt@.01(4.16) 
where the initial conditions are given in (LABEL:eq:LR.6). We call the complex diagonal matrix
\hb@xt@.01(4.17) 
the transfer process, because it gives the amplitudes of the forward going modes at positive ranges , in terms of the initial conditions at . The limit transfer process is described in the next proposition. It follows straight from [14, 15].
Proposition 4.1
As , converges weakly and in distribution to the diffusion Markov process . This process is complex and diagonal matrix valued, with entries solving the ItôSchrödinger equations
\hb@xt@.01(4.18) 
for , and initial conditions