Fourth-moment Analysis for Wave Propagation in the White-Noise Paraxial Regime
In this paper we consider the Itô-Schrödinger model for wave propagation in random media in the paraxial regime. We solve the equation for the fourth-order moment of the field in the regime where the correlation length of the medium is smaller than the initial beam width. As applications we prove that the centered fourth-order moments of the field satisfy the Gaussian summation rule, we derive the covariance function of the intensity of the transmitted beam, and the variance of the smoothed Wigner transform of the transmitted field. The second application is used to explicitly quantify the scintillation of the transmitted beam and the third application to quantify the statistical stability of the Wigner transform.
AMS subject classifications. 60H15, 35R60, 74J20.
Key words. Waves in random media, parabolic approximation, scintillation, Wigner transform.
In many wave propagation scenarios the medium is not constant, but varies in a complicated fashion on a scale that may be small compared to the total propagation distance. This is the case for wave propagation through the turbulent atmosphere, the earth’s crust, the ocean, and complex biological tissue for instance. If one aims to use transmitted or reflected waves for communication or imaging purposes it is important to characterize how such microstructure affects and corrupts the wave. Such a characterization is particularly important for modern imaging techniques such as seismic interferometry or coherent interferometric imaging that correlate wave field traces that have been strongly corrupted by the microstructure and use their space-time correlation function for imaging. The wave field correlations are second-order moments of the wave field and a characterization of the signal-to-noise ratio then involves a fourth-order moment calculation.
Motivated by the situation described above we consider wave propagation through time-independent media with a complex spatially varying index of refraction that can be modeled as the realization of a random process. Typically we cannot expect to know the index of refraction pointwise, but we may be able to characterize its statistics and we are interested in how the statistics of the medium affects the statistics of the wave field. In its most common form, the analysis of wave propagation in random media consists in studying the field solution of the scalar time-harmonic wave or Helmholtz equation
where is the free space homogeneous wavenumber and is a randomly heterogeneous index of refraction. Since the index of refraction is a random process, the field is also a random process whose statistical behavior can be characterized by the calculations of its moments. Even though the scalar wave equation is simple and linear, the relation between the statistics of the index of refraction and the statistics of the field is highly nontrivial and nonlinear. In this paper we consider a primary scaling regime corresponding to long-range beam propagation and small-scale medium fluctuations giving negligible backscattering. This is the so-called white-noise paraxial regime, as described by the Itô-Schrödinger model, which is presented in Section 2. This model is a simplification of the model (1.1) since it corresponds to an evolution problem, but yet in the regime that we consider it describes the propagated field in a weak sense in that it gives the correct statistical structure of the wave field. The Itô-Schrödinger model can be derived rigorously from (1.1) by a separation of scales technique in the high-frequency regime (see  in the case of a randomly layered medium and [23, 24, 25] in the case of a three-dimensional random medium). It models many situations, for instance laser beam propagation , time reversal in random media [5, 38], underwater acoustics , or migration problems in geophysics . The Itô-Schrödinger model allows for the use of Itô’s stochastic calculus, which in turn enables the closure of the hierarchy of moment equations [19, 29]. Unfortunately, even though the equation for the second-order moments can be solved, the equation for the fourth-order moments is very difficult and only approximations or numerical solutions are available (see [15, 28, 47, 50, 57] and [29, Sec. 20.18]).
Here, we consider a secondary scaling regime corresponding to the so-called scintillation regime and in this regime we derive explicit expressions for the fourth-order moments. The scintillation scenario is a well-known paradigm, related to the observation that the irradiance of a star fluctuates due to interaction of the light with the turbulent atmosphere. This common observation is far from being fully understood mathematically. However, experimental observations indicate that the statistical distribution of the irradiance is exponential, with the irradiance being the square magnitude of the complex wave field. Indeed it is a well-accepted conjecture in the physical literature that the statistics of the complex wave field becomes circularly symmetric complex Gaussian when the wave propagates through the turbulent atmosphere [51, 58], so that the irradiance is the sum of the squares of two independent real Gaussian random variables, which has chi-square distribution with two degrees of freedom, that is an exponential distribution. However, so far there is no mathematical proof of this conjecture, except in randomly layered media [18, Chapter 9]. The regime we consider here, which we refer to as the scintillation regime, gives results for the fourth-order moments that are consistent with the scintillation or Gaussian conjecture. We prove in Section 9 that the incoherent zero-mean wave field (i.e. the fluctuations of the wave field defined as the difference between the field and its expectation) has fourth-order moments that obey the Gaussian summation rule. As a result we can discuss the statistical character of the irradiance in detail in Section 10.
Certain functionals of the solution to the white-noise paraxial wave equation can be characterized in some specific regimes [3, 4, 12, 39]. An important aspect of such characterizations is the so-called statistical stability property which corresponds to functionals of the wave field becoming deterministic in the considered scaling regime. This is in particular the case in the limit of rapid decorrelation of the medium fluctuations (in both longitudinal and lateral coordinates). As shown in  the statistical stability also depends on the initial data and can be lost for very rough initial data even with a high lateral diversity as considered there. In [31, 32] the authors also consider a situation with rapidly fluctuating random medium fluctuations and a regime in which the so-called Wigner transform itself is statistically stable. The Wigner transform is known to be a convenient tool to analyze problems involving the Schrödinger equation [27, 43]. In Section 11 we are able to push through a detailed and quantitative analysis of the stability of this quantity using our results on the fourth-order moments. An important aspect of our analysis is that we are able to derive an explicit expression of the coefficient of variation of the smoothed Wigner transform as a function of the smoothing parameters, in the general situation in which the standard deviation can be of the same order as the mean. This is a realistic scenario, we are not deep into a statistical stabilization situation, but in a situation where the parameters of the problem give partly coherent but fluctuating wave functionals. Here we are for the first time able to explicitly quantify such fluctuations and how their magnitude can be controlled by smoothing of the Wigner transform. We believe that these results are important for the many applications where the smoothed Wigner transform appears naturally.
The outline of the paper is as follows: In Section 2 we introduce the Itô-Schrödinger model. In Section 3 we summarize our main results. In Sections 4-5 we describe the general equations for the moments of the field. In Section 6 we discuss the second-order moments. In Section 7 we introduce and analyze the fourth-order moments and the particular parameterization that will be useful to untangle these. In Section 8 we introduce the so-called scintillation regime where we can get an explicit characterization of the fourth-order moments via the main result of the paper presented in Proposition 8.1. Next we discuss three applications of the main result: In Section 9 we prove that the centered fourth-order moments satisfy the Gaussian summation rule, in Section 10 we compute the scintillation index, and in Section 11 we analyze the statistical stability of the smoothed Wigner transform.
2 The White-Noise Paraxial Model
Let us consider the time-harmonic wave equation with homogeneous wavenumber , random index of refraction , and source in the plane :
for and . Denote by the carrier wavelength (equal to ), by the typical propagation distance, and by the radius of the initial transverse source. The paraxial regime holds when the wavelength is much smaller than the radius , and when the propagation distance is smaller than or of the order of (the so-called Rayleigh length). The white-noise paraxial regime that we address in this paper holds when, additionally, the medium has random fluctuations, the typical amplitude of the medium fluctuations is small, and the correlation length of the medium fluctuations is larger than the wavelength and smaller than the propagation distance. In this regime the solution of the time-harmonic wave equation (2.1) can be approximated by 
where is the solution of the Itô-Schrödinger equation
with the initial condition in the plane :
Here the symbol stands for the Stratonovich stochastic integral and is a real-valued Brownian field over with covariance
The model (2.2) can be obtained from the scalar wave equation (2.1) by a separation of scales technique in which the three-dimensional fluctuations of the index of refraction are described by a zero-mean stationary random process with mixing properties: . The covariance function in (2.3) is then given in terms of the two-point statistics of the random process by
The covariance function is assumed to satisfy the following hypothesis:
The condition imposes that the Fourier transform
is continuous and bounded by Lebesgue dominated convergence theorem,
and it is also nonnegative by Bochner’s theorem (it is the power spectral density of a stationary process).
The condition then shows that ,
and therefore is continuous and bounded.
The white-noise paraxial model is widely used in the physical literature. It simplifies the full wave equation (2.1) by replacing it with the initial value-problem (2.2). It was studied mathematically in , in which the solution of (2.2) is shown to be the solution of a martingale problem whose -norm is preserved in the case . The derivation of the Itô-Schrödinger equation (2.2) from the three-dimensional wave equation in randomly scattering medium is given in .
3 Main Result and Quasi Gaussianity
Modeling with the white-noise paraxial model is often motivated by propagation through randomly heterogeneous media. The typical objective for such modeling is to describe some communication or imaging scheme, say with an object buried in the random medium. In many wave propagation and imaging scenarii the quantity of interest is given by a quadratic quantity of the field . For instance, in the time-reversal problems a wave field emitted by the source is recorded on an array, then time-reversed and re-propagated into the medium . The forward and time-reversed propagation paths give rise to a quadratic quantity in the field itself for the re-propagated field. Moreover, in important imaging approaches, in particular passive imaging techniques , the image is formed based on computing cross correlations of the field measured over an array again giving a quadratic expression in the field for the quantity of interest, the correlations. In a number of situations, in particular in optics, the measured quantity is an intensity, again a quadratic quantity in the field. As we explain in Section 6 the expected value of such quadratic quantities can in the white-noise paraxial regime be computed explicitly. In imaging applications this allows to compute the mean image and assess issues like resolution. However, it is important to go beyond this description and calculate the signal-to-noise ratio which requires to compute a fourth-order moment of the wave field. Despite the importance of the signal-to-noise ratio hitherto no rigorous result has been available that accomplishes this task. Indeed explicit expressions for the fourth moments has been a long standing open problem. This is what we push through in this paper. In the context of design of imaging techniques this insight is important to make proper balance in between noise and resolution in the image. We remark that in certain regimes one may be able to prove statistical stability, that is, that the signal-to-noise ratio goes to infinity in the scaling limit [38, 39]. The results we present here are more general in the sense that we can actually describe a finite signal-to-noise ratio and how the parameters of the problem determine this.
To summarize and explicitly articulate the main result regarding the fourth-order moment we consider first the first and second order moments of in (2.2) in the context when . We use the notations for the first and second-order moments
Note that is given explicitly in (6.9). For the second centered moment we use the notation:
Then, to obtain an expression for the fourth-order moment, one heuristic approach often used in the literature [13, 29] is to assume Gaussianity. Consider any complex circularly symmetric Gaussian process then we have  that the fourth-order moment can be expressed in terms of the second-order moments by the Gaussian summation rule as
If the centered field were a complex circularly symmetric Gaussian process, then the fourth-order moment of the field defined by:
This result is not correct in general. For instance, in the spot-dancing regime addressed in [9, 20, 21], the explicit calculation of the moments of all orders is carried out and exhibits non-Gaussian statistics, in particular, the intensity follows a Rice-Nakagami statistics. The spot-dancing regime is valid for a narrow initial beam, strong medium fluctuations, and short propagation distance:
with and the primed quantities of order one.
We show however in this paper that in the so-called scintillation regime the Gaussian summation rule (3.4) is valid . The scintillation regime is discussed in detail in Section 8, it is characterized by a wide initial beam, a long propagation distance, and weak medium fluctuations:
with and the primed quantities of order one. Moreover, in the scintillation regime, if the source spatial profile is Gaussian with radius :
Note that, in the scintillation regime, the field is partially coherent: the coherent field (i.e., the mean field ) has an amplitude which is of the same order as the standard deviation of the zero-mean incoherent field (i.e., the fluctuations of the field ). The surprising result that we report in this paper is that the incoherent field behaves like a random field with Gaussian statistics, as far as the fourth-order moments are concerned.
Finally in the strongly scattering scintillation regime when so that the mean field is vanishing and the field becomes completely incoherent, we have in fact:
These results can now be used to discuss a wide range of applications in imaging and wave propagation. The fourth moment is a fundamental quantity in the context of waves in complex media and the above result is the first rigorous derivation of it that makes explicit the particular scaling regime in which it is valid, moreover, when in fact the Gaussian assumption can be used.
In this paper we also discuss application to characterization of the scintillation in Section 10. The scintillation index describes the relative intensity fluctuations for the wave field. Despite being a fundamental physical quantity associated for instance with light propagation through the atmosphere, a rigorous derivation was not obtained before. We moreover give an explicit characterization of the signal to noise ratio for the Wigner transform in Section 11. The Wigner transform is a fundamental quadratic form of the field that is useful in the context of analysis of problems involving paraxial or Schrödinger equations, for instance time-reversal problems.
We remark finally that the results derived here can be useful in the analysis of ghost imaging experiments [7, 33, 44], enhanced focusing [40, 54, 55, 56] and super-resolution imaging problems [30, 36, 41], and intensity correlation [52, 53]. Results on this will be reported elsewhere.
Ghost imaging is a fascinating recent imaging methodology. It can be interpreted as a correlation-based technique since it gives an image of an object by correlating the intensities measured by two detectors, a high-resolution detector that does not view the object and a low-resolution detector that does view the object. The resolution of the image depends on the coherence properties of the noise sources used to illuminate the object, and on the scattering properties of the medium. This problem can be understood at the mathematical level by using the results presented in this paper.
Enhanced focusing refers to schemes for communication and imaging in a case where a reference signal propagating through the channel is available. Then this information can be used to design an optimal probe that focuses tightly at the desired focusing point. How to optimally design and analyze such schemes, given the limitations of the transducers and so on, can be analyzed using the moment theory presented in this paper. More generally Super resolution refers to the case where one tries to go beyond the classic diffraction limited resolution in imaging systems.
Intensity correlations is a recently proposed scheme for communication in the optical regime that is based on using cross corrections of intensities, as measured in this regime, for communication. This is a promising scheme for communication through relatively strong clutter. By using the correlation of the intensity or speckle for different incoming angles of the source one can get spatial information about the source. The idea of using the information about the statistical structure of speckle to enhance signaling is very interesting and corroborates the idea that modern schemes for communication and imaging require a mathematical theory for analysis of high-order moments.
The results derived in this paper have already opened the mathematical analysis of important imaging problems and we believe that many more problems than those mentioned here will benefit from the results regarding the fourth moments. In fact, enhanced transducer technology and sampling schemes allow for using finer aspects of the wave field involving second- and fourth-order moments and in such complex cases a rigorous mathematical analysis is important to support, complement, or actually disprove, statements based on physical intuition alone.
4 The Mean Field
In this section we give the expression of the mean field, that is, the first-order moment defined by (3.1). Using Itô’s formula for Hilbert space-valued processes  (the process takes values in ), we find that the function satisfies the damped Schrödinger equation in :
with the initial condition . This equation can be solved in the Fourier domain which gives
where is the Fourier transform of the initial field:
In this paper, unless mentioned explicitly, all integrals are over . The exponential damping of the mean field is noticeable, it can be physically explained by the random phase that the wave acquires as it propagates through the random medium. If the initial condition is the Gaussian profile (3.5) then we get
5 The General Moment Equations
The main tool for describing wave statistics are the finite-order moments. We show in this section that in the context of the Itô-Schrödinger equation (2.2) the moments of the field satisfy a closed system at each order [19, 29]. For , we define
for for . Note that here the number of conjugated terms equals the number of non-conjugated terms, otherwise the moments decay relatively rapidly to zero due to unmatched random phase terms associated with random travel time perturbations, as seen in the previous section. Using the stochastic equation (2.2) and Itô’s formula for Hilbert space-valued processes  (the process takes values in ), we find that the function satisfies the Schrödinger-type system in :
with the generalized potential
We introduce the Fourier transform
where is the Fourier transform (4.3) of the initial field and the operator is defined by
where we only write the arguments that are shifted. It turns out that the equation for the Fourier transform is easier to solve than the one for as we will see below.
6 The Second-Order Moments
The second-order moments play an important role, as they give the mean intensity profile and the correlation radius of the transmitted beam [16, 25], they can be used to analyze time reversal experiments [5, 38] and wave imaging problems [10, 11], and we will need them to compute the scintillation index of the transmitted beam and the variance of the Wigner transform. We describe them in detail in this section.
6.1 The Mean Wigner Transform
The second-order moments (3.1) satisfy the system :
starting from . The second-order moment is related to the mean Wigner transform defined by
that is the angularly-resolved mean wave energy density. The mean Wigner transform is in and . It is also bounded by . Using (6.1) we find that it satisfies the closed system
starting from , which is the Wigner transform of the initial field :
Eq. (6.3) has the form of a radiative transport equation for the wave energy density . In this context is the total scattering cross-section and is the differential scattering cross-section that gives the mode conversion rate.
By taking a Fourier transform in and an inverse Fourier transform in of Eq. (6.3):
we obtain a transport equation:
that can be solved and we find the following integral representation for :
where is defined in terms of the initial field as:
6.2 The Mutual Coherence Function
The mutual coherence function is defined by:
where is the mid-point and is the offset. It can be computed by taking the inverse Fourier transform of the expression (6.4):
and we find from (6.7) that the mutual coherence function has the form
7 The Fourth-Order Moments
We consider the fourth-order moment of the field, which is the main quantity of interest in this paper, and parameterize the four points in (5.1) in the special way:
In particular is the barycenter of the four points :
We denote by the fourth-order moment in these new variables:
In the variables the function satisfies the system in :
with the generalized potential
Note in particular that the generalized potential does not depend on the barycenter , and this comes from the fact that the medium is statistically homogeneous. If we assume that the source spatial profile is the Gaussian (3.5) with radius , then the initial condition for Eq. (7.4) is
The Fourier transform (in , , , and ) of the fourth-order moment is defined by: