## 1 Introduction

We analyze the notion of “field-field” cross correlations associated with scattered coda waves, observed at pairs of distinct receivers, to obtain an estimate of the Green's function with an emphasis on high-frequency body wave reflections.

As a model configuration for the crust, we consider a slab in which random medium fluctuations occur. The bottom of the slab is bounded by a deterministic discontinuity (a smooth reflector). We consider waves incident from above the slab, and place our receivers within the slab to study the corresponding cross correlations. Otherwise, in the entire configuration, the medium has a deterministic, smoothly varying component. In this general case, the incident wave is decomposed into wave packets. Each wave packet contains a particular scale. The decomposition is used to select a scale in relation to the fluctuations component of the medium in the slab. Through localization in phase space, the propagation and scattering of an incident wave packet can be described by a coupled system of paraxial wave equations written in curvilinear, boundary normal coordinates relative to the support of the wave packet in the deterministic component of the medium. Thus, in principle, each incident wave packet generates its own system of paraxial wave equations. The accuracy of this description can be proven to improve with increasingly finer scales.

In this article, we assume that the background medium is constant, that the slab is flat, and we consider a single wave packet in the limit of fine scales (high frequency). Indeed, we view the cross correlations in the context of parametrix constructions. The paraxial form of the system that we obtain in the limit allows for the use of Itô’s stochastic calculus (for Hilbert-space valued processes) to analyze the scattering due to the random fluctuations; indeed, it enables the closure of the hierarchy of moment equations.

The solution procedure of the coupled system of paraxial wave equations is based on an invariant embedding type approach, generating a transmission and a reflection operator capturing the scattering due to the random fluctuations in the medium. Thus, we arrive at a coupled system of Riccati equations for the two operator kernels. In the limit of fine scales in the sense of distributions, we then obtain a decoupled system of linear Itô-Schrödinger equations for the limiting transmission and reflection operator kernels, with a real-valued Brownian field; to be precise, only the moments of these kernels converge in the limit.

The solutions to the Itô-Schrödinger equations define the transmitted and backscattered fields, at least their statistics. The transmitted field is partly coherent; the backscattered field is weak and fully incoherent in the absence of a smooth (deterministic) reflector, while it is partly coherent in the presence of a reflector. We analyze the Wigner distributions of these fields. We obtain a description of the cross correlation between two points in terms of a blurring transformation of the paraxial random Green's function between these points. This transformation is statistically stable (in the sense that it does not depend on the realization of the random fluctuations of the medium) and it depends on the statistics of the fluctuations of the random medium.

In the past decade, the understanding of how cross-correlating diffuse fields recaptures the Green's function, has been an important topic of research (van Tiggelen, 2003). Cross correlating (diffuse) coda waves as discussed in Campillo and Paul (2003) resulted in the retrieval of surface waves observed at one station and excited at the other station. The idea of using ambient noise for the retrieval of a body-wave reflection response, in a planarly layered medium, dates back to Claerbout (1968). The retrieval of direct and reflected body waves using teleseismic (S-wave) coda was discussed in Tonegawa et al. (2009). The exploitation of a scattering medium in capturing the Green's function by field-field cross correlations was studied by Derode et al. (2003). However, the mathematical analysis of field-field cross correlations in this setting from the point of view of stochastic calculus has just begun. In our paper, we address a scattering regime in which the field is partly coherent and we aim at retrieving the Green's function for the particular realization of the random medium.

## 2 Scaling and assumptions

We consider acoustic waves propagating in 1+d spatial dimensions. The governing equations are

$$\rho (z,x)\frac{\partial u}{\partial t}+\nabla p=F,\text{\hspace{1em}}\text{\hspace{1em}}\frac{1}{K(z,x)}\frac{\partial p}{\partial t}+\nabla \cdot u=0\text{,}$$ | (1) |

**u**is the velocity field, ρ is the density of mass, and K is the bulk modulus of the medium; $(z,x)\in \mathbb{R}\times {\mathbb{R}}^{d}$ denotes the space coordinates. The source is modeled by the forcing term

**F**. We consider a configuration in which a random slab occupies the region

$${\Omega}_{r}=\{(z,x),x\in {\mathbb{R}}^{d},-L\le z\le 0\}\text{,}$$ |

_{r}, which vary rapidly in space. To simplify the analysis, the deterministic component of the medium in Ω

_{r}is assumed to be constant.

The medium is assumed to be matched at the top interface z = 0 (transparent boundary conditions). However, the interior interface can act as a smooth reflector and is part of the deterministic component of the medium. This model naturally captures primary reflections if the random fluctuations were absent. Indeed, we consider a mismatch at the bottom of the slab, z = –L < 0, which gives jump conditions at the interior interface. The deterministic component of the medium is thus given by

$$\frac{1}{K(z,x)}=\left\{\begin{array}{ll}{K}_{0}^{-1}\hfill & \text{if}\text{\hspace{0.28em}}z\le -L,\hfill \\ {K}_{1}^{-1}\left(1+\nu (z,x)\right)\hfill & \text{if}\text{\hspace{0.28em}}z\in (-L,0),\hfill \\ {K}_{1}^{-1}\hfill & \text{if}\text{\hspace{0.28em}}z\ge 0,\hfill \end{array}\right.$$ |

$$\rho (z,x)=\left\{\begin{array}{ll}{\rho}_{0}\hfill & \text{if}\text{\hspace{0.28em}}z\le -L,\hfill \\ {\rho}_{1}\hfill & \text{if}\text{\hspace{0.28em}}z\in (-L,0),\hfill \\ {\rho}_{1}\hfill & \text{if}\text{\hspace{0.28em}}z\ge 0,\hfill \end{array}\right.$$ |

**x**) models the medium fluctuations, with correlation length l

_{K}. The deterministic wavespeed for z > –L is ${c}_{1}=\sqrt{{K}_{1}/{\rho}_{1}}$, and for z < –L is ${c}_{0}=\sqrt{{K}_{0}/{\rho}_{0}}$.

The source, **F**, is located at z = z_{s} ≥ 0; we assume, here, that z_{s} = 0 and write $F(t,z,x)={f}_{s}(t,x)\delta (z-{z}_{s}){e}_{z}$, where **e**_{z} denotes the unit vector pointing in the z-direction, signifying a body force. We shall refer to waves propagating in the positive z direction as upgoing. The source generates downgoing waves which “propagate” through the random medium, are reflected by the interface at z = –L, and “propagate” up through the medium.

The source is now assumed to be of the form:

$$F(t,z,x)={\chi}_{\epsilon}(t,x)\delta (z){e}_{z}$$ | (2) |

The forcing function ${\chi}_{\epsilon}(t,x)=\chi ({\epsilon}^{-4}t,{\epsilon}^{-2}x)$ gives a wave packet oriented in the z direction of scale k, $\epsilon ={2}^{-k/4}$, corresponding with a dyadic parabolic decomposition of phase space. Here $\stackrel{\u02c6}{\chi}$ denotes a window function supported in a box in the $\left(\omega ,\text{\hspace{0.17em}}\kappa \right)$ Fourier domain (the support of $\stackrel{\u02c6}{\chi}$ is a finite distance away from the ω = 0 axis) and we use the Fourier transforms:

$$\chi (t,x)=\frac{1}{{(2\pi )}^{d+1}}\iint \stackrel{\u02c6}{\chi}(\omega ,\kappa ){e}^{-i(\omega t-k\cdot x)}\text{d}\omega \text{d}k=\frac{1}{2\pi}\int \stackrel{\u2323}{\chi}(\omega ,k){e}^{-i\omega t}\text{d}k$$ |

Thus, the transverse width, R_{0}, of the source function is of order ɛ^{2}. This scaling is consistent with the paraxial regime, which will be discussed in the next section, in as much as that it generates the wave solution up to an error of order ɛ^{2}. We assume that:

- • the correlation length or radius of the fluctuations, l
_{K}, is of the same order as R_{0}; this regime guarantees non-trivial interaction between the fluctuations of the medium and the waves; - • the propagation distance is of the order of L, which is of order 1; the ratio between the propagation distance and the correlation length of the fluctuations is of order ɛ
^{−2}.

The wave packet has a central frequency of order 2^{k}. We will establish results, which converge as k → ∞, signifying the high-frequency regime.

Henceforth we shall assume non-dimensionalized units chosen such that the deterministic bulk modulus K_{1} and density ρ_{1} in the region Ω_{r} are one and, hence, the deterministic wavespeed ${c}_{1}=\sqrt{{K}_{1}/{\rho}_{1}}$ and impedance ${Z}_{1}=\sqrt{{K}_{1}{\rho}_{1}}$ are equal to one. We assume that ${Z}_{0}=\sqrt{{K}_{0}{\rho}_{0}}\ne 1$. The medium fluctuations attain the scaled form

$$\frac{1}{K(z,x)}=\left\{\begin{array}{ll}{K}_{0}^{-1}\hfill & \text{if}\text{\hspace{0.28em}}z\le -L,\hfill \\ 1+{\epsilon}^{3}\nu (\frac{z}{{\epsilon}^{2}},\frac{x}{{\epsilon}^{2}})\hfill & \text{if}\text{\hspace{0.28em}}z\in (-L,0),\hfill \\ 1\hfill & \text{if}\text{\hspace{0.28em}}z\ge 0,\hfill \end{array}\right.$$ |

$$\rho (z,x)=\left\{\begin{array}{ll}{\rho}_{0}\hfill & \text{if}\text{\hspace{0.28em}}z\le -L,\hfill \\ 1\hfill & \text{if}\text{\hspace{0.28em}}z\in (-L,0),\hfill \\ 1\hfill & \text{if}\text{\hspace{0.28em}}z\ge 0,\hfill \end{array}\right.$$ |

$$C(z,x)=E[\nu ({z}^{\prime}+z,x\text{'}+x)\nu ({z}^{\prime},x\text{'})]\text{,}$$ | (3) |

$$D(x)={\int}_{-\infty}^{\infty}C(z,x)\text{d}z$$ | (4) |

We assume that v satisfies strong mixing conditions in z. The amplitude ɛ^{3} of the fluctuations has been chosen so as to obtain an effective limit of order one when ɛ goes to zero. That is, if the magnitude of the fluctuations is smaller than ɛ^{3}, then the wave would propagate as if the medium were homogeneous, while if the order of magnitude is larger, then the wave would not penetrate the slab at all. The scaling that we consider here, controlled by the choice of wave packet, corresponds to the partly coherent regime.

## 3 Coupled system of paraxial equations, and transmission and reflection operators

We eliminate the transverse components of the velocity field. Because both the medium and the source have transverse spatial variations at the scale ɛ^{2}, it is convenient to rescale the transverse coordinates accordingly, that is, ɛ^{2}x → x. In the region Ω_{r}, we then introduce the directional decomposition in the deterministic medium component,

$$p(t,z,{\epsilon}^{2}x)=\frac{1}{2\pi}\int \left({\stackrel{\u2323}{a}}^{\epsilon}(\omega ,z,x){e}^{i\omega \frac{z}{{\epsilon}^{4}}}+{\stackrel{\u2323}{b}}^{\epsilon}(\omega ,z,x){e}^{-i\omega \frac{z}{{\epsilon}^{4}}}\right){e}^{-i\omega \frac{t}{{\epsilon}^{4}}}\text{d}\omega \text{,}$$ | (5) |

The wave-packet source enables the use of the paraxial approximation in the deterministic component of the medium with the above-mentioned estimate, leading to the system for z ∈ (–L,0):

$$\frac{\partial {\stackrel{\u2323}{a}}^{\epsilon}}{\partial z}=\left[\frac{i\omega}{2\epsilon}\nu \left(\frac{z}{{\epsilon}^{2}},x\right)+\frac{i}{2\omega}{\Delta}_{x}\right]{\stackrel{\u2323}{a}}^{\epsilon}+{e}^{-2i\omega \frac{z}{{\epsilon}^{4}}}\left[\frac{i\omega}{2\epsilon}\nu \left(\frac{z}{{\epsilon}^{2}},x\right)+\frac{i}{2\omega}{\Delta}_{x}\right]{\stackrel{\u2323}{b}}^{\epsilon}\text{,}$$ | (6) |

$$\frac{\partial {\stackrel{\u2323}{b}}^{\epsilon}}{\partial z}=-{e}^{2i\omega \frac{z}{{\epsilon}^{4}}}\left[\frac{i\omega}{2\epsilon}\nu \left(\frac{z}{{\epsilon}^{2}},x\right)+\frac{i}{2\omega}{\Delta}_{x}\right]{\stackrel{\u2323}{a}}^{\epsilon}-\left[\frac{i\omega}{2\epsilon}\nu \left(\frac{z}{{\epsilon}^{2}},x\right)+\frac{i}{2\omega}{\Delta}_{x}\right]{\stackrel{\u2323}{b}}^{\epsilon}\text{,}$$ | (7) |

$${\stackrel{\u2323}{b}}^{\epsilon}(\omega ,z={0}^{-},x)=-\frac{1}{2}\stackrel{\u2323}{\chi}(\omega ,x)\text{,}$$ | (8) |

$${\stackrel{\u2323}{a}}^{\epsilon}(\omega ,z={(-L)}^{+},x)={R}_{0}{e}^{2i\omega \frac{L}{{\epsilon}^{4}}}{\stackrel{\u2323}{b}}^{\epsilon}(\omega ,z={(-L)}^{+},x)\text{,}$$ | (9) |

$${\stackrel{\u2323}{b}}^{\epsilon}(\omega ,{(-L)}^{-},x)={T}_{0}{\stackrel{\u2323}{b}}^{\epsilon}(\omega ,{(-L)}^{+},x)\text{,}$$ |

$${\stackrel{\u2323}{a}}^{\epsilon}(\omega ,{0}^{-},x)={\stackrel{\u2323}{a}}^{\epsilon}(\omega ,{0}^{-},x)+\frac{1}{2}\stackrel{\u2323}{\chi}(\omega ,x)$$ |

In the deterministic case, with v = 0, it can be proven that this system decouples with solutions accurate up to order ɛ^{2}.

We invoke an invariant imbedding approach to obtain the representation valid for –L < z < 0:

$${\stackrel{\u2323}{b}}^{\epsilon}(\omega ,{(-L)}^{-},x)={T}_{0}\int \text{\hspace{0.28em}}\stackrel{\u2323}{T}{}^{\epsilon}(\omega ,-L,z,x,x\text{'}){\stackrel{\u2323}{b}}^{\epsilon}(\omega ,z,x\text{'})\text{d}x\text{'}\text{,}$$ | (10) |

$${\stackrel{\u2323}{a}}^{\epsilon}(\omega ,z,x)={R}_{0}{e}^{2i\omega \frac{L}{{\epsilon}^{4}}}\int \text{\hspace{0.28em}}\stackrel{\u2323}{R}{}^{\epsilon}(\omega ,-L,z,x,x\text{'}){\stackrel{\u2323}{b}}^{\epsilon}(\omega ,z,x\text{'})\text{d}x\text{'}\text{,}$$ | (11) |

The formulation generalizes to curvilinear coordinates if the deterministic wave speed in the slab is no longer constant, and the rays are no longer straight, using techniques from microlocal analysis and harmonic analysis.

## 4 Itô-Schrödinger diffusion models for transmitted and backscattered fields

We centre according to the travel time associated with the deterministic medium component (constant wave speed, here) and define the transmitted and reflected pressure fields by

$${p}_{R}^{\epsilon}(s,x)\u2254p(2L+{\epsilon}^{4}s,{0}^{+},{\epsilon}^{2}x)-\frac{1}{2}\chi (\frac{2L}{{\epsilon}^{4}}+s,x)\text{,}$$ | (12) |

$${p}_{T}^{\epsilon}(s,x)\u2254p(L+{\epsilon}^{4}s,{(-L)}^{-},{\epsilon}^{2}x)$$ | (13) |

The field ${p}_{T}^{\epsilon}(s,x)$ is the field observed just below the bottom interface at z = (–L)^{−}; the field ${p}_{R}^{\epsilon}(s,x)$ is the field observed just above the top interface at z = 0^{+}:

$${p}_{R}^{\epsilon}(s,x)=\frac{1}{2\pi}\int {\stackrel{\u2323}{a}}^{\epsilon}(\omega ,{0}^{+},x){e}^{-i\omega s}\text{d}\omega \text{,}$$ | (14) |

$${p}_{T}^{\epsilon}(s,x)=\frac{1}{2\pi}\int {\stackrel{\u2323}{b}}^{\epsilon}(\omega ,{(-L)}^{-},x){e}^{-i\omega s}\text{d}\omega $$ | (15) |

These fields are now characterized via effective scaling limit models for the transmission and reflection operators:

**Proposition 4.1** The processes ${({p}_{T}^{\epsilon}(s,x))}_{s\in \mathbb{R},x\in {\mathbb{R}}^{d}}$, ${({p}_{R}^{\epsilon}(s,x))}_{s\in \mathbb{R},x\in {\mathbb{R}}^{d}}$ converge in distribution as ɛ → 0 in the space ${C}^{0}(\mathbb{R},{L}^{2}({\mathbb{R}}^{d},{\mathbb{R}}^{2}))\cap {L}^{2}(\mathbb{R},{L}^{2}({\mathbb{R}}^{d},{\mathbb{R}}^{2}))$ to the limit processes ${({p}_{T}(s,x))}_{s\in \mathbb{R},x\in {\mathbb{R}}^{d}}$, ${({p}_{R}(s,x))}_{s\in \mathbb{R},x\in {\mathbb{R}}^{d}}$ given by

$${p}_{R}(s,x)=-\frac{{R}_{0}}{4\pi}\iint \text{\hspace{0.28em}}\stackrel{\u2323}{R}(\omega ,-L,0,x,x\text{'})\stackrel{\u2323}{\chi}(\omega ,x\text{'})\text{d}x\text{'}{e}^{-i\omega s}\text{d}\omega \text{,}$$ | (16) |

$${p}_{T}(s,x)=-\frac{{T}_{0}}{4\pi}\iint \text{\hspace{0.28em}}\stackrel{\u2323}{T}(\omega ,-L,0,x,x\text{'})\stackrel{\u2323}{\chi}(\omega ,x\text{'})\text{d}x\text{'}{e}^{-i\omega s}\text{d}\omega \text{,}$$ | (17) |

$${L}^{2}(\mathbb{R},{L}^{2}({\mathbb{R}}^{d},{\mathbb{R}}^{2}))={L}^{2}(\mathbb{R}\times {\mathbb{R}}^{d},{\mathbb{R}}^{2})\text{.}$$ |

The operators

$${(\stackrel{\u2323}{R}(\omega ,-L,z,x,x\text{'}))}_{z\in [-L,0]}$$ |

$${(T(\omega ,-L,z,x,x\text{'}))}_{z\in [-L,0]}$$ |

$$\text{d}\stackrel{\u2323}{R}(\omega ,-L,z,x,x\text{'})=\frac{i}{2\omega}\left({\Delta}_{x}+{\Delta}_{x\text{'}}\right)\stackrel{\u2323}{R}(\omega ,-L,z,x,x\text{'})\text{d}z+\frac{i\omega}{2}\stackrel{\u2323}{R}(\omega ,-L,z,x,x\text{'})\circ \left(\text{d}B(z,x)+\text{d}B(z,x\text{'})\right)\text{,}$$ | (18) |

$$\text{d}\stackrel{\u2323}{T}(\omega ,-L,z,x,x\text{'})=\frac{i}{2\omega}{\Delta}_{x\text{'}}\stackrel{\u2323}{T}(\omega ,-L,z,x,x\text{'})\text{d}z+\frac{i\omega}{2}\stackrel{\u2323}{T}(\omega ,-L,z,x,x\text{'})\circ \text{d}B(z,x\text{'})\text{,}$$ | (19) |

$$\stackrel{\u2323}{R}(\omega ,-L,z=-L,x,{x}^{\prime})=\delta (x-{x}^{\prime}),\text{\hspace{1em}}\stackrel{\u2323}{T}(\omega ,-L,z=-L,x,{x}^{\prime})=\delta (x-{x}^{\prime})\text{.}$$ |

The symbol ° stands for the Stratonovich stochastic integral, and B(z, **x**) is a real-valued Brownian field with covariance

$$E[B({z}_{1},{x}_{1})B({z}_{2},{x}_{2})]=\mathrm{min}\{{z}_{1},{z}_{2}\}{C}_{0}({x}_{1}-{x}_{2})$$ | (20) |

Making use of the semigroup property of the effective operators we also find that the joint law for the direct arrival to the two points of observation located at (–L_{1}, ɛ^{2}**x**_{1}) and (–L_{2} ɛ^{2}**x**_{2}), with –L ≤ –L_{2} ≤ –L_{1} ≤ 0, can be characterized by

$$p({L}_{1}+{\epsilon}^{4}s,-{L}_{1},{\epsilon}^{2}{x}_{1})\sim -\frac{1}{4\pi}\iint \text{\hspace{0.28em}}\stackrel{\u2323}{T}(\omega ,-{L}_{1},0,{x}_{1},x\text{'})\stackrel{\u2323}{\chi}(\omega ,x\text{'})\text{d}x\text{'}{e}^{-i\omega s}\text{d}\omega \text{,}$$ | (21) |

$$p({L}_{2}+{\epsilon}^{4}s,-{L}_{2},{\epsilon}^{2}{x}_{2})\sim -\frac{1}{4\pi}\iint \text{\hspace{0.28em}}\stackrel{\u2323}{T}(\omega ,-{L}_{2},-{L}_{1},{x}_{2},x\mathrm{\text{'}\text{'}})\text{\hspace{0.28em}}\stackrel{\u2323}{T}(\omega ,-{L}_{1},0,x\mathrm{\text{'}\text{'}},x\text{'})\stackrel{\u2323}{\chi}(\omega ,x\text{'})\text{d}x\text{'}\text{d}x\mathrm{\text{'}\text{'}}{e}^{-i\omega s}\text{d}\omega $$ | (22) |

_{2}≤ –L

_{1}≤ 0, the operators $\text{\hspace{0.28em}}\stackrel{\u2323}{T}(\omega ,-L,-{L}_{2},x,x\text{'})$ and $\text{\hspace{0.28em}}\stackrel{\u2323}{R}(\omega ,-L,-{L}_{2},x,x\text{'})$ are statistically independent of $\text{\hspace{0.28em}}\stackrel{\u2323}{T}(\omega ,-{L}_{2},-{L}_{1},x,x\text{'})$ and $\text{\hspace{0.28em}}\stackrel{\u2323}{R}(\omega ,-{L}_{2},-{L}_{1},x,x\text{'})$, which we exploit in evaluating the cross-correlations. We remark here also that the transmission operator over the sub-slab (–L

_{1},0) appears in both expressions in (21)–(22) and it is this pairing that will lead to an expression for the cross correlations in terms of a statistically stable filter or transformation below. The general statistical properties of the operators $\stackrel{\u2323}{R}$, $\stackrel{\u2323}{T}$ were studied in Garnier and Sølna (2009a).

## 5 Characterization of cross correlations

Here, we characterize the “field-field” correlation function between the points (–L_{1}, ɛ^{2}**x**_{1}) and (–L_{2}, ɛ^{2}**x**_{2}). We assume here that $T\gg L$ and 0 < L_{1} < L_{2} < L. The field-field correlation function is given by

$${V}_{T}^{\epsilon}(\tau )={\int}_{0}^{T}p(t,-{L}_{1},{\epsilon}^{2}{x}_{1})p(t+\tau ,-{L}_{2},{\epsilon}^{2}{x}_{2})\text{d}t\text{.}$$ | (23) |

The configuration is illustrated in Fig. 1, we compute the cross correlation in between the two points at depth –L_{1} and –L_{2}. We will see that the wave field correlation function concentrates around specific time lags τ that correspond to travel times between the two observation points, and that the time extent of correlation function around these time lags is much smaller, of order ɛ^{4}, i.e., of the same order as the source pulse width.

Under the scaling assumptions of Section 2, we find that the correlations in (23) has leading contributions at four particular time lags:

$${V}_{T}^{\epsilon}(\pm ({L}_{2}-{L}_{1})+{\epsilon}^{4}s)/{\epsilon}^{4}\sim {V}_{t}^{\pm}(s)\text{,}$$ | (24) |

$${V}_{T}^{\epsilon}(\pm (2L-{L}_{1}-{L}_{2})+{\epsilon}^{4}s)/{\epsilon}^{4}\sim {V}_{r}^{\pm}(s)$$ | (25) |

Here, the amplitude scaling ɛ^{−4} corresponds to a re-scaling of the source time traces so that they have order one energy, but plays no significant role as the problem is linear. We shall here focus on the contribution ${V}_{t}^{+}$, which corresponds to correlation of wave components directly transmitted in between the points of observation (Fig. 2). The contribution ${V}_{t}^{+}$ is concentrated around time lag equal to +(L_{2}–L_{1}) and it comes from the correlation between the waves that propagate from the surface to the depth –L_{1} and then to the depth –L_{2}. The contribution ${V}_{t}^{-}$ is concentrated around time lag equal to –(L_{2}–L_{1}) and it comes from the correlation between the waves that have been reflected by the interface at z = –L and that propagate from this interface to the depth –L_{2} and then to the depth –L_{1}. We stress here that these wave components have been strongly affected by the multiple scattering in the medium and understanding how this affects the relation of the components V and Green's functions is our main objective. The contributions ${V}_{r}^{\pm}$ correspond to cross terms and they will be treated in detail elsewhere. We remark here that they give contributions concentrated around times lags that are larger than the ones of the first two contributions (since 2L–L_{1}–L_{2} > L_{2}–L_{1}) and they can be used for estimation of the depth L of the bottom interface.

The correlation component ${V}_{t}^{+}$ can be characterized in distribution in the scaling limit considered in this paper by the following expression using (21), (22) and the definition (23):

$${V}_{t}^{+}(s)=\frac{1}{2\pi}\iint {\Lambda}_{t}^{+}(\omega ,x;{x}_{1},{L}_{1})\text{\hspace{0.28em}}\times \text{\hspace{0.28em}}\stackrel{\u2323}{T}(\omega ,-{L}_{2},-{L}_{1},{x}_{2},{x}_{1}-x)\text{d}x{e}^{-i\omega s}\text{d}\omega \text{,}$$ | (26) |

$${\Lambda}_{t}^{+}(\omega ,x;{x}_{1},{L}_{1})=\frac{1}{4}\iint \text{\hspace{0.28em}}\stackrel{\u2323}{T}(\omega ,-{L}_{1},0,{x}_{1}-x,{y}_{2})\overline{\text{\hspace{0.28em}}\stackrel{\u2323}{T}(\omega ,-{L}_{1},0,{x}_{1},{y}_{1})}\stackrel{\u2323}{\chi}(\omega ,{y}_{2})\overline{\stackrel{\u2323}{\chi}(\omega ,{y}_{1})}\text{d}{y}_{1}\text{d}{y}_{2}$$ | (27) |

In the regime considered in this article, the filter ${\Lambda}_{t}^{+}$ defining the relevant transformation is self-averaging, see de Hoop and Sølna (2009), in the sense that

$${\Lambda}_{t}^{+}(\omega ,x;{x}_{1},{L}_{1})=\frac{1}{4}\iint E[\stackrel{\u2323}{T}(\omega ,-{L}_{1},0,{x}_{1}-x,{y}_{2})\overline{\text{\hspace{0.28em}}\stackrel{\u2323}{T}(\omega ,-{L}_{1},0,{x}_{1},{y}_{1})}]\stackrel{\u2323}{\chi}(\omega ,{y}_{2})\overline{\stackrel{\u2323}{\chi}(\omega ,{y}_{1})}\text{d}{y}_{1}\text{d}{y}_{2}$$ |

In order to characterize the filter ${\Lambda}_{t}^{+}$ we introduce the Wigner transform defined by

$${W}_{\omega}^{T}({L}_{1},x,x\text{'},\kappa ,\kappa \text{'})=\iint \text{\hspace{0.28em}}E[\text{\hspace{0.28em}}\stackrel{\u2323}{T}(\omega ,-{L}_{1},0,x+\frac{y}{2},\mathrm{x\text{'}}+\frac{y\text{'}}{2})\overline{\times \stackrel{\u2323}{T}(\omega ,-{L}_{1},0,x-\frac{y}{2},x\text{'}-\frac{y\text{'}}{2})}]{e}^{-i(\kappa \cdot y+\kappa \text{'}\cdot y\text{'})}\text{d}y\text{d}y\text{'}$$ |

The Wigner transforms can be shown to satisfy a set of transport equations that can be integrated, as shown in Garnier and Sølna (2009b), and we find the following integral representation for ${W}_{\omega}^{T}$:

$${W}_{\omega}^{T}(z,x,\mathrm{x\text{'}},\kappa ,\kappa \text{'})=\frac{1}{{(2\pi )}^{d}}\iint {e}^{-i(\kappa \text{'}+\kappa )\cdot a-i(\mathrm{x\text{'}}-x+\frac{\kappa}{\omega}z)\cdot b}{e}^{\frac{{\omega}^{2}}{4}{\int}_{0}^{z}D(a+\frac{b}{\omega}{z}^{\prime})-D(0)\text{d}{z}^{\prime}}\text{d}a\text{d}b$$ | (28) |

This gives then an integral expression for the filter. In order to get an explicit form for the filter and characterize the associated resolution scale, or filter support scale, we consider next a particular regime of relatively strong scattering or clutter.

### 5.1 The strongly cluttered regime

We introduce the correlation length $\ell $ and the standard deviation σ of the random fluctuations:

$$C(z,x)={\sigma}^{2}{C}_{0}(\frac{z}{\ell},\frac{x}{\ell})$$ |

With this representation we have:

$$D(x)={\sigma}^{2}\ell {D}_{0}(\frac{x}{\ell}),\text{\hspace{1em}}\text{\hspace{1em}}{D}_{0}(x)={\int}_{-\infty}^{\infty}{C}_{0}(z,x)\text{d}z$$ |

We assume also that the normalized autocorrelation function D_{0}(**x**) is at least twice differentiable at **x** = **0**, which corresponds to a smooth random medium.

We introduce the parameter, depending on ω,

$${\beta}_{L}(\omega )=L\frac{{\sigma}^{2}\ell {\omega}^{2}}{4}\text{,}$$ | (29) |

_{L}(ω) being large.

Second, we introduce the frequency-dependent resolution parameter:

$${\Sigma}_{\omega}^{2}({L}_{1})=\frac{16}{{\omega}^{2}\gamma {L}_{1}}\text{,}$$ | (30) |

$$\gamma =-\frac{{\sigma}^{2}}{d\ell}\Delta {D}_{0}(0)$$ |

Using the results of Garnier and Sølna (2009a), we find the following expression for the filter

$${\Lambda}_{t}^{+}(\omega ,x;{x}_{1},{L}_{1})=\frac{1}{4}{(\frac{6}{\pi \gamma {L}_{1}^{3}})}^{d/2}\mathrm{exp}(-\frac{|x{|}^{2}}{2{\Sigma}_{\omega}^{2}({L}_{1})})\int |\stackrel{\u2323}{\chi}(\omega ,{y}_{1}){|}^{2}\text{d}{y}_{1}$$ | (31) |

Thus, we have a situation in which a sharp filter and Green's function estimation have been enabled by the medium correlations. We remark that the expression (31) is valid provided $\left|{x}_{1}\right|\ll \sqrt{\gamma {L}_{1}^{3}}$, which means that the point (–L_{1}, ɛ^{2}**x**_{1}) is in the forward cone of wave energy (Garnier and Sølna, 2009a; de Hoop and Sølna, 2009). Note that this cone is wider than the usual deterministic light cone due to multiple scattering.

We have similarly

$${V}_{t}^{-}(s)=\frac{{R}_{0}^{2}}{2\pi}\iint {\Lambda}_{t}^{-}(\omega ,x;{x}_{2},{L}_{2},L)\text{\hspace{0.28em}}\times \text{\hspace{0.28em}}\stackrel{\u2323}{T}(\omega ,-{L}_{1},-{L}_{2},{x}_{1},{x}_{2}-x)\text{d}x{e}^{-i\omega s}\text{d}\omega \text{,}$$ |

$${\Lambda}_{t}^{-}(\omega ,x;{x}_{2},{L}_{2},L)=\frac{1}{4}\iint \text{\hspace{0.28em}}\stackrel{\u2323}{T}(\omega ,-{L}_{2},0,{y}_{4},{y}_{2})\text{\hspace{0.28em}}\times \text{\hspace{0.28em}}\overline{\stackrel{\u2323}{T}(\omega ,-{L}_{2},0,{y}_{3},{y}_{1})}$$ |

$$\text{\hspace{1em}}\text{\hspace{1em}}\times \text{\hspace{0.28em}}\stackrel{\u2323}{R}(\omega ,-L,-{L}_{2},{x}_{2}-x,{y}_{4})\overline{\text{\hspace{0.28em}}\stackrel{\u2323}{R}(\omega ,-L,-{L}_{2},{x}_{2},{y}_{3})}\stackrel{\u2323}{\chi}(\omega ,{y}_{2})\overline{\stackrel{\u2323}{\chi}(\omega ,{y}_{1})}\text{d}{y}_{1}\text{d}{y}_{2}$$ |

Note that, by reciprocity, we have

$${V}_{t}^{-}(s)=\frac{{\mathcal{R}}_{0}^{2}}{2\pi}\iint {\Lambda}_{t}^{-}(\omega ,x;{x}_{2},{L}_{2},L)\overline{\text{\hspace{0.28em}}\stackrel{\u2323}{T}(\omega ,-{L}_{2},-{L}_{1},{x}_{2}-x,{x}_{1})}\text{d}x{e}^{-i\omega s}\text{d}\omega \text{,}$$ | (32) |

We can also get a simple expression of the filter in the strongly cluttered regime. Self-averaging implies that

$${\Lambda}_{t}^{-}(\omega ,x;{x}_{2},{L}_{2},L)=\frac{1}{4}\iint E[\text{\hspace{0.28em}}\stackrel{\u2323}{T}(\omega ,-{L}_{2},0,{y}_{4},{y}_{2})\overline{\text{\hspace{0.28em}}\stackrel{\u2323}{T}(\omega ,-{L}_{2},0,{y}_{3},{y}_{1})}]$$ | (33) |

$$\text{\hspace{1em}}\times E[\text{\hspace{0.28em}}\stackrel{\u2323}{R}(\omega ,-L,-{L}_{2},{x}_{2}-x,{y}_{4})\overline{\text{\hspace{0.28em}}\stackrel{\u2323}{R}(\omega ,-L,-{L}_{2},{x}_{2},{y}_{3})}]\stackrel{\u2323}{\chi}(\omega ,{y}_{2})\overline{\stackrel{\u2323}{\chi}(\omega ,{y}_{1})}\text{d}{y}_{1}\text{d}{y}_{2}\text{d}{y}_{3}\text{d}{y}_{4}$$ |

## 6 Discussion

We presented an analysis for partly coherent and partly incoherent body waves generated by a (teleseismic) wave packet remotely incident on a slab (the crust) containing a medium consisting of a deterministic component and a random field, based on the dyadic parabolic decomposition of phase space coupled to the scaling of the random fluctuations. The deterministic component consists, here, of a planarly layered medium, but can be generalized to contain conormal singularities (discontinuities) combined with smooth wave speed variations. The wave packets provide a frame to represent and study partly coherent wave fields.

To obtain information about the deterministic medium component, one needs to consider “field-field” cross correlations. We showed that these cross correlations are characterized by a transformation (blurring filter) of the random paraxial Green's function, between the points at which the fields are taken, the transformation containing information about the statistics of the random fluctuations. If the points are taken purely transverse to the propagation direction of the wave packet in the deterministic component, the blurring significantly increases, which is consistent with the usual stationary phase arguments. In general our results fully characterize the cross correlation function in the regime in which the source and random medium scaling is such that the waves remain concentrated in phase space.