Acoustic wave propagation through solar granulation: Validity of effective-medium theories, coda waves

Context. The frequencies, lifetimes, and eigenfunctions of solar acoustic waves are affected by turbulent convection, which is random in space and in time. Since the correlation time of solar granulation and the periods of acoustic waves ($\sim$5 min) are similar, the medium in which the waves propagate cannot a priori be assumed to be time independent. Aims. We compare various effective-medium solutions with numerical solutions in order to identify the approximations that can be used in helioseismology. For the sake of simplicity, the medium is one dimensional. Methods. We consider the Keller approximation, the second-order Born approximation, and spatial homogenization to obtain theoretical values for the effective wave speed and attenuation (averaged over the realizations of the medium). Numerically, we computed the first and second statistical moments of the wave field over many thousands of realizations of the medium (finite-amplitude sound-speed perturbations are limited to a 30 Mm band and have a zero mean). Results. The effective wave speed is reduced for both the theories and the simulations. The attenuation of the coherent wave field and the wave speed are best described by the Keller theory. The numerical simulations reveal the presence of coda waves, trailing the coherent wave packet. These late arrival waves are due to multiple scattering and are easily seen in the second moment of the wave field. Conclusions. We find that the effective wave speed can be calculated, numerically and theoretically, using a single snapshot of the random medium (frozen medium); however, the attenuation is underestimated in the frozen medium compared to the time-dependent medium. Multiple scattering cannot be ignored when modeling acoustic wave propagation through solar granulation.


Introduction
Solar seismic waves interact with small-scale convective motions near the solar surface via a wave scattering process, affecting their properties (e.g., propagation speed, frequency, amplitude, and phase). As the e-folding lifetime of solar granulation is comparable to the period of the waves, the medium may not be assumed to be frozen. Furthermore, the spatial spectrum of convection encompasses all scales, including those that are comparable to the wavelengths of p and f modes.
Most approaches that have been proposed so far assume a separation of scales between the waves and the medium. Often the wave period is assumed to be much smaller than the time scale of the evolution of convective flows (Brown 1984;Delache & Fossat 1988;Rosenthal et al. 1999). Murawski & Roberts (1993a,b), using this assumption, derived a model for the scattering of the f mode by granulation using the binary collision approximation (Howe 1971) and found a mode frequency reduction due to the scattering as well as a large attenuation, which compare favorably with observations (Duvall et al. 1998). Other authors assume that the wave period, or the wavelength, is much larger than the temporal, or the spatial, scale of convection, which allows one to apply homogenization techniques (Hanasoge et al. 2013;Bhattacharya et al. 2015).
Numerical simulations provide a useful means to study the interaction of seismic waves with convection (e.g., Ball et al. 2016;Houdek et al. 2017;Schou & Birch 2020). Turbulent convection has an indirect effect on the waves through a change in the average medium (e.g., via a turbulent pressure term) and, in addition, it affects the physics of wave propagation and attenuation via a scattering process.
Here, we study this problem under a highly simplified setup. We consider a one-dimensional steady medium that contains sound-speed perturbations over a finite region, but other than that it is uniform. There is no a priori separation of scales in space nor in time between the incoming wave packet and the medium. For relative sound-speed perturbations of a significant amplitude (e.g., 5% and above) multiple scattering plays a significant role in the redistribution of wave energy. We compare our numerical simulations with theoretical approximations, which are easy to implement in this context.
Since the medium is random in both time and space, we study the effect of the medium on the waves in a statistical sense by computing the first and second moments of the quantities of interest (e.g., the wave field) over many realizations. From the expectation value of the wave field, also known as the coherent or ballistic wave field, we can extract the attenuation and the effective wave speed for example. For the variance of the wave field, we can extract information about the distribution of backward-and forward-scattered energy. This includes latearrival fluctuations due to multiply-scattered (coda) waves.
Article number, page 1 of 13 arXiv:2010.01174v2 [astro-ph.SR] 13 Oct 2020 A&A proofs: manuscript no. aanda We state the problem in Section 2 and explain the numerical implementation in Section 3. Various effective medium theories are reviewed in Section 4. We present our results in Section 5, and discuss them in Section 6.

The random medium
We consider a uniform one-dimensional background with sound speed c 0 = 10 km/s, a value of the same order of magnitude as the sound speed at the solar surface. We perturb the medium by adding locally a space-and time-dependent random fluctuation: This is shown in Fig. 1, where the filled circles symbolize the fluctuation. In Eq. (1), δc has a zero mean so that c = c 0 where angle brackets denote an expectation value.
- The sound-speed perturbation is specified through the autocorrelation where we assume a separation between time and space. The value of is at most 0.1 in our simulations. The random medium can equivalently be characterized by its power spectrum In time, we choose an exponential profile where τ is the e-folding lifetime. For granulation, we have τ ≈ 400 s (e.g., Title et al. 1989). The temporal power spectrum is Lorentzian, In space, we consider two different types of profile. The first choice is an exponential medium (hereafter Medium 1), which will enable us to carry out approximations analytically: For a granulation-like medium, it is reasonable to choose a = 1 Mm. In Fourier space, The second choice (hereafter Medium 2) is a spatial power spectrum of the form (e.g., Baran 2013) where C = πβ α+1 /Γ(α + 1) is a normalization factor such that the spatial autocorrelation function equals 1 at x = 0.The parameters α and β can be tuned to obtain a power spectrum that peaks at the desired wavenumber. Here we fix α = 1 and β = 6.7 × 10 −4 R where R = 696 Mm, such that the spatial power spectrum peaks at kR = 1500. In real space, for α = 1, we have The two power spectra and their corresponding autocorrelation functions are shown in Fig. 2. From the knowledge of the power spectrum P(k, ω), we can compute a realization of the sound speed perturbations as follows: where N(ω, k) is realization of a complex Gaussian random variable with zero mean and unit variance (the real and the imaginary parts are independent). To ensure that δc(x, t) is real, we have N(k, ω) = N * (−k, −ω). This way to proceed is based on the assumptions of stationarity and horizontal spatial homogeneity of the medium (e.g., Gizon & Birch 2004).

The wave equation
The displacement ξ of acoustic waves is given by (Lynden-Bell & Ostriker 1967) Here, we have ignored gravity, rotation, damping as well as any background flows. This equation has been derived in a background medium where the parameters ρ and c are independent of time. However, Legendre (2003) showed that this formulation remains valid for a time-varying medium. Taking the divergence of Eq. (11) and denoting φ = ∇ · ξ, we obtain In this paper, we assume that the density is constant and consider the following 1D acoustic wave equation We implement two numerical codes. The first code is a timedomain code to study the propagation of the wavepacket through a time-dependent random medium, based on Eq. (13). As initial condition, we inject at location x 0 < X a wave packet of central frequency ω 0 and frequency width σ: where As shown in the schematics of Fig. 1, the incoming wave packet first travels in the +x direction in the homogeneous medium, experiences scattering inside the perturbed medium, then comes out (outgoing wave packet) and propagates in the +x direction in the homogeneous medium. Part of the wave packet is backscattered and travels in the −x direction. The simulation box is large enough so that the wave packet is not affected by the computational boundaries at x = 0 and x = x max . The second code is a frequency-domain code to study the wave field in a frozen medium (τ → ∞). For a sound speed that does not depend on time, we can take the temporal Fourier transform of Eq. (13) to obtain the wave equation in the frequency domain, i.e. the Helmholtz equation with Dirichlet boundary condition at x = 0, and the Sommerfeld outgoing radiation condition The tilde denotes the temporal Fourier transform.

Characterizing the wave field
The wave field is affected randomly by the perturbations. The statistical effects can however be studied by looking at the moments of the wave field, i.e. by doing some averages over the realizations of the random medium.
In particular, the coherent wave field is attenuated because each wave packet travels in a different random realization of the medium and is deformed in a different way. This damping is related to the lifetime of the average acoustic wave. The coherent wave field also propagates with a different velocity than c 0 , depending on frequency, called the effective wave speed.
An approximate representation of the coherent wave field inside the perturbed medium is therefore where the effective wave number is and the spatial attenuation is The effective wave speed is defined by Similarly, we define k 0 (ω), the wave number for an unperturbed wave field, such that We want to solve the (simplified) problem of acoustic wave scattering numerically and find which approximations work to retrieve the coherent wave field. In particular, we check whether we can get rid of the time dependence and assume a frozen medium. Furthermore, we want to investigate the phenomenon of multiple scattering due to the finite-amplitude perturbations. This is more easily done by looking at the second moment (variance) of the wave field. It contains information that is otherwise zeroed out by doing a mere average. By doing so, the coda waves, which trail the ballistic wave packet and are often studied in seismology, can be readily observed.

Numerical scheme to solve for φ(x, t)
In order to solve numerically Eq. (13), we use an explicit finitedifference scheme of second order. We choose ω 0 /2π = 3 mHz and σ/2π = 1 mHz so that the frequency range of study is 1 to 5 mHz, which is a reasonable choice for solar acoustic waves. The wave packet is initially at x 0 = 100 Mm, while x max = 200 Mm and t max = 10000 s. We set X = 120 Mm and L = 30 Mm. The resolutions for the simulations are ∆x = 50 km and ∆t = 2.5 s, so that c 0 ∆t/∆x = 0.5 < 1.
An example of time-domain simulation with medium 1 is shown in the online movie and in Fig. 3. The wave packet begins to be perturbed when it enters the random medium. Most of the signal is transmitted forward roughly in the form of a wave packet (ballistic wave packet). Small oscillations trail that signal, propagating either forward or backward. Once out of the perturbation, the shape of the wave packet is not modified anymore.

Numerical scheme to solve forφ(x, ω)
The code uses a second-order discretization scheme with a spatial resolution δx = 4 km. A tridiagonal system is inverted with the tridiagonal matrix algorithm (Thomas algorithm). c 0 , x max , X and L are the same as for the time-domain code. The resolution is done for frequencies between 1 and 5 mHz.

Measuring the attenuation
Following Aki & Richards (2002), after propagating between two points x 1 and x 2 (x 2 > x 1 ) in an attenuating medium, a plane wave is damped by a factor The spatial attenuation k i could be measured from the amplitude difference between the incoming and the outgoing wave packets. However, this method leads to artifacts due to boundary effects occurring at the edges of the random medium. Therefore, we rather consider the wave packet inside the perturbation. We take x 1 = 126 Mm and x 2 = 144 Mm, each point being 6 Mm away from the edge of the perturbation. As shown in Fig. 4 The power of the signal has been attenuated during the propagation from x 1 to x 2 . At each frequency, we then fit a first order polynomial to the natural logarithm of the norm of the Fourier component in order to retrieve the decay coefficient. We fit its slope between the vertical dotted lines, corresponding to [x 1 , x 2 ]. The vertical dashed lines delimit the location of the perturbation in time (for a wave packet propagating at c 0 ) and in space. In the figure, the Fourier components have been normalized so that they have the same amplitude before entering the perturbation.

Measuring the effective wave speed
We first apply a temporal Fourier transform to φ . Then we fit to Re( φ (x, ω) ) (we could have chosen the imaginary part arbitrarily), in the perturbed region, an exponentially decreasing oscillatory function where the decay rate has been determined via the method to measure the attenuation from the previous section. More precisely, we fit where A, x s , and c eff are the free parameters, with x s being a phase shift. A similar fit is done on φ 0 to take numerical dispersion into account. This is illustrated in Fig. 5.

Effective medium theories
Depending on the values of the parameters, in particular k 0 a and , different theories can be used to compute the effective parameters c eff and k i . For a medium whose spatial scale is much less Table 1. Theories used in this paper for the effective wave speed and attenuation in a frozen medium as → 0, and their range of validity. The Keller and the Born theories are for medium 1.

Theory
Validity range k i c eff Keller 1964 1 Notes. ( †) Approximation for k 0 a > 1 of Rytov et al. (1989b) who made the derivation for a Gaussian autocorrelation function and single scattering. ( ‡) The dominant term is that of Keller for small perturbations (see for example Fig. F.1). # See Rytov et al. (1989a). than the wavelength (λ a, regime of Rayleigh scattering), the homogenization method is appropriate and gives an effective sound speed c eff = c 0 (1 − 3/2 2 ) (derivation for a frozen medium in Appendix D). On the other hand, for small wavelengths (λ a, small-angle scattering regime), the geometrical optics approach is relevant and implies c eff = c ray = c 0 (1 − 2 ) (derivation for a frozen medium in Appendix E). These two approaches give an effective wave speed which is independent of frequency and of the power spectrum of the perturbation but do not provide any attenuation. Another caveat is that in our model, we are in the regime of Mie scattering or large-angle scattering (Aki & Wu 1988) because λ 3.3a. The wave number can therefore not be considered large nor small compared to 1 and other theories may be required.
We explore two other derivations in the case of small perturbations ( 1). In this regime, two methods are used: the Keller solution, derived for a frozen (Appendix A) or time-dependent medium (Appendix B); the Born second-order approximation for a frozen medium (derived in Appendix C). The Keller and the Born solutions converge toward the same values for small perturbations ( = 0.01 for instance). However, at = 0.1, the Born solution is very different from the Keller one, and it is not possible to fit a function of the form of Eq.(19). We come back to this point in Section 5.1.
If k 0 a 1 or k 0 a 1, the effective wave speed obtained from the Keller theory converges toward the results of the ho-mogenization and geometrical optics approaches, respectively. Table 1 summarizes the expressions for a frozen medium 1 (for which the analytical expressions can be easily derived) as → 0. We check the agreement of these theories with our simulations in Appendix F. The Keller and the Born second-order theories are consistent and predict that k i is proportional to k 2 0 while the effective wave speed difference is essentially independent of k 0 provided that k 0 a 1. For comparison purposes, we note that Bourret (1963) and Sato et al. (2012) (pp. 214-220) found that for a 3D frozen medium with an exponential autocorrelation (i.e., similar to medium 1), the attenuation is proportional to k 4 0 for k 0 a 1 and to k 2 0 for k 0 a 1 while the behavior of the effective wave speed is in qualitative agreement. In particular, the effective wave speed is always less than c 0 . On the other hand, van der Baan (2001) used the wave localization theory to make use of so-called self-averaging quantities; he derived the effective medium using one realization of a one-dimensional perturbation in density and bulk modulus. He found that the attenuation coefficient tends toward a constant value at high k 0 a, while the effective wave speed difference is positive, in agreement with previous studies (e.g., Müller et al. 1992). The sign of the difference c eff − c 0 and the dependence of the attenuation on frequency therefore seem to depend strongly on the equation that is solved.

Results and comparisons
We compute the properties of the effective medium using the procedure explained in Sections 3.3 and 3.4, with = 0.1 and a = 1 Mm. The 1-σ error bars shown later on on the attenuation and wave speed measurements are obtained from ten sets of 10 4 realizations.

Coherent wave field
We first reconstruct the theoretical coherent wave field obtained when using the various theories. For all of them (except the Born theory that provides directly the wave field), we assume the form written in Eq. (19). In Fig. 6, we plot Re( φ (3 mHz, x) ) (x ∈ [X, X+L]) for a frozen medium 1. Clearly the Keller approximation does the best job at approximating the true (numerical) solution. As mentioned before, neither the homogenization technique nor the geometrical optics makes an attenuation emerge. The Born solution is a good approximation on about half of the random medium at this frequency, before it becomes out of phase with the numerical solution while its amplitude also starts to increase. The discrepancy is worse and arises earlier in the medium for higher frequencies. Thus the Born approximation, although similar to the Keller approximation when → 0, performs poorly for a 10% perturbation in a medium of size L = 30 Mm.  Fig. 7. Attenuation of the coherent wave packet vs frequency for media 1 and 2, after propagation through a band of perturbed medium. The 1D theory from Keller is overplotted in dashed lines. 1-σ error bars are shown. Fig. 7 shows the measured attenuation for simulations with τ = 400 s and τ → ∞. The case τ = 1 day (not shown on the plot) lies within the error bars of the curve for the frequency code, which is to be expected as the typical time scales involved (the period of the wave, about 5 minutes, and the time it takes for it to travel through the medium, about 1 h) are much less than one day. We superimpose the attenuation that one expects from the time-dependent Keller theory.

Attenuation
The attenuation by medium 1 is an increasing function of frequency, with a value of about 1.5% of the wave number at 3 mHz for τ = 400 s. The ratio k i /k 0 is a linear function of frequency, meaning that k i is quadratic, as expected from the Keller theory. For medium 2 however, the attenuation is not quadratic. It reaches about 0.5% of the wave number at 3 mHz, which is smaller than the medium 1 value by a factor 3. The smaller attenuation values are caused by the lack of power toward low wave numbers in the spectrum of the perturbation: the absence of large scales in the perturbation means that the incoherence between the realizations of the wave packets occurs preferentially on small scales, thereby decreasing the overall broadening of the wave packet. In medium 2, the ratio k i /k 0 stabilizes above 5 mHz for τ = 400 s, while it reaches a maximum at about 3 mHz for τ → ∞. This may indicate that there is a preferred scale of damping of the coherent wave field. Fig. 8 shows the effective wave speed computed with for τ = 400 s and τ → ∞, as well as the time-dependent Keller theory, the (frozen) spatial homogenization solution and the (frozen) geometrical optics solution. Like for the attenuation, the case τ = 1 day (not shown) lies within the error bars of the curve for the frequency code. The effective wave speed is less than the unperturbed sound speed c 0 . This is due in part to waves being scattered back and forth, contributing to the overall transmitted signal but at a later time than the unperturbed wave. The second reason is the delay experienced by forward-scattered waves. Indeed, in the regime of geometrical optics (λ/a 1) where scattering occurs essentially forward, the effective wave speed is given by the geometric velocity c ray = c −1 −1 < c 0 .

Effective wave speed
The effective wave speed in medium 1 is an increasing function of frequency, with a shift from c 0 by about −0.7% at 3 mHz for τ = 400 s. The Keller theory is in relative agreement for low frequencies ( f ≤ 2 mHz) but it predicts a constant wave speed at higher frequencies. On the other hand, the measured effective wave speed in medium 2 clearly changes from the homogenized velocity c h at 1 mHz to the geometric velocity c ray at 5 mHz. We note a remarkable agreement at all frequencies between the simulations and the Keller theory for medium 2.

Variance of wave field
The mean of the perturbation is zero, therefore looking at the coherent wave field may not be enough to directly detect multiple scattering because one would only see oscillations mixed within the noise. In the regime of strong perturbations, the coherent part would vanish and only the fluctuating part would remain, solely accessible via second order moments. One can for instance look at the envelope of the signal by studying the variance of the wave field.
As shown in Fig. 9 and in the online movie, it is composed of three parts: a peak corresponding to the variance of the ballistic wave packet, coda waves (late-arriving waves) propagating forward, and coda waves propagating backward. The forward-propagating coda results from waves back-scattered an even number of times in the perturbed medium. The backwardpropagating coda forms a plateau of width 2L and results from single back-scattering. In geophysics, a connection has been made between the functional form of the coda in time domain and the complexity of the scattering medium (e.g., Sato et al. 2012).
We decompose the domain in three regions (before, after and in the random medium) and integrate spatially the variance over each of these three regions at t m = 8500 s, i.e. after the coherent wave packet went through the random medium and just after the plateau of back-scattered signal went out of it: It gives us a measurement of the variance that, respectively, has been back-scattered, transmitted or is still trapped in the slab at this particular time. For medium 2, the back-scattered variance makes up for about 50% of the total variance for τ = 400 s, and 75% for τ = 1 day. The reason for these high amounts is that the spectrum of medium 2 peaks at small scales, therefore more back-scattering takes place than for instance in medium 1 where these values become respectively 20% and 15%.

Dependence on correlation time of the medium
Calculations of an effective medium are easier to carry when the perturbation is frozen because one can work directly in the frequency domain. Therefore, we study here how the effective parameters k i and c eff depend on the correlation time of the medium. Fig. 10 shows the relative errors in the attenuation, e k i , and in the effective wave speed difference, e c , between a given correlation time and the τ → ∞ case at 2, 3 and 4 mHz, for medium 2: e k i being generally positive, the attenuation is underestimated by the frozen-medium approximation. Our understanding is that since the power of the perturbation mostly lies at high wave numbers, the attenuation mostly comes from the small-scale incoherence between the realizations of the wave packets. Therefore, there must be two regimes: one at small values of τ where the attenuation increases with τ, and one at greater values of τ where the attenuation decreases, because persisting scatterers start to create less small-scale incoherence, so less attenuation. The transition between the two regimes corresponds to a resonance, located according to the theory at about τ = 195 s, τ = 180 s and τ = 135 s at 2, 3 and 4 mHz. On the other hand, e c being negative, the approximation overestimates the decrease in effective wave speed, because longer-lived features are better "seen" by the wave packets. The decrease is therefore a monotonic function of τ, with its asymptotic value at τ → ∞ only determined by the value of the ratio of the wave number over the typical size of the scatterer. The error is frequency-dependent and, on average over the three central frequencies, is 29% (respectively −5%) for the attenuation (respectively the effective wave speed difference) at τ = 400 s. As for the variance, we assume τ = 1 day ∞. This is justified as the propagation time in the random medium of length 30 Mm is about 1 h 1 day. We compute therefore e bsc (τ) = (E bsc (τ) − E bsc (1 day))/E bsc (1 day), (30) e out (τ) = (E out (τ) − E out (1 day))/E out (1 day), (31) e tr (τ) = (E tr (τ) − E tr (1 day))/E tr (1 day).
The relative errors at τ = 400 s are then about 30%, 80% and 290% for the back-scattered, trapped and outgoing variance, respectively. Hence it appears that for medium 2, the variance is more sensitive to the correlation time than the coherent wave field, and that the back-scattered coda is less sensitive than the rest of the variance. The relative errors for medium 1 at τ = 400 s are presented for comparison purposes in Table 2. In this case, the frozenmedium approximation overestimates the attenuation. Most of the power is indeed located at large scales, so the attenuation is mostly caused by the large-scale incoherence between the realizations (shifts of the wave packets), which triggers a broadening and damping of the coherent wave packet. Therefore, the impact of scattering is larger if the scatterers persist while the wave packets propagate through them than if the scatterers evolve in time. On the other hand, the frozen-medium approximation still overestimates, albeit by a larger amount, the decrease in effective wave speed. The error for both quantities does not depend much on frequency, and is about −25% for the attenuation and −19% for the effective wave speed. For the back-scattered coda, the error increases to 46%. Fig. 10. Relative error on the attenuation (top) and the effective wave speed difference (bottom) at 2, 3 and 4 mHz (medium 2). The error is between the quantities at τ and at τ → ∞. The dashed lines are the predictions from the time-dependent Keller theory. Fig. 11. Relative error in the variance integrated in space before, in and after the random medium, at t = 8500 s. The error is between the quantities at τ and at τ = 1 day.

Accuracy of the theories
All theories predict a decrease in the effective wave speed. The effective wave speed and the attenuation of the coherent wave field are best described by the Keller approximation. The Born second-order solution, although consistent with the Keller solution for small perturbations, performs poorly for larger am- Table 2. Relative error at τ = 400 s for the measured quantities for media 1 and 2. e k i and e c are averaged over the three central frequencies.
Medium 1 Medium 2 Notes. For the coherent wave field, the errors are computed using both the temporal and the frequency codes. They are averaged over the three central frequencies 2, 3 and 4 mHz. On the other hand, since we study the variance in time domain, the errors for this quantity are computed using only the temporal code.
plitudes, therefore it may not be suited for the study of acoustic wave scattering by solar granulation unless it is on small distances (< 30 Mm). The homogenization technique and the geometrical optics do not model the attenuation of the coherent wave field. However they correctly represent the decrease in wave speed for low and high frequencies, respectively.

Validity of the frozen-medium approximation
It is more convenient to study acoustic wave propagation in the frequency domain, but this is easily doable only when the coefficients of the wave equation do not depend on time, i.e. when one can use a snapshot of the random medium. As summarized in Table 2, we find that for medium 2, the attenuation is underestimated by the frozen-medium approximation by 29% at the frequencies of interest for the Sun. As for the effective wave speed difference, which is an important quantity since it is directly related to the helioseismic travel times, it is overestimated by 5%. The greater error for k i seemingly arises from the presence of a resonance of the function k i (τ) at a correlation time close to that of granulation, while the effective wave speed does not exhibit such a feature. We note that the relative error in c eff − c 0 is similar to that of k i in medium 1, when the power of the perturbation is distributed at low scales. The frozen-medium approximation underestimates the variance of the amplitude of back-scattered coda waves by about 30%.

Detectability of coda waves
The numerical simulations show the emergence of coda waves, which are an interesting effect of multiple scattering present both in single realizations of the wave field and in its variance, but not in the coherent wave field. Coda waves are seen trailing the ballistic wave packet, and also as late arrival back-scattered waves (in one dimension). In helioseismology, acoustic waves are measured via the two-point cross-covariance function of the solar oscillations. Therefore, in order to identify coda waves in the Sun, one needs to study the statistical variance of this crosscovariance function.

Appendix A: Keller approximation: Time-independent random medium
Starting from a time-independent random medium c(x), we can take the Fourier transport of the wave equation: The autocorrelation written in Eq.
(2) can be simplified to For clarity, we drop the argument ω in the expression ofφ. Keller (1964) considers an unbounded spatially random medium and assumes statistical homogeneity, isotropy and stationarity. The calculation could be generalized to the case of a localized perturbation, however we follow the original derivation. It does accurately model our problem since the amplitude attenuation and the effective wave speed shift arise because of the perturbed region. Therefore only the boundary effects are not taken into account. Keller made the first part of his derivation in time-domain, using the fact that the Green's function for the 3D wave equation is essentially a delta function, which simplifies the calculation. In 1D however, the Green's function is related to the Heaviside step function. We shall first derive the Keller solution in frequency domain for a frozen medium, then generalize in Appendix B to the solution in time domain.
The wave equation given by Eq. (A.1) can be written as The unperturbed equation, assuming a constant background sound speed, is The corresponding Green's function G 0 , solution of L 0G0 (x, x ) = δ(x − x ) where δ is the Dirac delta function, is where k 0 = ω/c 0 . Keller has shown that one can find a new wave equation for the coherent wave field under the form We assume that the coherent wave field also satisfies a wave equation with a complex wave number k so that .11) In this case, Therefore, On the other hand, We can define the complex wave number by Since the autocorrelation function of the perturbation depends here only on the difference x − x, I(x) = I. In the smallperturbation approximation, one can also replace k by k 0 in the right-hand term, to get finally We note that it is possible to keep k in the right-hand side, one then has to solve a biquadratic complex equation. Here we only use the approximation. In this paper, we used in one case an exponential correlation function (A.19) where ζ = x − x. In this case ∂ ζ f 1 (ζ) = −sign(ζ) f 1 (ζ)/a and ∂ 2 ζ f 1 (ζ) = f 1 (ζ)/a 2 − 2 a δ(ζ), so that Thus k 2 = k 2 0 + 2 k 2 0 3 − 4(k 0 a) 2 1 + 4(k 0 a) 2 +2i 2 k 2 0 k 0 a + k 0 a 1 + 4(k 0 a) 2 . (A.21) This formula gives the damping Im(k) = k i of the coherent wave φ and the effective wave speed ω/Re(k) = c eff of the medium. For medium 2, we evaluate the integral numerically.

Appendix B: Keller approximation: Time-dependent random medium
Here, we extend the previous analysis to a time-dependent random medium c(x, t). We rewrite the problem as follows: The associated Green's function, solution of where Θ is the Heaviside step function. With these new operators, writing the wave field as The calculations are similar to those for the time-independent random medium. Replacing again k by k 0 in the O( 2 ) terms, one gets for medium 1 We have demonstrated here the possibility to develop a timedependent theory given the knowledge of the power spectrum (or autocorrelation function) of the perturbation. We note that here too, the solution for medium 2 presented in the corpus is evaluated numerically.

Appendix C: Second-order Born approximation
Another theory is the second-order Born approximation, which we derive here for a time-independent random medium c(x). It is similar to the Keller theory, but one does not look for an effective wave equation satisfied by the mean wave field. Instead, one writes the mean wave field as a series up to a certain order, each term being proportional to a power of . Using the same notations for the operators as in Appendix A, denotingφ 0 the unperturbed wave field andφ 1 the correction such thatφ =φ 0 +φ 1 , the 1st-order Born approximation reads φ =φ 0 −L −1 0L 1φ0 + O( 2 ). (C.1) Taking the average, one gets φ =φ 0 + O( 2 ). This means that we have to go down to the second order: which, averaged, gives We can compute L 1L −1 0L 1 φ 0 and L 2 φ 0 easily because these are mostly equations A.10 and A.15 replacing φ (x) = e ikx bỹ φ 0 (x) = e ik 0 x . One finally needs to applyL −1 0 which is a convolution by the Green's function. In order to converge, the integration requires a compact support. To model the localization of the perturbation between X and X + L, we introduce the window function w(x) = Θ(x − X) − Θ(x − (X + L)) (C.4) wherex = (x + x )/2, so that δc(x)δc(x ) = 2 e −|ζ|/a w(x). (C.5) The approximate solution in [X, X + L] is φ (x) φ 0 (x) 1 + 2 3 2 ik 0 a − (k 0 a) 2 + (k 0 a) 2 2ik 0 a − 1 (x − X) , (C.6) which, since 1, can be written (omitting a phase term) in the form φ (x) e ik(x−X) where k has the same expression as for the Keller theory (Eq. (A.21)). To this level of approximation, the effective k does not depend on L.

Appendix D: Spatial homogenization
In order to perform the spatial homogenization for a timeindependent random medium c(x), we consider the variable ψ = c 2 φ, (D.1) which is solution of Multiplying the equation by ∂ t ψ and integrating over space, then applying an integration by parts, we find that