text:
# Temporal fluctuations of waves in weakly nonlinear disordered media1

## I Introduction

Scattering11footnotetext: Accepted for publication in Physical Review E. You can contact the author at http://www.ilc.msu.su/skipetr/ . of waves in disordered media has proved to be a nontrivial topic possessing intriguing and still not completely understood features. Accordingly to the strength of disorder, one observes a variety of propagation regimes ranging from ballistic transport, through single scattering and wave diffusion, to the Anderson localization. In this paper we are interested in the regime of wave diffusion, corresponding to a relatively strong disorder, which is, however, still largely insufficient to bring the system to the localization transition ($`k\mathrm{}1`$, where $`k`$ is a wave number in the medium, and $`\mathrm{}`$ is a mean free path).

It is well-known, that multiple scattering of coherent wave in a disordered medium results in a complicated spatial intensity distribution $`I(𝐫,t)`$ known as a “speckle pattern”. The speckle pattern is highly irregular and appears random to the eye. It is now well established that the speckle pattern exhibits large intensity fluctuations $`\delta I(𝐫,t)^2I(𝐫,t)^2`$, originating from the interference of partial waves arriving at $`𝐫`$ with completely randomized phases. Here the angular brackets $`\mathrm{}`$ denote ensemble averaging, and $`\delta I(𝐫,t)=I(𝐫,t)I(𝐫,t)`$. Besides, the speckle pattern possesses nontrivial long-range spatial correlation $`C_{\delta I}(𝐫,\mathrm{\Delta }𝐫)=\delta I(𝐫\mathrm{\Delta }𝐫/2,t)\delta I(𝐫+\mathrm{\Delta }𝐫/2,t)`$ even for $`\mathrm{\Delta }r>\mathrm{}`$. This correlation is due to interaction of diffusing modes. If the points $`𝐫\pm \mathrm{\Delta }𝐫/2`$ are far enough from the boundaries of the medium, $`C_{\delta I}1/\mathrm{\Delta }r`$. Recently, it has been shown that in a particular case of a point source of waves embedded inside a disordered medium, there exists an infinite-range contribution to $`C_{\delta I}(𝐫,\mathrm{\Delta }𝐫)`$ originating from scattering events taking place in the immediate neighborhood of the source. This contribution is highly sensitive to the short-distance properties of disorder, as well as to the source size and shape.

If the scatterers in the medium are allowed to move, $`I(𝐫,t)`$ fluctuates with time, and the statistics of these fluctuations is also a subject of active research. A wave propagating in a disordered, multiple-scattering medium undergoes a large number of scattering events, and hence the scattered intensity is highly sensitive to displacements of scatterers. Consequently, the decay of the intensity autocorrelation function $`C_{\delta I}(𝐫,\tau )=\delta I(𝐫,t)\delta I(𝐫,t+\tau )`$ is considerably faster than in the single-scattering case. Recently, long-range autocorrelation function of intensity fluctuations has been measured, and the existence of the universal conductance fluctuations (analogous to that in disordered conductors) has been demonstrated for optical waves. Theoretical analysis of the temporal correlation function of multiple-scattered waves has been extended to amplifying disordered media, as well as to the case of intense incident waves producing flows of scatterers in the disordered medium. An additional contribution to $`C_{\delta I}(𝐫,\tau )`$, originating from scattering in the immediate neighborhood of source and/or detector, decaying much slower than all the previously known contributions, is predicted to exist.

High sensitivity of multiple-scattering speckle patterns to scatterer motion gave rise to a new technique for studying the scatterer dynamics in disordered, turbid media, the so-called “diffusing-wave spectroscopy” (DWS). The latter is now widely applied in concentrated colloidal suspensions, foams, emulsions, granular and biological media. Besides, the DWS has been extended to macroscopically heterogeneous turbid media, providing a tool for imaging of dynamic heterogeneities and visualization of scatterer flows in the bulk of the medium. A generalization of DWS technique has been also accomplished for anisotropic disordered media. Recently, the DWS approach has been extended to nonergodic turbid media.

The above-mentioned, extensive studies of temporal fluctuations of multiple-scattered waves, as well as the numerous application of DWS, are all restricted to linear disordered media. In general, little information is available on the subject of multiple scattering in nonlinear disordered media. Meanwhile, the question concerning the way in which the nonlinearity affects the multiple-scattering speckle pattern still remains open and continue to attract the researches. Considerable efforts have been made to understand the phenomena of coherent backscattering in disordered media with Kerr-type nonlinearity: a narrow dip has been predicted to appear on the top of the backscattering peak. Weak localization effects are shown to exist in the radiation of second harmonic and difference frequency, though their experimental observation failed. Also studied with account for disorder is the optical phase conjugation. More recently, correlations in transmission and reflection coefficients of second harmonic waves have been investigated both theoretically and experimentally, and the angular correlation functions of fundamental wave in a disordered medium with Kerr-type nonlinearity have been calculated. Despite the fact that theoretical description of wave scattering in nonlinear media is complicated by the simultaneous presence of both disorder and nonlinearity, the standard diagram technique for impurity scattering has been extended to the case of disordered medium with nonlinearity of Kerr type.

Very recently, it has been shown that the speckle pattern resulting from the multiple scattering of coherent wave in a nonlinear disordered medium with Kerr-type nonlinearity, should be extremely sensitive to changes of scattering potential, i.e. much more sensitive than the linear speckle pattern. This high sensitivity has been explained by the multiplicity of solutions of nonlinear wave equation. The multiplicity of solutions has been then shown to lead to the temporal instability of the multiple-scattering speckle pattern in nonlinear medium, resulting in spontaneous fluctuations of scattered wave with time. An important prediction of Ref.  is that the nonlinearity should exceed some threshold value for the instability to develop. The threshold value is principally determined by the absorption length $`L_a`$, or by the sample size $`L`$, if $`L<L_a`$. The striking feature of the established result is that the threshold value of nonlinearity tends to zero in an unbounded medium without absorption. Purely elastic, unbounded nonlinear multiple-scattering systems are therefore always unstable. The physical origin of the instability is easy to understand. Nonlinearity modifies the phases of partial waves propagating in the medium. The phase modifications are proportional to the intensity $`I(𝐫,t)`$ and affect the mutual interference of partial waves. As it is this interference which is responsible for $`I(𝐫,t)`$, a sort of feedback establishes in the medium. A small modification of $`I(𝐫,t)`$ causes modifications of phases of partial waves which, in their turn, produce changes of $`I(𝐫,t)`$, and so on. It is well-known that nonlinear wave systems with sufficiently strong, positive feedback may become unstable. As an example, we cite a family of nonlinear optical systems with two-dimensional feedback, where spontaneous formation of complicated spatial structures is observed. Despite the absence of disorder, such systems exhibit transition to seemingly chaotic dynamics with increasing nonlinearity. An analogy can be drawn between the nonlinear optical systems with two-dimensional feedback and nonlinear disordered media by considering the scattering as a (three-dimensional) feedback mechanism. In the case of disordered media, however, the feedback is of random nature and it is therefore hopeless to expect regular spatial structures to form. Meanwhile, the instability can manifest itself in spontaneous fluctuations of speckle pattern. In order to clarify the issue of instability of speckle patterns in nonlinear disordered media, we consider the following questions:

* Can the multiple scattering provide a positive feedback mechanism for waves propagating in a nonlinear disordered medium?
* If “yes”, how strong the nonlinearity should be for the instability to develop?

A general announcement of our principal answers to the above questions has been given in our recent Letter. In the present paper, we discuss and justify the assumptions and approximations underlying our conclusions, provide the missing details of calculations, and give a comprehensive discussion of results. Also developed and discussed is the perturbation approach to the calculation of the temporal autocorrelation function of multiple-scattered wave in a nonlinear disordered medium. It is important that the validity condition of the perturbation theory coincides with the condition for the instability threshold as obtained by using the self-consistent approach. In addition, we give a detailed consideration to an experimentally important case of moving scatterers, when the decrease of the time autocorrelation function is due to a combined effect of spontaneous and scatterer-motion-induced fluctuations of the speckle pattern.

The remainder of the paper is arranged as follows. In Sec. II, we introduce the nonlinear wave equation, and discuss how the path integral approach can be applied for its analysis. We also formulate the basic models and approximations used throughout the paper. Section III is devoted to linear disordered media. In this section, we provide the expressions for the spatio-temporal intensity correlation functions. Although correlations of multiple-scattered waves in linear media are well studied nowadays, we present their first, to our knowledge, treatment with a simultaneous account for absorption, boundary conditions at the sample surface, and scatterer motion. The results of Sec. III serve as a base for further calculations. In Sec. IV, we present a calculation of dephasing of waves in a nonlinear disordered medium. Our calculation takes into account the fluctuations of the local refractive index due to nonlinear effects, as well as the long-range spatial correlation of this fluctuations. Three “nonlinear” contributions to the dephasing are identified in addition to the usual, “linear” term originating directly from the motion of scatterers. Further, in Sec. V we develop a perturbation theory for calculation of the temporal autocorrelation function of a multiple-scattered wave, and show its failure at short correlation times, for sufficiently weak absorption. A condition of validity of the perturbation theory is established by comparing the linear and nonlinear contributions to the dephasing found in Sec. IV. Section VI presents an alternative, self-consistent approach to the calculation of the temporal autocorrelation of scattered wave. Development of self-consistent theory requires some additional assumptions which are also discussed in this section. In Sec. VII, the main results of our self-consistent approach are presented and discussed. The multiple-scattering speckle pattern is shown to exhibit spontaneous fluctuations even in the absence of scatterer motion, which we interpret as a signature of its instability. A comparison of self-consistent and perturbative results is given, and the condition of the speckle pattern instability is shown to coincide with the condition of validity of the perturbation theory. Finally, concluding remarks are presented in Sec. VIII. In order to maintain the text of the paper readable, we have chosen to collect the technical details of calculations in three appendices. Appendix A is devoted to the derivation of the field-field spatio-temporal correlation function. In Appendix B we compute the spatio-temporal long-range intensity correlation function. Appendix C provides the details of calculations of path distributions $`\rho _s(𝐫)`$ and $`\rho _s(𝐫,𝐫^{})`$ defined in Sec. IV.

## II Wave equation and path integrals

We consider a scalar monochromatic wave of frequency $`\omega `$ propagating in a random medium with Kerr-type nonlinearity. The wave amplitude $`\psi (𝐫,t)`$ obeys a nonlinear wave equation:

$`\left\{^2+k_0^2\left[\epsilon _0^{}+i\epsilon _0^{\prime \prime }+\delta \epsilon (𝐫,t)+\epsilon _2\left|\psi (𝐫,t)\right|^2\right]\right\}\psi (𝐫,t)`$ (1)
$`=0.`$ (2)
Here $`k_0`$ is the free-space wave number, $`\epsilon _0=\epsilon _0^{}+i\epsilon _0^{\prime \prime }`$ is the average (complex) dielectric function, $`\delta \epsilon (𝐫,t)`$ is the fluctuating part of the dielectric function, $`\epsilon _2`$ is the nonlinear susceptibility (the two latter quantities are assumed to be real). Eq. (2) is valid only if $`\delta \epsilon (𝐫,t)+\epsilon _2\left|\psi (𝐫,t)\right|^2`$ does not change significantly on the time scale of $`\omega ^1`$. The expression in the square brackets of Eq. (2) can be considered as some “effective” dielectric function of the medium. General analysis of Eq. (2) for arbitrary relation between various terms comprising this function constitutes a formidable task, and is not a purpose of this paper. We assume the following hierarchy:

$`\epsilon _2\left|\psi (𝐫,t)\right|^2\delta \epsilon (𝐫,t),\text{ }\epsilon _2\left|\psi (𝐫,t)\right|^2\epsilon _0^{},\text{ }\epsilon _0^{\prime \prime }\epsilon _0^{}.`$ (3)
In other words, we assume that the role of nonlinearity is less significant than that of disorder, and that absorption is weak allowing multiple scattering of waves in the medium. It is then convenient to define the effective refractive index $`n_0=(\epsilon _0^{})^{1/2}`$, the absorption length $`\mathrm{}_a=n_0/(k_0\epsilon _0^{\prime \prime })`$, and the nonlinear coefficient $`n_2=\epsilon _2/(2n_0)`$, which determines the nonlinear correction to the (linear) refractive index of the medium: $`n(𝐫,t)=n_0+n_2I(𝐫,t)`$, where $`I(𝐫,t)=\left|\psi (𝐫,t)\right|^2`$ is the wave intensity.

In this paper, we study the fluctuations of the solution of Eq. (2) with time $`t`$. In a linear medium ($`\epsilon _2=0`$), these fluctuations can only be due to random fluctuations of $`\delta \epsilon (𝐫,t)`$ with time. The fluctuations of $`\psi (𝐫,t)`$ are commonly characterized by the autocorrelation function $`C_\psi (𝐫,\tau )=\psi (𝐫,t)\psi ^{}(𝐫,t+\tau )`$. We assume that this autocorrelation function is independent of $`t`$ which implies that for fixed $`𝐫`$, $`\psi (𝐫,t)`$ represents a stationary random process (this is obviously true if the sample geometry and the source distribution do not change with time, and if $`\delta \epsilon (𝐫,t)`$ is a stationary random process). Me take $`\delta \epsilon (𝐫,t)`$ to be a Gaussian random field with zero mean and the correlation function $`C_{\delta \epsilon }(\mathrm{\Delta }𝐫,\tau )=\delta \epsilon (𝐫\mathrm{\Delta }𝐫/2,t)\delta \epsilon (𝐫+\mathrm{\Delta }𝐫/2,t+\tau )`$. For a medium composed of point-like scatterers undergoing Brownian motion with a diffusion coefficient $`D_B`$,

$`C_{\delta \epsilon }(\mathrm{\Delta }𝐫,\tau )={\displaystyle \frac{4\pi /(k^4\mathrm{})}{(4\pi D_B\tau )^{3/2}}}\mathrm{exp}\left({\displaystyle \frac{\mathrm{\Delta }r^2}{4D_B\tau }}\right),`$ (4)
where $`k=k_0n_0`$, and the mean free path $`\mathrm{}\mathrm{}_a`$ is introduced (a weak scattering limit $`k\mathrm{}1`$ is assumed). A natural time scale for scattering of waves in the medium described by Eq. (4) is set by the characteristic time needed for a scatterer to move a distance of order of the wavelength: $`\tau _0=(4k^2D_B)^1`$. From here on, we will be interested in short correlation times $`\tau \tau _0`$.

In the linear case ($`\epsilon _2=0`$), several approaches have been elaborated to analyze Eq. (2). We mention the diagrammatic techniques, theory of radiative transfer, and the method of path integrals. The three approaches are known to give equivalent results for $`C_\psi `$ at $`\tau \tau _0`$. In the present paper, we adopt the method of path integrals that was originally proposed in the framework of quantum electrodynamics, but later has been successfully used in various areas of physics, and, in particular, for the analysis of wave scattering problems. The method is based on the fact that the solution $`\psi (𝐫,t)`$ of the wave equation (2) can be written in a form of a functional integral, with integration performed over all possible trajectories (paths) going from the source of waves to $`𝐫`$. Since in the weak scattering limit ($`k\mathrm{}1`$) different trajectories can be considered independent, it appears that the correlation function $`C_\psi `$ reduces to the following integral:

$`C_\psi (𝐫,\tau )=I_0{\displaystyle _0^{\mathrm{}}}P(𝐫,s)\mathrm{exp}\left[{\displaystyle \frac{1}{2}}\mathrm{\Delta }\phi ^2(\tau )_s\right]𝑑s,`$ (5)
where $`I_0`$ is the average intensity in a nonabsorbing medium, $`P(𝐫,s)`$ is a weight coefficient of paths of length $`s`$, and $`\mathrm{\Delta }\phi ^2(\tau )_s`$ denotes the squared phase difference $`\mathrm{\Delta }\phi (t,\tau )=\phi (t+\tau )\phi (t)`$, averaged over various realizations of disorder, and over all possible paths of the same length $`s`$. From here on, we denote such an averaging by $`\mathrm{}_s`$. Note that $`\mathrm{\Delta }\phi (\tau )_s=0`$ for the model of Brownian point-like scatterers. Meanwhile,

$`\mathrm{\Delta }\phi ^2(\tau )_s^{(0)}={\displaystyle \frac{\tau }{\tau _0}}\times {\displaystyle \frac{s}{\mathrm{}}},`$ (6)
where the superscript $`(0)`$ denotes the linear case. It is worth noting that $`\mathrm{\Delta }\phi ^2(\tau )_s^{(0)}`$ does not depend neither on the sample geometry, nor on the source and detector positions. Its value is only determined by the scatterer dynamics (through the single-scattering correlation time $`\tau _0`$), and the path length $`s`$. In contrast, $`P(𝐫,s)`$ can only be calculated if the sample geometry, source distribution, and detector position $`𝐫`$ are specified. In what follows, we restrict our analysis to a semi-infinite medium occupying the half-space $`z>0`$, and illuminated by a plane monochromatic wave incident at $`z=0`$. For $`s\mathrm{}`$, $`P(𝐫,s)`$ becomes

$`P(𝐫,s)=\left({\displaystyle \frac{3z^2}{4\pi \mathrm{}s^3}}\right)^{1/2}\mathrm{exp}\left({\displaystyle \frac{3z^2}{4s\mathrm{}}}{\displaystyle \frac{s}{\mathrm{}_a}}\right).`$ (7)
Once the field correlation $`C_\psi (𝐫,\tau )`$ is known, the autocorrelation of intensity fluctuations $`C_{\delta I}(𝐫,\tau )=\delta I(𝐫,t)\delta I(𝐫,t+\tau )`$ can be found applying the factorization approximation: $`C_{\delta I}(𝐫,\tau )=|C_\psi (𝐫,\tau )|^2`$.

Combining together Eqs. (5)–(7), one obtains the normalized autocorrelation function of multiple-scattered wave in a semi-infinite disordered medium:

$`g_1^{(L)}(𝐫,\tau )={\displaystyle \frac{C_\psi (𝐫,\tau )}{C_\psi (𝐫,0)}}=\mathrm{exp}\left\{\left[\alpha \left(\tau \right){\displaystyle \frac{\mathrm{}}{L_a}}\right]{\displaystyle \frac{z}{\mathrm{}}}\right\},`$ (8)
where the superscript $`(L)`$ denotes the linear case, $`\alpha ^2(\tau )=3\tau /(2\tau _0)+\mathrm{}^2/L_a^2`$, and $`L_a=(\mathrm{}\mathrm{}_a/3)^{1/2}\mathrm{}`$. For the diffusely reflected wave, we assume $`z\mathrm{}`$ and get

$`g_1^{(L)}(\mathrm{},\tau )g_1^{(L)}(\tau )=\mathrm{exp}\left\{\alpha \left(\tau \right)+{\displaystyle \frac{\mathrm{}}{L_a}}\right\}.`$ (9)
From here on, we will use $`g_1(𝐫,\tau )`$ to denote the normalized autocorrelation function at a point $`𝐫`$ inside the medium, while $`g_1(\tau )`$ — for the normalized autocorrelation function of diffusely reflected wave. The superscripts $`(L)`$ and $`(NL)`$ will be used to distinguish between linear and nonlinear cases.

Now we turn to the nonlinear medium. Strictly speaking, the method of path integrals cannot be applied for the analysis of Eq. (2), once $`\epsilon _20`$. The failure of the path-integral technique follows from the fact that this approach relies on the superposition principle which is not valid for waves in nonlinear media. However, if the nonlinearity is weak \[which is ensured by the first two inequalities of Eq. (3)\], Eq. (5) is still approximately valid provided that its main ingredients $`P(𝐫,s)`$ and $`\mathrm{\Delta }\phi ^2(\tau )_s`$ are computed with account for nonlinear effects. To simplify such a calculation, we assume that the nonlinearity is sufficiently weak to validate the following two assumptions:

* Propagation of waves in a weakly nonlinear disordered medium is diffusive with a mean free path $`\mathrm{}`$ unaffected by the nonlinearity. This implies that nonlinear refraction is negligible at distances of order $`\mathrm{}`$, and consequently, that $`\mathrm{\Delta }n^2k\mathrm{}1`$, where $`\mathrm{\Delta }n=n_2I_0`$, and $`I_0`$ is the average intensity in the absence of absorption. This assumption is an alternative formulation of the fact that the role of nonlinearity is much less significant than that of disorder \[see also the first two inequalities of Eq. (3)\].
* Intensity of the third harmonic remains always much smaller than the intensity of the fundamental wave. This implies either that $`\psi (𝐫,t)`$ is considered as a complex quantity \[in this case, Eq. (2) is a nonlinear Schrödinger equation, $`\left|\psi (𝐫,t)\right|^2`$ is time-independent, and the third harmonic is not generated at all\], or that the medium has a sufficient degree of dispersion for the phase matching condition to be violated: $`\left|\mathrm{\Delta }k\right|\mathrm{}1`$ with $`\mathrm{\Delta }k=k_33k`$ and $`k_3`$ being the wave number at frequency $`3\omega `$.

Assumption (i) allows us to consider the path distribution $`P(𝐫,s)`$ being unaffected by nonlinearity. The only object to be recalculated with account for nonlinear effects is then $`\mathrm{\Delta }\phi ^2(\tau )_s`$. Before going into an explicit calculation of $`\mathrm{\Delta }\phi ^2(\tau )_s`$, we devote the next section to a brief derivation of some important results for linear medium.

## III Correlations in a linear medium

As indicated above, we consider a monochromatic plane wave incident at the surface $`z=0`$ of a semi-infinite medium occupying the half-space $`z>0`$. The average intensity at $`z>\mathrm{}`$ can be then found in the diffusion approximation: $`I(𝐫,t)=I_0\mathrm{exp}(z/L_a)`$. The spatio-temporal correlation function of the field is given by a solution of the Bethe-Salpeter equation (see Appendix A for details of the calculation):

$`C_\psi (𝐫,\mathrm{\Delta }𝐫,\tau )`$ $`=`$ $`\psi (𝐫\mathrm{\Delta }𝐫/2,t)\psi ^{}(𝐫+\mathrm{\Delta }𝐫/2,t+\tau )`$ (10)
$`=`$ $`I_0{\displaystyle \frac{\mathrm{sin}(k\mathrm{\Delta }r)}{k\mathrm{\Delta }r}}\mathrm{exp}\left\{{\displaystyle \frac{\mathrm{\Delta }r}{2\mathrm{}}}\alpha \left(\tau \right){\displaystyle \frac{z}{\mathrm{}}}\right\}.`$ (11)
Let us now consider the correlation functions of intensity fluctuations. In addition to $`\delta I(𝐫,t)`$, which is the deviation of intensity from its average value, it is convenient to define $`\mathrm{\Delta }I(𝐫,t,\tau )=I(𝐫,t+\tau )I(𝐫,t)`$, which is the change of the local intensity during the time interval $`\tau `$. While $`\delta I(𝐫,t)=\mathrm{\Delta }I(𝐫,t,\tau )=0`$, the correlation functions $`C_{\delta I}(𝐫,\mathrm{\Delta }𝐫,\tau )`$ and $`C_{\mathrm{\Delta }I}(𝐫,\mathrm{\Delta }𝐫,\tau )=\mathrm{\Delta }I(𝐫\mathrm{\Delta }𝐫/2,t,\tau )\mathrm{\Delta }I(𝐫+\mathrm{\Delta }𝐫/2,t,\tau )`$ for $`\mathrm{\Delta }r<\mathrm{}`$ can be found in the factorization approximation:

$`C_{\delta I}(𝐫,\mathrm{\Delta }𝐫,\tau )=\left|C_\psi (𝐫,\mathrm{\Delta }𝐫,\tau )\right|^2,`$ (12)
$`C_{\mathrm{\Delta }I}(𝐫,\mathrm{\Delta }𝐫,\tau )=2\left[C_{\delta I}(𝐫,\mathrm{\Delta }𝐫,0)C_{\delta I}(𝐫,\mathrm{\Delta }𝐫,\tau )\right].`$ (13)
Both correlation functions (12) and (13) decrease exponentially with $`\mathrm{\Delta }r/\mathrm{}`$, and thus become negligible for $`\mathrm{\Delta }r>\mathrm{}`$. Intensity correlation persists, however, even for two points separated by a distance $`\mathrm{\Delta }r>\mathrm{}`$. This correlation is due to the diffusive nature of wave transport in the medium and can be found either using the Langevin approach or applying diagrammatic methods. We give the details of calculations in Appendix B, the final results are

$`C_{\delta I}(𝐫,\mathrm{\Delta }𝐫,\tau )={\displaystyle \frac{3}{(k\mathrm{})^2}}I_0^2{\displaystyle _0^{\mathrm{}}}𝑑KK`$ (14)
$`\times Q[K,\sqrt{K^2+\mathrm{}^2/L_a^2},{\displaystyle \frac{z}{\mathrm{}}},{\displaystyle \frac{\mathrm{\Delta }z}{\mathrm{}}},\alpha \left(\tau \right)]`$ (15)
$`\times J_0\left(K\mathrm{\Delta }R/\mathrm{}\right),`$ (16)
$`C_{\mathrm{\Delta }I}(𝐫,\mathrm{\Delta }𝐫,\tau )={\displaystyle \frac{6}{(k\mathrm{})^2}}I_0^2{\displaystyle _0^{\mathrm{}}}𝑑KK`$ (17)
$`\times \mathrm{\Delta }Q[K,\sqrt{K^2+\mathrm{}^2/L_a^2},{\displaystyle \frac{z}{\mathrm{}}},{\displaystyle \frac{\mathrm{\Delta }z}{\mathrm{}}},\alpha \left(\tau \right),\alpha (0)]`$ (18)
$`\times J_0\left(K\mathrm{\Delta }R/\mathrm{}\right).`$ (19)
Here we use the cylindrical coordinates: $`𝐫=\{𝐑,z\}`$, $`J_0`$ is the Bessel function of zeroth order, the function $`Q`$ is defined in Appendix B, and

$`\mathrm{\Delta }Q[\mathrm{},\alpha \left(\tau \right),\alpha (0),]`$ $`=`$ $`Q[\mathrm{},\alpha (0)]`$ (20)
$``$ $`Q[\mathrm{},\alpha \left(\tau \right)].`$ (21)
Due to a rather complicated structure of the function $`Q`$ (see Appendix B), further calculations can be done only approximately. In the case of $`z\pm \mathrm{\Delta }z/2,\mathrm{\Delta }r\mathrm{}/\alpha (\tau )`$, we get

$`C_{\delta I}(𝐫,\mathrm{\Delta }𝐫,\tau )`$ $``$ $`{\displaystyle \frac{3}{4}}{\displaystyle \frac{I_0^2}{(k\mathrm{})^2}}{\displaystyle \frac{\mathrm{}}{z_g}}\left[{\displaystyle \frac{1}{\chi }}{\displaystyle \frac{1}{\sqrt{1+\chi ^2}}}\right]`$ (22)
$`\times `$ $`\mathrm{exp}\left[2\alpha \left(\tau \right){\displaystyle \frac{z}{\mathrm{}}}\right],`$ (23)
where $`z_g`$ is the geometrical average of $`z`$-coordinates of the two points for which the correlation is computed: $`z_g=\sqrt{z^2\mathrm{\Delta }z^2/4}`$, and $`\chi =\mathrm{\Delta }r/(2z_g)`$. For $`\alpha (\tau )=0`$, Eq. (23) is exact. If $`\chi <1`$, the correlation behaves essentially as $`1/\mathrm{\Delta }r`$, while for $`\chi >1`$ it becomes proportional to $`z_g^2/\mathrm{\Delta }r^3`$. For the correlation function of Eq. (19) we find:

$`C_{\mathrm{\Delta }I}(𝐫,\mathrm{\Delta }𝐫,\tau )`$ (24)
$`\{\begin{array}{cc}4\left[\alpha \left(\tau \right)\alpha (0)\right](z/\mathrm{})C_{\delta I}(𝐫,\mathrm{\Delta }𝐫,0),\hfill & \\ 2\alpha (\tau )(z/\mathrm{})1\hfill & \\ 2C_{\delta I}(𝐫,\mathrm{\Delta }𝐫,0),2\alpha (\tau )(z/\mathrm{})1\hfill & \end{array}`$ (25)
## IV Dephasing of waves in a nonlinear medium

Consider a single wave path of length $`s`$ going from the source of waves to some point $`𝐫`$. The phase acquired by a wave traveling along such a path can be written as

$`\phi (t)={\displaystyle _0^s}k_0n[𝐫(s_1),t]𝑑s_1,`$ (26)
where the integration is along the path, and $`n(𝐫,t)=n_0+n_2I(𝐫,t)`$. The squared difference $`\mathrm{\Delta }\phi (t,\tau )=\phi (t+\tau )\phi (t)`$, averaged over various realizations of disorder, and over all possible paths of length $`s`$, is found directly from Eq. (26):

$`\mathrm{\Delta }\phi ^2(\tau )_s={\displaystyle \underset{j=0}{\overset{3}{}}}\mathrm{\Delta }\phi ^2(\tau )_s^{(j)},`$ (27)
where the four contributions corresponding to $`j=0,\mathrm{},3`$ originate from different physical processes. Below, we give explicit expressions of these terms and discuss their origin.

The first term in Eq. (27), $`\mathrm{\Delta }\phi ^2(\tau )_s^{(0)}`$ is the linear term given by Eq. (6). The next three terms, namely the terms corresponding to $`j=1,2,3`$, are absent in the linear case and only appear because of nonlinear nature of wave interaction with the medium. Explicit expressions for these terms are

$`\mathrm{\Delta }\phi ^2(\tau )_s^{(1)}`$ $`=`$ $`{\displaystyle \frac{2n_2}{n_0\mathrm{}}}{\displaystyle \frac{\tau }{\tau _0}}{\displaystyle _0^s}I(𝐫)𝑑s_1_s`$ (28)
$`=`$ $`{\displaystyle \frac{2n_2}{n_0}}{\displaystyle \frac{\tau }{\tau _0}}I(𝐫)_s{\displaystyle \frac{s}{\mathrm{}}},`$ (29)
$`\mathrm{\Delta }\phi ^2(\tau )_s^{(2)}`$ $`=`$ $`{\displaystyle \frac{\pi }{n_0}}k_0n_2^2{\displaystyle _0^s}C_{\mathrm{\Delta }I}(𝐫,0,\tau )𝑑s_1_s`$ (30)
$`=`$ $`{\displaystyle \frac{\pi }{n_0}}k_0\mathrm{}n_2^2C_{\mathrm{\Delta }I}(𝐫,0,\tau )_s{\displaystyle \frac{s}{\mathrm{}}},`$ (31)
$`\mathrm{\Delta }\phi ^2(\tau )_s^{(3)}`$ $`=`$ $`k_0^2n_2^2{\displaystyle _0^s}{\displaystyle _0^s}C_{\mathrm{\Delta }I}(𝐫,\mathrm{\Delta }𝐫,\tau )𝑑s_1𝑑s_2_s=`$ (32)
$`=`$ $`(k_0\mathrm{})^2n_2^2C_{\mathrm{\Delta }I}(𝐫,\mathrm{\Delta }𝐫,\tau )_s\left({\displaystyle \frac{s}{\mathrm{}}}\right)^2.`$ (33)
Here the integrations are assumed along wave paths of length $`s\mathrm{}`$ \[in Eq. (33), both integrals are along the same path\]. Eq. (31) originates from the short-range correlation of intensity fluctuations \[see Eqs. (11)–(13)\]:

$`\mathrm{\Delta }\phi ^2(\tau )_s^{(2)}`$ $`=`$ $`k_0^2n_2^2C_{\mathrm{\Delta }I}(𝐫,0,\tau )_s{\displaystyle _0^s}𝑑s{\displaystyle _{\mathrm{}}^{\mathrm{}}}d(\mathrm{\Delta }s)`$ (35)
$`\left[{\displaystyle \frac{\mathrm{sin}(k\mathrm{\Delta }s)}{k\mathrm{\Delta }s}}\right]^2\mathrm{exp}\left({\displaystyle \frac{\mathrm{\Delta }s}{\mathrm{}}}\right),`$
where the wave path is assumed to be ballistic at distances shorter than $`\mathrm{}`$. Eq. (35) reduces to (31) for $`k\mathrm{}1`$. Next, the term given by Eq. (33) is due to the long-range correlations of intensity fluctuations \[see Eq. (19)\]. The averages entering into the right-hand sides of Eqs. (29)–(35) are

$`I(𝐫)_s={\displaystyle d^3𝐫\rho _s(𝐫)I(𝐫)},`$ (36)
$`C_{\mathrm{\Delta }I}(𝐫,0,\tau )_s={\displaystyle d^3𝐫\rho _s(𝐫)C_{\mathrm{\Delta }I}(𝐫,0,\tau )},`$ (37)
$`C_{\mathrm{\Delta }I}(𝐫,\mathrm{\Delta }𝐫,\tau )_s={\displaystyle d^3𝐫_1d^3𝐫_2\rho _s(𝐫_1,𝐫_2)}`$ (38)
$`C_{\mathrm{\Delta }I}(𝐫,\mathrm{\Delta }𝐫,\tau ),`$ (39)
where the integrations are over the volume of the disordered medium, and $`𝐫_{1,2}=𝐫\pm \mathrm{\Delta }𝐫/2`$. In Eqs. (36)–(39), $`\rho _s(𝐫)`$ is the probability density for a path of length $`s`$, going from the source of waves to the point $`𝐫`$, to pass through the vicinity of $`𝐫`$, and $`\rho _s(𝐫_1,𝐫_2)`$ is the probability density for the path to pass consequently through the vicinities of $`𝐫_1`$ and $`𝐫_2`$. These two “path distributions” should be calculated with account for a particular geometry of disordered sample and source of waves. Once the geometry is fixed, the calculation is straightforward. For the case of a plane wave incident upon a semi-infinite disordered medium, the calculations of $`\rho _s(𝐫)`$ and $`\rho _s(𝐫_1,𝐫_2)`$ are presented in Appendix C.

Let us discuss briefly the physical origins of nonlinear contributions to the dephasing given by Eqs. (29)–(33). $`\mathrm{\Delta }\phi ^2(\tau )_s^{(1)}`$ describes the change of the effective wave number in the nonlinear medium: $`k=k_0n_0k(𝐫)=k_0[n_0+n_2I(𝐫,t)]`$. This contribution can be either positive or negative, depending on the sign of $`n_2`$, but its absolute value is always much smaller than $`\mathrm{\Delta }\phi ^2(\tau )_s^{(0)}`$, as long as $`\mathrm{\Delta }nn_0`$. $`\mathrm{\Delta }\phi ^2(\tau )_s^{(1)}`$ can therefore only cause a small correction to the linear correlation function. The next contribution, $`\mathrm{\Delta }\phi ^2(\tau )_s^{(2)}`$, originates from fluctuations of the local intensity, while $`\mathrm{\Delta }\phi ^2(\tau )_s^{(3)}`$ is due to the long-range correlation of these fluctuations. An important difference between the linear term (6), the first nonlinear term (29), and the two last nonlinear terms (31), (33), is that the latter terms do not depend explicitly on $`\tau _0`$. The terms given by Eqs. (31) and (33) are determined by the intensity fluctuations, and not by the scatterer displacements. This might seem to be a meaningless statement, as the intensity fluctuations are, in their turn, caused by the scatterer motion. The important point is that the scatterer motion is not the only possible reason for the fluctuations of intensity with time. Weak, spontaneous fluctuations of $`I(𝐫,t)`$ (due to thermal fluctuations of various parameters, vibrations in the experimental setup, fluctuations of the incident wave, etc.) are inevitable in real physical systems. Eqs. (31) and (33) provide a mechanism for this spontaneous and generally weak fluctuations to affect the dephasing $`\mathrm{\Delta }\phi ^2(\tau )_s`$ and, consequently, the temporal correlation function of scattered wave.

## V Perturbation theory

As stated in the title, the present paper is devoted to weakly nonlinear disordered media. We limit ourselves to a weak nonlinearity, as otherwise the problem becomes too involved. Above, we have already mentioned that we assume $`\mathrm{\Delta }n^2k\mathrm{}1`$, and that this condition allows us to consider the transport of average intensity to remain unaffected by the nonlinearity \[assumption (i) of Sec. II\]. This allows us to use “linear” results, $`I(𝐫)=I_0\mathrm{exp}(z/L_a)`$ and Eqs. (C7), (C16) of Appendix C, for $`I(𝐫)`$, $`\rho _s(𝐫)`$, and $`\rho _s(𝐫_1,𝐫_2)`$ in Eqs. (36)–(39). It seems then natural to assume that $`C_{\mathrm{\Delta }I}(𝐫,0,\tau )`$ and $`C_{\mathrm{\Delta }I}(𝐫,\mathrm{\Delta }𝐫,\tau )`$ are also close to their linear values. We can therefore replace these correlation functions in Eqs. (37) and (39) by the expressions found in Sec. III. Then, making use of Eqs. (12)–(19), and performing necessary integrations, we obtain from Eqs. (29)–(33):

$`\mathrm{\Delta }\phi ^2(\tau )_s^{(1)}`$ $`=`$ $`2\mathrm{\Delta }n{\displaystyle \frac{\tau }{\tau _0}}\left\{1H\left[\alpha (0)\sqrt{{\displaystyle \frac{s}{12\mathrm{}}}}\right]\right\}{\displaystyle \frac{s}{\mathrm{}}},`$ (40)
$`\mathrm{\Delta }\phi ^2(\tau )_s^{(2)}`$ $`=`$ $`2\pi k_0\mathrm{}\mathrm{\Delta }n^2\{H\left[\alpha (\tau )\sqrt{{\displaystyle \frac{s}{3\mathrm{}}}}\right]`$ (41)
$``$ $`H\left[\alpha (0)\sqrt{{\displaystyle \frac{s}{3\mathrm{}}}}\right]\}{\displaystyle \frac{s}{\mathrm{}}},`$ (42)
$`\mathrm{\Delta }\phi ^2(\tau )_s^{(3)}`$ $`=`$ $`6\mathrm{\Delta }n^2S[\alpha (\tau )\sqrt{s/\mathrm{}},\alpha (0)\sqrt{s/\mathrm{}}]`$ (43)
$`\times `$ $`\left({\displaystyle \frac{s}{\mathrm{}}}\right)^{3/2}.`$ (44)
Here $`H(x)=\sqrt{\pi }x\mathrm{exp}(x^2)[1\mathrm{Erf}(x)]`$ and $`n_0=1`$ is assumed for simplicity. The function $`S(u,v)`$ in Eq. (44) is

$`S(u,v)`$ $`=`$ $`9{\displaystyle _0^{\mathrm{}}}𝑑RR{\displaystyle _0^{\mathrm{}}}𝑑KK{\displaystyle _0^{\mathrm{}}}𝑑z{\displaystyle _0^{2z}}d(\mathrm{\Delta }z)`$ (45)
$`\times `$ $`f(z,\mathrm{\Delta }z,R)\mathrm{\Delta }Q(K,\sqrt{K^2+v^2},z,\mathrm{\Delta }z,u,v)`$ (46)
$`\times `$ $`J_0(KR),`$ (47)
$`f(z,\mathrm{\Delta }z,R)`$ $`=`$ $`{\displaystyle \frac{2z+\sqrt{R^2+\mathrm{\Delta }z^2}}{\sqrt{R^2+\mathrm{\Delta }z^2}}}`$ (48)
$`\times `$ $`\mathrm{exp}\left[{\displaystyle \frac{3}{4}}\left(2z+\sqrt{R^2+\mathrm{\Delta }z^2}\right)^2\right]`$ (49)
$``$ $`{\displaystyle \frac{2z+\sqrt{R^2+4z^2}}{\sqrt{R^2+4z^2}}}`$ (50)
$`\times `$ $`\mathrm{exp}\left[{\displaystyle \frac{3}{4}}\left(2z+\sqrt{R^2+4z^2}\right)^2\right].`$ (51)
Unfortunately, integrations in Eq. (47) cannot be performed in the general case. We find, however, the following approximate results:

$`S(u,v)\{\begin{array}{cc}(uv),\hfill & uv1,v1\text{,}\hfill \\ 1,\hfill & uv>1,v1\text{,}\hfill \\ (uv)v^3,\hfill & uv1,v>1\text{,}\hfill \\ v^3,\hfill & uv>1,v>1\text{.}\hfill \end{array}`$ (52)
Here numerical factors of order unity are omitted before each of the four asymptotic expressions.

While approximate, the above results enable one to compute the temporal autocorrelation function $`g_1^{(NL)}(\tau )`$ of diffusely reflected wave numerically using Eqs. (5), (7), and (27):

$`g_1^{(NL)}(\tau )=W\left[\mathrm{\Delta }\phi ^2(\tau )_s\right]/W[0],`$ (53)
$`W\left[\mathrm{\Delta }\phi ^2(\tau )_s\right]`$ (54)
$`={\displaystyle _0^{\mathrm{}}}P(\mathrm{},s)\mathrm{exp}\left[{\displaystyle \frac{1}{2}}\mathrm{\Delta }\phi ^2(\tau )_s\right]𝑑s.`$ (55)
We present the results of the calculation in Fig. 1 for fixed $`\mathrm{\Delta }n`$, $`k_0\mathrm{}`$, and the values of the inverse absorption length $`\mathrm{}/L_a`$ indicated near each curve. The “linear” correlation functions, corresponding to the same three values of the absorption length, and to $`\mathrm{\Delta }n=0`$, are shown by dashed lines. For weak absorption ($`\mathrm{}/L_a=0`$, $`10^3`$), the initial (i.e. at short correlation delay times $`\tau `$) decrease of the “nonlinear” autocorrelation function is much faster than that of the linear one. Hence, our perturbation approach fails at short times. Indeed, the results of this section are based on the assumption that $`C_{\mathrm{\Delta }I}(𝐫,0,\tau )`$ and $`C_{\mathrm{\Delta }I}(𝐫,\mathrm{\Delta }𝐫,\tau )`$ in the nonlinear medium are close to their values in the linear one. Applying the factorization approximation, we find that this implies that $`1|g_1^{(NL)}(\tau )|^21|g_1^{(L)}(\tau )|^2`$. Although this condition seems to hold well for sufficiently large $`\tau `$, it can be violated at small $`\tau `$, where, as follows from Fig. 1, $`1|g_1^{(NL)}(\tau )|^2`$ can become much larger than $`1|g_1^{(L)}(\tau )|^2`$, if absorption is sufficiently weak.

To estimate the region of validity of our perturbation approach, we require that the linear contribution to the dephasing $`\mathrm{\Delta }\phi ^2(\tau )_s`$ given by Eq. (6), should be considerably greater than the sum of nonlinear contributions \[Eqs. (40)–(44)\]. As the longest path length contributing to the integral of Eqs. (5) and (55) is $`s\mathrm{}/\alpha ^2`$, we obtain the conditions of validity of the perturbation theory in the form:

$`\mathrm{\Delta }n^2\alpha (\tau )^2\left[k_0\mathrm{}+\alpha (\tau )^1\right]1.`$ (56)
Since $`\alpha (\tau )`$ is an increasing function of its argument, and $`\alpha (0)=\mathrm{}/L_a`$, condition $`\mathrm{\Delta }n^2(L_a/\mathrm{})^2[k_0\mathrm{}+L_a/\mathrm{}]1`$ ensures Eq. (56) at any $`\tau `$. It is the case for the upper curve of Fig. 1, corresponding to $`\mathrm{}/L_a=5\times 10^3`$. For such a strong absorption, the perturbation theory is valid at any $`\tau `$, and the nonlinear autocorrelation function is close to the linear result. By contrast, if $`\mathrm{\Delta }n^2(L_a/\mathrm{})^2[k_0\mathrm{}+L_a/\mathrm{}]>1`$, the perturbation approach can be applied only for sufficiently long correlation times $`\tau >\tau _c`$, where the critical time $`\tau _c`$ is determined by Eq. (56). For the lower curve of Fig. 1, corresponding to $`\mathrm{}/L_a=0`$, we find $`(\tau _c/\tau _0)^{1/2}2\times 10^3`$. It is worthwhile to note that in the absence of absorption, even an infinitely small nonlinearity ($`\mathrm{\Delta }n0`$) suffices for the perturbation theory to fail at $`\tau /\tau _00`$. Our condition of validity of the perturbation approach (56) is consistent with the result of Spivak and Zyuzin, who have shown that the perturbation analysis of the sensitivity of speckle pattern in a nonlinear disordered medium to changes of scattering potential fails for $`\mathrm{\Delta }n^2(L/\mathrm{})^3>1`$, where $`L\mathrm{}`$ is the typical size of the medium.

## VI Self-consistent analysis

As demonstrated in the previous section, the perturbation theory fails to describe the temporal autocorrelation function of wave diffusely reflected from a nonlinear medium for $`\tau <\tau _c`$, where $`\tau _c`$ is defined by Eq. (56). To calculate $`g_1^{(NL)}(\tau )`$ at such short times, one has to use a non-perturbative, self-consistent analysis. However, the possibility of performing such an analysis is considerably limited by the mathematical complexity of the considered problem. To make the self-consistent analysis possible, we adopt the following two additional assumptions:

* The statistics of a wave field scattered in a weakly nonlinear disordered medium, is close to Gaussian. Consequently, the factorization approximation holds in a weakly nonlinear medium: $`C_{\delta I}(𝐫,\tau )|C_\psi (𝐫,\tau )|^2`$.
* The functional form of the $`𝐫`$-dependence of $`g_1^{(NL)}(𝐫,\tau )`$ is the same as that of $`g_1^{(L)}(𝐫,\tau )`$: $`g_1^{(NL)}(𝐫,\tau )=\mathrm{exp}[\beta (\tau )z/\mathrm{}]`$, where $`\beta (\tau )`$ is some unknown function, which can depend not only on $`\tau `$ but also on other parameters of the problem (namely, on $`\mathrm{}/L_a`$, $`k_0\mathrm{}`$, $`\mathrm{\Delta }n`$).

Strictly speaking, the above assumptions define a sort of perturbation theory, but now we do not limit the values of deviations of intensity and field correlation functions from their values in the linear case. Instead, we assume that the nonlinearity does not cause significant modifications of the statistics of scattered waves \[assumption (iii)\] and of the functional form of the field correlation function \[assumption (iv)\]. Note that now the autocorrelation function $`g_1^{(NL)}(𝐫,\tau )`$ can deviate significantly from $`g_1^{(L)}(𝐫,\tau )`$. Condition (iv) fixes the functional form of this deviation, but implies no constraints on its absolute value.

Obviously, both the assumptions (iii) and (iv) require the nonlinearity to be weak. Assumption (iii) is justified under the same conditions as (i) (see Sec. II), since the Gaussian statistics of the total scattered wave field $`\psi (𝐫,t)`$ is a consequence of the complete randomization of phases of partial waves arriving at $`𝐫`$. The reason for the randomization is that the typical distance $`\mathrm{}`$ between individual scattering events in a multiple-scattering sequence is much larger than the wavelength ($`k\mathrm{}1`$). Obviously, such a mechanism of phase randomization is equally effective in both linear and weakly nonlinear media, as long as (i) holds.

To justify the ansatz of assumption (iv), we apply Eq. (5) and write $`g_1^{(NL)}`$ as

$`g_1^{(NL)}(𝐫,\tau )`$ $``$ $`\left({\displaystyle \frac{3z^2}{4\pi \mathrm{}}}\right)^{1/2}\mathrm{exp}\left({\displaystyle \frac{z}{L_a}}\right)`$ (57)
$`\times `$ $`{\displaystyle _0^{s_1(\tau )}}{\displaystyle \frac{1}{s^{3/2}}}\mathrm{exp}\left({\displaystyle \frac{3z^2}{4\mathrm{}s}}\right)𝑑s,`$ (58)
where we assumed that $`s/\mathrm{}_a+(1/2)\mathrm{\Delta }\phi ^2(\tau )_s`$ increases monotonically with $`s`$ and becomes of order unity at $`s=s_1(\tau )s_0=z^2/(2\mathrm{})`$. After performing integration, Eq. (58) can be approximately rewritten as $`\mathrm{exp}(z/L_a)\{1z[3/(\pi \mathrm{}s_1)]^{1/2}\}\mathrm{exp}[\beta (\tau )z/\mathrm{}]`$ with $`\beta z/\mathrm{}z(\mathrm{}s_1)^{1/2}z/L_a1`$. In the opposite limit of $`\beta z/\mathrm{}1`$, $`g_1^{(NL)}(𝐫,\tau )`$ vanishes and the functional form of its $`z`$-dependence is of no importance. Anyway, it will be seen from the following that the exact functional form of the $`𝐫`$-dependence of $`g_1(𝐫,\tau )`$ is not of crucial importance, since $`g_1`$ is integrated over the whole medium during the calculation of the correlation function of diffusely reflected wave.

Making use of Eqs. (29)–(33) and relying on the assumption (iii), (iv), we recalculate the nonlinear contributions to the dephasing $`\mathrm{\Delta }\phi ^2(\tau )_s`$. Due to the assumption (iv), the result can be obtained from Eqs. (40)–(44) by a simple substitution: $`\alpha (\tau )\beta (\tau )+\alpha (0)`$. Recalling that $`\alpha (0)=\mathrm{}/L_a`$, we obtain:

$`\mathrm{\Delta }\phi ^2(\tau )_s^{(1)}`$ $`=`$ $`2\mathrm{\Delta }n{\displaystyle \frac{\tau }{\tau _0}}\left[1H\left({\displaystyle \frac{\mathrm{}}{L_a}}\sqrt{{\displaystyle \frac{s}{12\mathrm{}}}}\right)\right]{\displaystyle \frac{s}{\mathrm{}}},`$ (59)
$`\mathrm{\Delta }\phi ^2(\tau )_s^{(2)}`$ $`=`$ $`2\pi k_0\mathrm{}\mathrm{\Delta }n^2\{H\left[(\beta +{\displaystyle \frac{\mathrm{}}{L_a}})\sqrt{{\displaystyle \frac{s}{3\mathrm{}}}}\right]`$ (60)
$``$ $`H\left[{\displaystyle \frac{\mathrm{}}{L_a}}\sqrt{{\displaystyle \frac{s}{3\mathrm{}}}}\right]\}{\displaystyle \frac{s}{\mathrm{}}},`$ (61)
$`\mathrm{\Delta }\phi ^2(\tau )_s^{(3)}`$ $``$ $`6\mathrm{\Delta }n^2`$ (62)
$`\times `$ $`\{\begin{array}{cc}\beta \times (s/\mathrm{})^2,\hfill & \\ \beta \sqrt{s/\mathrm{}}1,(\mathrm{}/L_a)\sqrt{s/\mathrm{}}1\hfill & \\ (s/\mathrm{})^{3/2},\hfill & \\ \beta \sqrt{s/\mathrm{}}>1,(\mathrm{}/L_a)\sqrt{s/\mathrm{}}1\hfill & \\ \beta \times (L_a/\mathrm{})^3\times \sqrt{s/\mathrm{}},\hfill & \\ \beta \sqrt{s/\mathrm{}}1,(\mathrm{}/L_a)\sqrt{s/\mathrm{}}>1\hfill & \\ (L_a/\mathrm{})^3,\hfill & \\ \beta \sqrt{s/\mathrm{}}>1,(\mathrm{}/L_a)\sqrt{s/\mathrm{}}>1\hfill & \end{array}`$ (63)
As follows from Eqs. (59)–(63), the nonlinear dephasing depends on the unknown function $`\beta `$. Since the temporal autocorrelation function of diffusely reflected wave $`g_1^{(NL)}(\tau )=\mathrm{exp}[\beta (\tau )]`$ is determined by $`\beta `$ as well, Eqs. (53) and (55) allow us to formulate a self-consistent equation for $`\beta `$:

$`\mathrm{exp}\left[\beta \left(\tau \right)\right]=F\left[\beta \left(\tau \right)\right],`$ (64)
where $`F(\beta )=W\left[\mathrm{\Delta }\phi ^2(\tau )_s\right]/W[0]`$, $`W[\mathrm{}]`$ is defined by Eq. (55), and $`\mathrm{\Delta }\phi ^2(\tau )_s`$ is a sum of terms given by Eqs. (6), (59)–(63). Eq. (64) is the main result of our self-consistent analysis. Although the functional form of $`F(\beta )`$ is rather complicated, and Eq. (64) cannot be solved analytically, the numerical solution is straightforward and can be carried out for any values of $`\mathrm{\Delta }n`$, $`\mathrm{}/L_a`$, $`k_0\mathrm{}`$ and $`\tau `$. Eq. (64) can be considered as a self-consistent equation for the autocorrelation function of diffusely reflected wave, since $`\beta (\tau )`$ and $`g_1^{(NL)}(\tau )`$ are directly related. It is worthwhile to note that in the absence of nonlinearity ($`\mathrm{\Delta }n=0`$), Eq. (64) yields $`\beta (\tau )=\alpha (\tau )\mathrm{}/L_a`$, and hence Eq. (9) is recovered for $`g_1^{(L)}(\tau )`$.

## VII Results and discussion

We start the analysis of Eq. (64) from the case of immobile scatterers, taking a limit of $`\tau /\tau _00`$. We denote the autocorrelation functions corresponding to this limit by $`g_1^{(L),(NL)}(0^+)`$ in order to distinguish them from $`g_1^{(L),(NL)}(0)`$, which correspond to $`\tau =0`$. Obviously, $`g_1^{(L)}(0^+)=g_1^{(L)}(0)=1`$. In a nonlinear medium, Eqs. (61) and (63) can still contribute to the dephasing even for $`\tau /\tau _00`$. These contributions are insensitive to the sign of $`\mathrm{\Delta }n`$. A corresponding value of $`\beta (0^+)`$ and, consequently, of the field autocorrelation function $`g_1^{(NL)}(0^+)`$ can be found by solving Eq. (64) numerically. In Fig. 2, we plot the left-hand and the right-hand sides of Eq. (64) for fixed $`\left|\mathrm{\Delta }n\right|=10^4`$, $`k_0\mathrm{}=100`$, and for several values of $`\mathrm{}/L_a`$. If absorption is strong ($`\mathrm{}/L_a3\times 10^3`$ for considered $`\left|\mathrm{\Delta }n\right|`$ and $`k_0\mathrm{}`$), Eq. (64) has a unique solution $`\beta (0^+)=0`$ which corresponds to $`g_1^{(NL)}(0^+)=1`$. However, a second solution $`\beta (0^+)>0`$ appears for sufficiently weak absorption ($`\mathrm{}/L_a<3\times 10^3`$). The point of appearance of the second solution is a bifurcation point of Eq. (64). To choose the solution realizable in a physical system, we note that the first solution \[$`\beta (0^+)=0`$\] exists only for $`\tau /\tau _0=0`$, and disappears at finite $`\tau /\tau _0`$, since $`F(0)<1`$ for $`\tau /\tau _0>0`$. This solution is therefore inaccessible by continuity, and “unstable” with respect to small scatterer displacements. The physically realizable solution should represent the limit of $`\beta (\tau )`$ for $`\tau /\tau _00`$, which is given by the second solution of Eq. (64). It is therefore this solution which one expects to be realized in a real physical system.

The fact that Eq. (64) can have a positive solution $`\beta (0^+)>0`$ is a very important issue, since $`\beta (0^+)>0`$ leads to $`g_1^{(NL)}(0^+)=\mathrm{exp}[\beta (0^+)]<1`$. A value of the temporal autocorrelation function which is less than unity is commonly associated with temporal fluctuations of scattered waves. In the considered case, however, the reason for these fluctuations is not the motion of scatterers, as the limit of $`\tau /\tau _00`$ corresponds to immobile scatterers. The fluctuations are spontaneous and represent a clear signature of instability of the multiple-scattering speckle pattern.

Despite a rather complicated structure of the function $`F(\beta )`$ in Eq. (64), a relation between the parameters of the problem corresponding to the onset of the speckle pattern instability can be found analytically. As follows from Fig. 2, the initial \[at $`\beta (0^+)=0`$\] decay of $`F[\beta (0^+)]`$ should be faster than the decay of $`\mathrm{exp}[\beta (0^+)]`$, for the second solution of Eq. (64) to appear. A surface in a three-dimensional space of the problem parameters $`\mathrm{\Delta }n`$, $`\mathrm{}/L_a`$, $`k_0\mathrm{}`$, separating the stable \[$`\beta (0^+)=0`$\] and unstable \[$`\beta (0^+)>0`$\] regions, is therefore given by the equation

$`{\displaystyle \frac{}{\beta (0^+)}}F[\beta (0^+)]|_{\beta (0^+)=0}=1.`$ (65)
Recalling that $`L_a/\mathrm{}1`$, $`k_0\mathrm{}1`$ is assumed, we obtain from Eq. (65):

$`p=\mathrm{\Delta }n^2\left({\displaystyle \frac{L_a}{\mathrm{}}}\right)^2\left[k_0\mathrm{}+{\displaystyle \frac{L_a}{\mathrm{}}}\right]1,`$ (66)
where we introduce a control parameter $`p`$, and a numerical factor of order unity is omitted. If $`p<1`$, the multiple-scattering speckle pattern is stable \[$`g_1^{(NL)}(0^+)=1`$\], while for $`p>1`$, an instability shows up leading to $`g_1^{(NL)}(0^+)<1`$. A striking feature of Eq. (66) is that $`p`$ can become larger than unity even for very small $`\left|\mathrm{\Delta }n\right|`$, provided that the extensive parameter $`L_a/\mathrm{}`$ is large enough. Our condition of the speckle pattern stability $`p<1`$ agrees with the condition of validity of the perturbation theory developed in Sec. V \[Eq. (56)\], evaluated at $`\tau /\tau _0=0`$. This readily explains the failure of the perturbation theory for short correlation times and weak absorption: perturbation approach is not suitable for description of unstable regimes. Moreover, the fact that the condition of validity of the perturbation theory and the speckle pattern stability condition agrees indicates that the additional assumptions (iii) and (iv) of Sec. VI are not essential for obtaining unstable regimes.

To illustrate our self-consistent theoretical framework, we solve Eq. (64) numerically for $`\tau /\tau _00`$, and plot the resulting temporal autocorrelation function of diffusely reflected wave, $`g_1^{(NL)}(0^+)`$, in Fig. 3. As discussed above, $`g_1^{(NL)}(0^+)<1`$ corresponds to spontaneous fluctuations of scattered waves which is a manifestation of the speckle pattern instability. It follows from Fig. 3, that in accordance with Eq. (66), an infinitely small $`\left|\mathrm{\Delta }n\right|`$ is sufficient to make the speckle pattern unstable in the absence of absorption, while a certain threshold degree of nonlinearity is required to destabilize the speckle pattern in a dissipative medium. In the absence of absorption, Eq. (64) always has two solutions: $`\beta (0^+)=0`$ and $`\beta (0^+)>0`$, corresponding to $`g_1^{(NL)}(0^+)=1`$ and $`g_1^{(NL)}(0^+)<1`$, respectively. As discussed above, it is the second solution, shown by a dashed line in Fig. 3, which is the physically realizable one. In contrast, if $`\mathrm{}/L_a0`$, $`\left|\mathrm{\Delta }n\right|`$ should be greater than some threshold value for the second solution $`\beta (0^+)>1`$ to appear. Threshold values of $`\left|\mathrm{\Delta }n\right|`$ following from Eq. (66) are shown in Fig. 3 by arrows.

For a nonabsorbing, elastically scattering medium ($`\mathrm{}/L_a=0`$), the value of $`g_1^{(NL)}(0^+)`$ can be estimated analytically. Indeed, at $`\tau /\tau _00`$ the principal contribution to $`\mathrm{\Delta }\phi ^2(\tau )_s`$ is given by $`\mathrm{\Delta }\phi ^2(\tau )_s^{(2)}`$ for $`s/\mathrm{}(k_0\mathrm{})^2`$, and by $`\mathrm{\Delta }\phi ^2(\tau )_s^{(3)}`$ for $`s/\mathrm{}(k_0\mathrm{})^2`$. This allows us to put $`\mathrm{\Delta }\phi ^2(\tau )_s\mathrm{\Delta }\phi ^2(\tau )_s^{(2)}`$ if $`s/\mathrm{}<(k_0\mathrm{})^2`$, and $`\mathrm{\Delta }\phi ^2(\tau )_s\mathrm{\Delta }\phi ^2(\tau )_s^{(3)}`$ if $`s/\mathrm{}>(k_0\mathrm{})^2`$. Integration in Eq. (55) can be then carried out, and Eq. (64) is easily solved, yielding

$`g_1^{(NL)}(0^+)`$ (67)
$`1\{\begin{array}{cc}2\left|\mathrm{\Delta }n\right|^{2/3},\hfill & \left|\mathrm{\Delta }n\right|<(k_0\mathrm{})^{3/2}\text{,}\hfill \\ 3\left|\mathrm{\Delta }n\right|\sqrt{k_0\mathrm{}},\hfill & \left|\mathrm{\Delta }n\right|>(k_0\mathrm{})^{3/2}\text{.}\hfill \end{array}`$ (68)
This result agrees well with the numerical solution of Eq. (64) presented in Fig 3.

The physical origin of the instability of speckle pattern for $`p>1`$ can be revealed by realizing that the system “coherent wave + nonlinear disordered medium” has a positive three-dimensional feedback. In a nonlinear medium, an infinitely small perturbation of the wave intensity $`I(𝐫,t)`$ produces a change of the local refractive index, which alters the phases of waves propagating in the medium and, consequently, affects their mutual interference. Since it is this interference which determines $`I(𝐫,t)`$, the loop of the feedback is closed. For $`p1`$, the feedback is sufficiently strong to compensate for the (diffusive on average) spreading of the initial intensity perturbation, and the speckle pattern $`I(𝐫,t)`$ is unstable. It is worthwhile to note that unstable regimes are not exceptional in nonlinear wave systems and, in particular, in optical systems (see, e.g., Refs. ).

Our Eq. (64) is in no way limited to the case of $`\tau /\tau _00`$, and can be used to compute the temporal autocorrelation function of diffusely reflected wave at any $`\tau /\tau _0>0`$. In the latter case, one should take into account all the four dephasing terms given by Eqs. (6), (59)–(63). The results of the numerical solution of Eq. (64) are shown in Fig. 4. For weak absorption ($`\mathrm{}/L_a=0,10^3`$), $`p>1`$ and the speckle pattern is unstable \[$`g_1^{(NL)}(0^+)<1`$\]. As $`\tau /\tau _0`$ increases, the difference between the “nonlinear” and “linear” curves becomes less pronounced. For strong absorption ($`\mathrm{}/L_a=5\times 10^3`$), $`p`$ becomes less than unity and stability of the speckle pattern is recovered \[$`g_1^{(NL)}(0^+)=1`$\]. We remind that the temporal autocorrelation function $`g_1^{(L)}(\tau )`$, corresponding to a linear medium, always equals $`1`$ for $`\tau /\tau _0=0`$, as shown by dashed lines in Fig. 4.

It is instructive to compare the results obtained using the perturbation theory of Sec. V and the self-consistent approach of Sec. VI. Such a comparison is shown in Fig. 5. The two upper curves, corresponding to relatively strong absorption ($`\mathrm{}/L_a=5\times 10^3`$, $`p<1`$), are almost indistinguishable, which means that for $`p<1`$, the perturbation theory works very well. In contrast, for $`p>1`$ (see the two lower curves of Fig. 5, corresponding to $`\mathrm{}/L_a=0`$), the perturbation and the self-consistent curves are close only at the right from the dotted vertical line, showing the minimum time $`\tau _c`$ at which the perturbation theory is valid \[see Eq. (56)\]. For $`\tau \tau _c`$, the perturbation and the self-consistent results disagree significantly, which confirms our conclusion about the failure of the perturbation approach at short correlation delay times.

## VIII Conclusion

We are now in a position to answer the two central questions formulated in the introductory section:

* The phenomenon of multiple scattering is capable of providing a positive feedback for a coherent wave propagating in a nonlinear disordered medium.
* The onset of the speckle pattern instability occurs when the control parameter

$`p=\mathrm{\Delta }n^2\left({\displaystyle \frac{L_a}{\mathrm{}}}\right)^2\left[k_0\mathrm{}+{\displaystyle \frac{L_a}{\mathrm{}}}\right]`$ (69)
becomes of order or larger than unity. The speckle pattern is stable for $`p<1`$.

The instability of the multiple-scattering speckle pattern manifests itself in spontaneous fluctuations of the scattered wave field and intensity. The following features are characteristic for the development of the instability. First, the development of the instability is independent of the sign of nonlinearity. This is not common for nonlinear waves since the instability is often due to self-focusing phenomena, which only occur for $`n_2>0`$. The instability of waves in a disordered nonlinear medium has nothing to do with the self-focusing, and occurs at relatively weak nonlinearities, when the self-focusing can be neglected. Second, in the absence of absorption, the speckle pattern is unstable for any (even infinitely small) value of the nonlinearity strength $`\left|\mathrm{\Delta }n\right|`$, while in a dissipative medium $`\left|\mathrm{\Delta }n\right|`$ should exceed a certain threshold value for the instability to show up. Finally, the instability results in a value of the autocorrelation function $`g_1^{(NL)}(\tau )`$ of scattered wave, which is smaller than $`1`$ for $`\tau /\tau _00`$ (i.e., in the absence of scatterer motion). This is a clear signature of spontaneous fluctuations of the multiple-scattered speckle pattern, and should be observable in experiments.

###### Acknowledgements.

The author is indebted to R. Maynard for numerous discussions and continuous support of this work.

## A Spatio-temporal filed correlation in a linear medium

In this appendix, we provide a derivation of the spatio-temporal correlation function, $`C_\psi (𝐫,\mathrm{\Delta }𝐫,\tau )=\psi (𝐫\mathrm{\Delta }𝐫/2,t)\psi ^{}(𝐫+\mathrm{\Delta }𝐫/2,t+\tau )`$, of a random field $`\psi (𝐫,t)`$ in the bulk of disordered medium. Starting from the linear wave equation \[Eq. (2) with $`\epsilon _2=0`$\], we obtain the Bethe-Salpeter equation in the form:

$`C_\psi (𝐫,\mathrm{\Delta }𝐫,\tau )=\psi (𝐫\mathrm{\Delta }𝐫/2,t)\psi ^{}(𝐫+\mathrm{\Delta }𝐫/2,t)`$ (A1)
$`+{\displaystyle 𝑑𝐫_a𝑑𝐫_b𝑑𝐫_c𝑑𝐫_d\overline{G}(𝐫\mathrm{\Delta }𝐫/2,𝐫_a)}`$ (A2)
$`\times \overline{G}^{}(𝐫+\mathrm{\Delta }𝐫/2,𝐫_b)U(𝐫_a,𝐫_b,𝐫_c,𝐫_d)`$ (A3)
$`\times C_\psi [(𝐫_c+𝐫_d)/2,𝐫_d𝐫_c,\tau ],`$ (A4)
where the integrations are over the volume of disordered medium, $`\overline{G}`$ is the average Green function of the linear wave equation, and $`U`$ is the irreducible four-point vertex. Far enough from the medium boundaries, we can replace $`\overline{G}`$ by its value in the infinite medium:

$`\overline{G}(𝐫_1,𝐫_2)`$ (A5)
$`={\displaystyle \frac{1}{4\pi \left|𝐫_1𝐫_2\right|}}\mathrm{exp}\left[\left(ik{\displaystyle \frac{1}{2\mathrm{}}}\right)\left|𝐫_1𝐫_2\right|\right].`$ (A6)
Now we assume that point scatterers in the medium undergo Brownian motion, and that the correlation function of the dielectric function fluctuations is given by Eq. (4). For $`\tau \tau _0`$, we can neglect the $`\tau `$-dependence of $`U`$ in Eq. (A4), which in the limit of $`k\mathrm{}1`$ becomes

$`U(𝐫_a,𝐫_b,𝐫_c,𝐫_d)`$ (A7)
$`={\displaystyle \frac{4\pi }{\mathrm{}}}\delta (𝐫_a𝐫_c)\delta (𝐫_b𝐫_d)\delta (𝐫_a𝐫_b).`$ (A8)
Eq. (A4) then reduces to

$`C_\psi (𝐫,\mathrm{\Delta }𝐫,\tau )={\displaystyle \frac{4\pi }{\mathrm{}}}{\displaystyle 𝑑𝐫_a\overline{G}(𝐫\mathrm{\Delta }𝐫/2,𝐫_a)}`$ (A9)
$`\times \overline{G}^{}(𝐫+\mathrm{\Delta }𝐫/2,𝐫_a)C_\psi (𝐫_a,0,\tau ),`$ (A10)
where we assumed the coherent field $`\psi (𝐫,t)`$ to be negligible.

Now we remind that $`C_\psi (𝐫_a,0,\tau )C_\psi (𝐫_a,\tau )`$ is a slow-varying function of $`𝐫_a`$. We can therefore pull it out of the integral in Eq. (A10), taking its value at $`𝐫`$. Performing the remaining integral, we find

$`C_\psi (𝐫,\mathrm{\Delta }𝐫,\tau )={\displaystyle \frac{\mathrm{sin}(k\mathrm{\Delta }r)}{k\mathrm{\Delta }r}}\mathrm{exp}\left({\displaystyle \frac{\mathrm{\Delta }r}{2\mathrm{}}}\right)C_\psi (𝐫,\tau ).`$ (A11)
This equation relates the spatio-temporal correlation of the field with its purely temporal correlation. Eq. (A11) holds for any sample geometry and source distribution, far enough from boundaries and sources, and for $`\tau \tau _0`$. In the case of a plane wave incident at the surface $`z=0`$ of a semi-infinite disordered medium, occupying the $`z>0`$ half-space,

$`C_\psi (𝐫,\tau )=I_0\mathrm{exp}\left[\alpha (\tau )(z/\mathrm{})\right]`$ (A12)
with $`\alpha ^2(\tau )=3\tau /(2\tau _0)+(\mathrm{}/L_a)^2`$, leading to Eq. (11) in the main text. In the factorization approximation (12), the square of Eq. (A11) gives the short-range correlation function of intensity fluctuations.

## B Long-range spatio-temporal correlation of intensity fluctuations in a linear medium

To calculate the spatio-temporal correlation function of intensity fluctuations, $`C_{\delta I}(𝐫,\mathrm{\Delta }𝐫,\tau )=\delta I(𝐫\mathrm{\Delta }𝐫/2,t)\delta I(𝐫+\mathrm{\Delta }𝐫/2,t+\tau )`$, for $`\mathrm{\Delta }r>\mathrm{}`$ we generalize the Langevin approach, which has been initially developed for calculation of purely spatial correlations. We assume that the spatio-temporal correlation function of Langevin random sources is determined by the short-range correlation function of intensity fluctuations. In the factorization approximation, the latter is given by the square of the field-field correlation function \[Eq. (A11)\]. Hence, the generalized Langevin equations read

$`D_p\left[^21/L_a^2\right]\delta I(𝐫,t)=\mathrm{div}𝐣_{ext}(𝐫,t),`$ (B1)
$`j_{ext}^{(i)}(𝐫\mathrm{\Delta }𝐫/2,t)j_{ext}^{(m)}(𝐫+\mathrm{\Delta }𝐫/2,t+\tau )`$ (B2)
$`={\displaystyle \frac{1}{3}}\delta _{im}{\displaystyle \frac{c^22\pi \mathrm{}}{k^2}}\left|C_\psi (𝐫,\tau )\right|^2\delta (\mathrm{\Delta }𝐫),`$ (B3)
where $`D_p=c\mathrm{}/3`$ is the diffusion coefficient, and $`c`$ is the speed of wave in the medium. In the case of a plane wave incident upon a boundary $`z=0`$ of a semi-infinite disordered medium, it is convenient to make a two-dimensional Fourier transform of Eq. (B1) in the $`\{x,y\}`$ plane. Then the equations corresponding to $`\delta I(𝐊_1,z_1,t_1)`$ and $`\delta I(𝐊_2,z_2,t_2)`$ are multiplied, and the product is ensemble averaged. Further transformations, which are equivalent to those discussed in Ref. , yield

$`\delta I(𝐊_1,z_1,t_1)\delta I(𝐊_2,z_2,t_2)=I_0^2{\displaystyle \frac{6\pi }{k^2\mathrm{}}}\delta (𝐊_1𝐊_2)`$ (B4)
$`{\displaystyle _0^{\mathrm{}}}dz^{}[(𝐊_\mathrm{𝟏}𝐊_2)G(p,z_1,z^{})G(p,z_2,z^{})`$ (B5)
$`+{\displaystyle \frac{}{z^{}}}G(p,z_1,z^{}){\displaystyle \frac{}{z^{}}}G(p,z_2,z^{})\left]\mathrm{exp}\right[2\alpha (\tau ){\displaystyle \frac{z^{}}{\mathrm{}}}],`$ (B6)
where $`p^2=𝐊^2+1/L_a^2`$, and

$`G(p,z_1,z^{})`$ $`=`$ $`{\displaystyle \frac{1}{p}}\mathrm{sinh}(p\times \mathrm{min}\{z_1,z^{}\})`$ (B7)
$`\times `$ $`\mathrm{exp}(p\times \mathrm{max}\{z_1,z^{}\})`$ (B8)
is the Green function of Eq. (B1). Evaluating Eq. (B6) and transforming the result back to the real space, after lengthy but straightforward algebra we obtain

$`C_{\delta I}(𝐫,\mathrm{\Delta }𝐫,\tau )`$ $`=`$ $`{\displaystyle \frac{3}{(k\mathrm{})^2}}I_0^2{\displaystyle _0^{\mathrm{}}}𝑑KK`$ (B9)
$`\times `$ $`Q[K,\sqrt{K^2+\mathrm{}^2/L_a^2},{\displaystyle \frac{z}{\mathrm{}}},{\displaystyle \frac{\mathrm{\Delta }z}{\mathrm{}}},\alpha (\tau )]`$ (B10)
$`\times `$ $`J_0\left(K{\displaystyle \frac{\mathrm{\Delta }R}{\mathrm{}}}\right),`$ (B11)
where

$`Q(K,p,\zeta ,\mathrm{\Delta }\zeta ,\alpha )=\mathrm{exp}\left(2p\zeta \right)\times [{\displaystyle \frac{2\alpha ^2+K^2p^2}{4\alpha ^34\alpha p^2}}`$ (B12)
$`[(\alpha ^2p^2)(K^2+p^2)+\alpha ^2(K^2+p^2)\mathrm{cosh}(2p\zeta _1)`$ (B13)
$`+\alpha p(K^2+p^2)\mathrm{sinh}(2p\zeta _1)]/\left[4\alpha \mathrm{exp}(2\alpha \zeta _1)(\alpha ^2p^2)p^2\right]]`$ (B14)
$`+\left(1+{\displaystyle \frac{K^2}{p^2}}\right){\displaystyle \frac{\mathrm{sinh}(p\zeta _1)\mathrm{sinh}(p\zeta _2)}{2\mathrm{exp}\left[2(\alpha +p)\zeta _2\right](\alpha +p)}}`$ (B15)
$`+{\displaystyle \frac{\mathrm{exp}(p\zeta _2)\mathrm{sinh}(p\zeta _1)}{4\alpha p^2(\alpha +p)}}`$ (B16)
$`\times \{p(p^2+2\alpha pK^2)[\mathrm{cosh}(p\zeta _2)\mathrm{exp}[(2\alpha +p)\zeta _2]`$ (B17)
$`\mathrm{cosh}(p\zeta _1)\mathrm{exp}[(2\alpha +p)\zeta _1]]`$ (B18)
$`+(p^32\alpha K^2pK^2)[\mathrm{sinh}(p\zeta _2)\mathrm{exp}[(2\alpha +p)\zeta _2]`$ (B19)
$`\mathrm{sinh}(p\zeta _1)\mathrm{exp}[(2\alpha +p)\zeta _1]]\},`$ (B20)
where $`\zeta _1=\zeta \mathrm{\Delta }\zeta /2`$ and $`\zeta _2=\zeta +\mathrm{\Delta }\zeta /2`$.

Equations (B11), (B20) cannot be evaluated in a general form. In the limits of $`\alpha (\tau )\zeta ,\alpha (\tau )\mathrm{\Delta }R/\mathrm{}1`$ we can, however, approximately replace $`Q`$ in Eq. (B11) by

$`Q(K,p,\zeta ,\mathrm{\Delta }\zeta ,\alpha )`$ $``$ $`{\displaystyle \frac{1}{2K}}\left[\mathrm{exp}(K\mathrm{\Delta }\zeta )\mathrm{exp}(2K\zeta )\right]`$ (B21)
$`\times `$ $`\mathrm{exp}[2\alpha (\tau )\zeta ],`$ (B22)
which yields Eq. (23) in the main text.

## C Density of wave paths in a disordered medium

Propagation of waves in a disordered medium can be interpreted in terms of partial waves traveling along various paths inside the medium. The spatial distribution of such paths and their relative weights depend on the scattering properties of the medium, and on the geometry of the sample. In the case of multiple scattering, the simplest and, at the same time, sufficiently accurate model of wave propagation is the diffusion model. According to this model, wave paths in the medium coincide with trajectories of Brownian particles. The probability $`G(𝐫_1,𝐫_2,s)`$ for a path of length $`s`$ to pass from $`𝐫_1`$ to $`𝐫_2`$ is then given by a solution of the diffusion equation, which in the absence of absorption reads:

$`{\displaystyle \frac{d}{ds}}G(𝐫_1,𝐫_2,s){\displaystyle \frac{1}{3}}\mathrm{}^2G(𝐫_1,𝐫_2,s)=\delta (𝐫_1𝐫_2)\delta (s),`$ (C1)
where $`\mathrm{}`$ is the mean free path. Commonly used boundary conditions for Eq. (C1) consist in putting $`G=0`$ at open boundaries and $`_nG=0`$ at reflecting boundaries of the sample (where $`_n`$ denotes the normal derivative of $`G`$). $`G(𝐫_1,𝐫_2,s)`$ is called the Green function, or the propagator. For a semi-infinite medium occupying the half-space $`z>0`$, one finds

$`G(𝐫_1,𝐫_2,s)`$ $`=`$ $`\left({\displaystyle \frac{3}{4\pi \mathrm{}s}}\right)^{3/2}\{\mathrm{exp}[{\displaystyle \frac{3}{4\mathrm{}s}}(\mathrm{\Delta }R^2+\mathrm{\Delta }z^2)]`$ (C2)
$``$ $`\mathrm{exp}[{\displaystyle \frac{3}{4\mathrm{}s}}(\mathrm{\Delta }R^2+Z^2)]\},`$ (C3)
where cylindrical coordinates are used: $`𝐫_i=\{𝐑_i,z_i\}`$, $`\mathrm{\Delta }𝐑=𝐑_1𝐑_2`$, $`\mathrm{\Delta }z=z_1z_2`$, and $`Z=z_1+z_2`$.

Following Ref. , we introduce $`\rho _s(𝐫_1,𝐫_2,𝐫_3)`$, the density distribution of paths of length $`s`$, as a number of visits of a given site $`𝐫_2`$ inside $`d^3𝐫_2`$ in the ensemble of paths of length $`s`$ starting at $`𝐫_1`$ and ending at $`𝐫_3`$, over the total length of the ensemble distinct paths:

$`\rho _s(𝐫_1,𝐫_2,𝐫_3)`$ $`=`$ $`{\displaystyle \frac{1}{sG(𝐫_1,𝐫_3,s)}}{\displaystyle _0^s}𝑑pG(𝐫_1,𝐫_\mathrm{𝟐},p)`$ (C4)
$`\times `$ $`G(𝐫_2,𝐫_3,sp).`$ (C5)
$`\rho _s(𝐫_1,𝐫_2,𝐫_3)`$ describes the probability density for a path of a given length $`s`$, starting at $`𝐫_1`$ and ending at $`𝐫_3`$, to pass through $`𝐫_2`$. This quantity is normalized:

$`{\displaystyle d^3𝐫_2\rho _s(𝐫_1,𝐫_2,𝐫_3)}=1,`$ (C6)
where the integration is performed over the volume of disordered medium.

As the Green function $`G`$ is known \[Eq. (C3)\], the calculation of $`\rho _s(𝐫_1,𝐫_2,𝐫_3)`$ is straightforward. For diffusely reflected paths, assuming that the first and the last scattering events take place at $`z=\mathrm{}`$, we obtain

$`\rho _s(\mathrm{},𝐫,\mathrm{})\rho _s(𝐫)={\displaystyle \frac{1}{A}}{\displaystyle \frac{6z}{\mathrm{}s}}\mathrm{exp}\left({\displaystyle \frac{3z^2}{\mathrm{}s}}\right),`$ (C7)
where $`A\mathrm{}`$ is the surface of the semi-infinite medium. Equation (C7) defines the probability density for a diffusely reflected path of length $`s`$ to pass through a vicinity of some point $`𝐫=\{x,y,z\}`$.

Generalizing definition of $`\rho _s`$, we define the probability density for a path of length $`s`$ starting at $`𝐫_1`$ and ending at $`𝐫_4`$ to pass consequently through $`𝐫_2`$ and $`𝐫_3`$:

$`\rho _s(𝐫_1,𝐫_2,𝐫_3,𝐫_4)`$ $`=`$ $`{\displaystyle \frac{2}{s^2G(𝐫_1,𝐫_4,s)}}{\displaystyle _0^s}𝑑p{\displaystyle _0^{sp}}𝑑q`$ (C9)
$`G(𝐫_1,𝐫_\mathrm{𝟐},p)G(𝐫_2,𝐫_3,q)`$
$`\times `$ $`G(𝐫_3,𝐫_4,spq).`$ (C10)
The normalization of Eq. (C10) is

$`{\displaystyle d^3𝐫_2d^3𝐫_3\rho _s(𝐫_1,𝐫_2,𝐫_3,𝐫_4)}=1.`$ (C11)
For a semi-infinite medium, we get

$`\rho _s(\mathrm{},𝐫,𝐫^{},\mathrm{})`$ $``$ $`\rho _s(𝐫,𝐫^{})={\displaystyle \frac{1}{A^2}}{\displaystyle \frac{9}{2\pi \mathrm{}^2s^2}}`$ (C12)
$`\times `$ $`\{{\displaystyle \frac{Z+\sqrt{\mathrm{\Delta }R^2+\mathrm{\Delta }z^2}}{\sqrt{\mathrm{\Delta }R^2+\mathrm{\Delta }z^2}}}`$ (C13)
$`\times `$ $`\mathrm{exp}\left[{\displaystyle \frac{3}{4\mathrm{}s}}\left(Z+\sqrt{\mathrm{\Delta }R^2+\mathrm{\Delta }z^2}\right)^2\right]`$ (C14)
$``$ $`{\displaystyle \frac{Z+\sqrt{\mathrm{\Delta }R^2+Z^2}}{\sqrt{\mathrm{\Delta }R^2+Z^2}}}`$ (C15)
$`\times `$ $`\mathrm{exp}[{\displaystyle \frac{3}{4\mathrm{}s}}(Z+\sqrt{\mathrm{\Delta }R^2+Z^2})^2]\},`$ (C16)
where $`𝐫=\{𝐑,z\}`$, $`\mathrm{\Delta }𝐑=𝐑𝐑^{}`$, $`\mathrm{\Delta }z=zz^{}`$, $`Z=z+z^{}`$.

Although Eqs. (C7) and (C16) have been found for a nonabsorbing medium, it is easy to show that these results hold in the presence of spatially-homogeneous absorption as well. This stems from the fact that the attenuation of wave in a homogeneously absorbing medium depends only on the path length, while being independent of the path geometry. As a consequence, the Green function \[Eq. (C3)\] should be multiplied by a factor $`\mathrm{exp}(s/\mathrm{}_a)`$, where $`\mathrm{}_a`$ is the absorption length. This factor, however, disappears after the substitution of the Green function (C3) in Eqs. (C5) and (C10). Consequently, $`\rho _s(𝐫)`$ and $`\rho _s(𝐫,𝐫^{})`$ are independent of $`\mathrm{}_a`$ and remain unchanged.