Paraxial coupling of propagating modes in three-dimensional waveguides with random boundaries

Paraxial coupling of propagating modes in three-dimensional waveguides with random boundaries

Liliana Borcea111Computational and Applied Mathematics, Rice University, Houston, TX 77005. and Josselin Garnier222Laboratoire de Probabilités et Modèles Aléatoires & Laboratoire Jacques-Louis Lions, Université Paris VII, Site Chevaleret, 75205 Paris Cedex 13, France.

We analyze long range wave propagation in three-dimensional 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 mode-dependent scales: the scattering mean free path, the cross-range 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 three-dimensional 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 cross-range 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


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,


The bottom of the waveguide is assumed rigid


and we take a pressure release boundary condition at the perturbed top boundary


Perturbed means that the boundary has small fluctuations around the mean depth ,


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 two-dimensional 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.

Figure 1.1: Schematic of the problem setup. The system of coordinates has range origin at the source. The rigid bottom boundary is assumed flat and the pressure release top boundary has fluctuations around the value . The cross-range and the range are unbounded, that is .

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 mode-dependent scattering mean free path, which is the distance over which the modes lose coherence; the mode-dependent decoherence length, which is the cross-range offset over which the mode amplitudes decorrelate; and the mode-dependent 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


with Fourier coefficients satisfying a separable problem for the Helmholtz equation


with boundary conditions


and outgoing radiation conditions at . The Fourier transform of the source


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,


The decomposition is in the orthonormal basis of the eigenfunctions of the self-adjoint differential operator in ,

with eigenvalues that are simple [16].

To simplify the analysis, we assume in this paper that the wave speed is homogeneous


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


The eigenvalues are


where is the wavenumber, and only the first of them are non-negative


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 two-dimensional Helmholtz equation


with outgoing, radiation conditions at . The evanescent components solve


with decay condition at Here we introduced the coefficients of the source profile in the basis of the eigenfunctions


and the mode wavenumbers


We assume that none of the vanishes in the bandwidth, so that there are no standing waves. That is to say,


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 cross-range profile of the source is larger than the wavelength. The source generates a quasi-plane wave, with slowly varying envelope satisfying a Schrödinger-like equation.

Explicitly, we assume that the source is of the form


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 cross-range scales, similar to , and at range scale, similar to the Rayleigh length. We rename the field in this scaling as


The Fourier coefficients of (LABEL:eq:poeps) are given by the scaled version of (LABEL:eq:HOM.4)


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



for . The evanescent components are obtained similarly from (LABEL:eq:HOM.11)



for . These modes are exponentially damped and can be neglected.

In summary, the paraxial approximation of the wave field is given by


It is a superposition of forward going modes with complex valued amplitudes given by (LABEL:eq:HOM.15), and solving the paraxial equations


with initial conditions


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:

  1. The transverse width of the source and the central wavelength satisfy


    as in the previous section. The correlation length of the boundary fluctuations is similar to ,


    so that there is a non-trivial interaction between the boundary fluctuations and the wavefield. Here is the scaled order-one correlation length defined below.

  2. The scale of the propagation distance is much larger than . More precisely,


    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 high-frequency scaling assumption (LABEL:eq:sc.3) ensures that the propagation distance is similar to the Rayleigh length.

  3. 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


The process is bounded, zero-mean, stationary and mixing, meaning in particular that its covariance is integrable111More 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


and we denote by its integral over ,


Our assumption on the differentiability of implies that is four times differentiable. Note that is the maximum of the integrated covariance , so we have


We define the scaled square amplitude and correlation length of the boundary fluctuations through the equations


3.1 Change of coordinates

We introduce the change of coordinates from to , with


It straightens the boundary to , for any and . The pressure field in the new coordinates is denoted by


It satisfies the simple boundary conditions


and the partial differential equation


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


The higher-order 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)


We define the forward and backward going wave mode amplitudes and by


so that the complex valued amplitudes of the propagating modes can be written as

Definition (LABEL:eq:def.ab) implies that


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


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


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 ,


We dropped the higher-order 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


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,


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


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


with scaled range independent of . The source directivity in the range direction suggests observing the wavefield on a cross-range scale that is smaller than that in range. We choose it as


with scaled cross-range 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


and obtain from (LABEL:eq:LR.1)-(LABEL:eq:LR.3) that they satisfy the scaled system


for , with initial conditions


and end conditions


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


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


follows from


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


The initial conditions are given in (LABEL:eq:LR.6).

We obtain from (LABEL:eq:LR.5) that and have the block form


where the bar denotes complex conjugation. The blocks are diagonal, with entries




The entries of the diagonal blocks depend only on the mode indices and the frequency, via . The entries of the off-diagonal 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


with complex, diagonal blocks and .

4.3 The diffusion limit

The limit of as is a multi-dimensional 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 off-diagonal 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


where the initial conditions are given in (LABEL:eq:LR.6). We call the complex diagonal matrix


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


for , and initial conditions