Issue 
A&A
Volume 643, November 2020



Article Number  A168  
Number of page(s)  13  
Section  The Sun and the Heliosphere  
DOI  https://doi.org/10.1051/00046361/202039201  
Published online  20 November 2020 
Acoustic wave propagation through solar granulation: Validity of effectivemedium theories, coda waves^{⋆}
^{1}
Max Planck Institute for Solar System Research, JustusvonLiebigWeg 3, 37077 Göttingen, Germany
email: poulier@mps.mpg.de
^{2}
GeorgAugustUniversität, FriedrichHundPlatz 1, 37077 Göttingen, Germany
Received:
14
August
2020
Accepted:
29
September
2020
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 (∼5 min) are similar, the medium in which the waves propagate cannot a priori be assumed to be time independent.
Aims. We compare various effectivemedium 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 secondorder 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 (finiteamplitude soundspeed 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 ballistic 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 timedependent medium. Multiple scattering cannot be ignored when modeling acoustic wave propagation through solar granulation.
Key words: Sun: oscillations / Sun: granulation / waves / scattering / Sun: helioseismology
Movies associated to Figs. 3 and 9 are available at https://www.aanda.org
© P.L. Poulier et al. 2020
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Open Access funding provided by Max Planck Society.
1. Introduction
Solar seismic waves interact with smallscale convective motions near the solar surface via a wave scattering process, affecting their properties (e.g., propagation speed, frequency, amplitude, and phase). As the efolding 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 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 onedimensional steady medium that contains soundspeed 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 soundspeed 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 forwardscattered energy. This includes latearrival fluctuations due to multiplyscattered (coda) waves.
We state the problem in Sect. 2 and explain the numerical implementation in Sect. 3. Various effective medium theories are reviewed in Sect. 4. We present our results in Sect. 5, and discuss them in Sect. 6.
2. Statement of the problem
2.1. The random medium
We consider a uniform onedimensional background with sound speed c_{0} = 10 km s^{−1}, 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 timedependent 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.
Fig. 1.
Schematics of the problem. 
The soundspeed 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 efolding 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 granulationlike 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 k R_{⊙} = 1500. In real space, for α = 1, we have
The two power spectra and their corresponding autocorrelation functions are shown in Fig. 2.
Fig. 2.
Panel a: power spectrum as a function of the adimensional wave number k R_{⊙}. Panel b: power spectrum as a function of frequency. Panel c: spatial autocorrelation. Panel d: temporal autocorrelation. The vertical dotted lines are drawn at the values of the correlation parameters chosen for medium 1, namely τ = 400 s and a = 1 Mm. 
From the knowledge of the power spectrum P(k, ω), we can compute a realization of the sound speed perturbations as follows:
where 𝒩(ω, k) is a 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 𝒩(k, ω) = 𝒩^{*}(−k, −ω). This way to proceed is based on the assumptions of stationarity and horizontal spatial homogeneity of the medium (e.g., Gizon & Birch 2004).
2.2. The wave equation
The displacement ξ of acoustic waves is given by (LyndenBell & 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 timevarying 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 a wave packet through a timedependent 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 frequencydomain 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.
2.3. 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 finiteamplitude 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.
3. Numerical methods
3.1. 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} = 10 000 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 timedomain 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.
Fig. 3.
Top: wave packet propagation through a realization of a random medium (medium 2, with ϵ = 0.1 and τ = 400 s) located between the vertical dashed lines at different time steps. Middle: average over 10 000 realizations. Bottom: square root of the variance of the wave field. See the movie online. 
3.2. Numerical scheme to solve for
The code uses a secondorder 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 timedomain code. The resolution is done for frequencies between 1 and 5 mHz.
3.3. 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 e^{−ki(x2 − x1)}. 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, we take the temporal Fourier transform of ⟨ϕ(x, t)⟩ where x ∈ [x_{1}, x_{2}]. 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.
Fig. 4.
Measuring the attenuation with the temporal code. Top: Coherent wave packet at x_{1} = 126 Mm (blue) and x_{2} = 144 Mm (red). Bottom: Natural log of the spectrum of the coherent wave packet at different frequencies. 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. 
3.4. Measuring the effective wave speed
We first apply a temporal Fourier transform to ⟨ϕ⟩. Then we fit to (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.
Fig. 5.
Measuring the effective wave speed. Top: Coherent wave packet (red curve) experiencing a travel time shift compared to the unperturbed wave packet (blue curve). The vertical dashed lines represent the arrival times at x_{2}. Bottom: Fit of a decaying cosine to the real part of the temporal Fourier transform in the random region, at ω/2π = 1 mHz. 
4. 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 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, smallangle 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 that 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 largeangle 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.
Theories used in this paper for the effective wave speed and attenuation in a frozen medium as ϵ → 0, and their range of validity.
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 timedependent medium (Appendix B); the Born secondorder 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 Sect. 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 homogenization 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 secondorder theories are consistent and predict that k_{i} is proportional to 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) found that for a 3D frozen medium with an exponential autocorrelation (i.e., similar to medium 1), the attenuation is proportional to for k_{0}a ≪ 1 and to 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 socalled selfaveraging quantities; he derived the effective medium using one realization of a onedimensional 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.
5. Results and comparisons
We compute the properties of the effective medium using the procedure explained in Sects. 3.3 and 3.4, with ϵ = 0.1 and a = 1 Mm. The 1σ error bars shown later on the attenuation and wave speed measurements are obtained from ten sets of 10^{4} realizations.
5.1. 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 (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. 6.
Theoretical solutions compared with the coherent wave field inside the random medium from the numerical simulation, for a frozen medium 1, at 3 mHz. Keller: blue. Born: red. Homogenization: orange. Geometrical optics: dashed green. The black crosses are the numerical solution. 
5.2. Attenuation
Figure 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 min, 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 timedependent Keller theory.
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. 
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.
5.3. Effective wave speed
Figure 8 shows the effective wave speed computed with for τ = 400 s and τ → ∞, as well as the timedependent 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 forwardscattered 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}.
Fig. 8.
Speed of the coherent wave packet vs frequency for media 1 and 2, after propagation through a perturbed medium. The 1D theory adapted from Keller (dashed lines) as well as the prediction from the homogenization theory and the geometrical optics are shown. 1 − σ error bars are shown. 
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.
5.4. 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 (latearriving waves) propagating forward, and coda waves propagating backward. The forwardpropagating coda results from waves backscattered an even number of times in the perturbed medium. The backwardpropagating coda forms a plateau of width 2L and results from single backscattering. 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).
Fig. 9.
Square root of the variance of the wave field as a function of position at a given time t = 8500 s. Top: medium 1. Bottom: medium 2. See the movie online. 
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 backscattered signal went out of it:
It gives us a measurement of the variance that, respectively, has been backscattered, transmitted or is still trapped in the slab at this particular time. For medium 2, the backscattered 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 backscattering takes place than for instance in medium 1 where these values become respectively 20% and 15%.
5.5. 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.
Figure 10 shows the relative errors in the attenuation, e_{ki}, 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:
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 timedependent Keller theory. 
e_{ki} being generally positive, the attenuation is underestimated by the frozenmedium approximation. Our understanding is that since the power of the perturbation mostly lies at high wave numbers, the attenuation mostly comes from the smallscale 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 smallscale 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 longerlived 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 frequencydependent and, on average over the three central frequencies, is 29% (respectively −5%) for the attenuation (respectively the effective wave speed difference) at τ = 400 s.
Relative error at τ = 400 s for the measured quantities for media 1 and 2.
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
As shown in Fig. 11, the relative errors at τ = 400 s are then about 30%, 80% and 290% for the backscattered, 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 backscattered coda is less sensitive than the rest of the variance.
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. 
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 largescale 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 frozenmedium 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 backscattered coda, the error increases to 46%.
6. Discussion
6.1. 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 secondorder solution, although consistent with the Keller solution for small perturbations, performs poorly for larger amplitudes, 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.
6.2. Validity of the frozenmedium 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 frozenmedium 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 frozenmedium approximation underestimates the variance of the amplitude of backscattered coda waves by about 30%.
6.3. 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 backscattered waves (in one dimension). In helioseismology, acoustic waves are measured via the twopoint crosscovariance 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.
Acknowledgments
We thank Aaron C. Birch for useful discussions and comments. PLP is a member of the International Max Planck Research School (IMPRS) for Solar System Science at the University of Göttingen. The computational resources were provided by the German Data Center for SDO through grant 50OL1701 from the German Aerospace Center (DLR).
References
 Aki, K., & Richards, P. G. 2002, Quantitative Seismology, 2nd edn. (University Science Books) [Google Scholar]
 Aki, K., & Wu, R. S. 1988, Scattering and Attenuation of Seismic Waves, Part I (Springer) [CrossRef] [Google Scholar]
 Ball, W. H., Beeck, B., Cameron, R. H., & Gizon, L. 2016, A&A, 592, A159 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baran, O. A. 2013, Adv. Astron. Space Phys., 3, 89 [Google Scholar]
 Bhattacharya, J., Hanasoge, S., & Antia, H. M. 2015, ApJ, 806, 246 [NASA ADS] [CrossRef] [Google Scholar]
 Bourret, R. C. 1963, Appl. Sci. Res. Sect. A, 12, 223 [CrossRef] [Google Scholar]
 Brown, T. M. 1984, Science, 226, 687 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Delache, P., & Fossat, E. 1988, in Seismology of the Sun and SunLike Stars, ed. E. J. Rolfe, ESA SP, 286, 671 [Google Scholar]
 Duvall, T. L., Jr, Kosovichev, A. G., & Murawski, K. 1998, ApJ, 505, L55 [NASA ADS] [CrossRef] [Google Scholar]
 Gizon, L., & Birch, A. C. 2004, ApJ, 614, 472 [NASA ADS] [CrossRef] [Google Scholar]
 Hanasoge, S. M., Gizon, L., & Bal, G. 2013, ApJ, 773, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Houdek, G., Trampedach, R., Aarslev, M. J., & ChristensenDalsgaard, J. 2017, MNRAS, 464, L124 [NASA ADS] [CrossRef] [Google Scholar]
 Howe, M. S. 1971, J. Fluid Mech., 45, 785 [NASA ADS] [CrossRef] [Google Scholar]
 Keller, J. B. 1964, Proceedings of Symposia in Applied Mathematics (Providence, R.I.: American Mathematical Society) [Google Scholar]
 Legendre, G. 2003, Ph.D. Thesis, Universite paris VI [Google Scholar]
 LyndenBell, D., & Ostriker, J. P. 1967, MNRAS, 136, 293 [Google Scholar]
 Müller, G., Roth, M., & Korn, M. 1992, Geophys. J. Int., 110, 29 [CrossRef] [Google Scholar]
 Murawski, K., & Roberts, B. 1993a, A&A, 272, 595 [NASA ADS] [Google Scholar]
 Murawski, K., & Roberts, B. 1993b, A&A, 272, 601 [NASA ADS] [Google Scholar]
 Papanicolaou, G. C., & Varadhan, S. R. S. 1982, Diusion with random coefficients, eds. G. Kallianpur, P. R. Krishnaiah, & J. K. Ghosh (North Holland) [Google Scholar]
 Rosenthal, C. S., ChristensenDalsgaard, J., Nordlund, Å., Stein, R. F., & Trampedach, R. 1999, A&A, 351, 689 [NASA ADS] [Google Scholar]
 Rytov, S. M., Kravtsov, Y. A., & Tatarskii, V. I. 1989a, Priniciples of Statistical Radiophysics, 3. Elements of Random Fields (Springer) [Google Scholar]
 Rytov, S. M., Kravtsov, Y. A., & Tatarskii, V. I. 1989b, Principles of Statistical Radiophysics, 4. Wave Propagation Through Random Media (Springer) [Google Scholar]
 Sato, H., Fehler, M. C., & Maeda, T. 2012, Seismic Wave Propagation and Scattering in the Heterogeneous Earth, 2nd edn. (Springer) [CrossRef] [Google Scholar]
 Schou, J., & Birch, A. C. 2020, A&A, 638, A51 [CrossRef] [EDP Sciences] [Google Scholar]
 Title, A. M., Tarbell, T. D., Topka, K. P., et al. 1989, ApJ, 336, 475 [NASA ADS] [CrossRef] [Google Scholar]
 van der Baan, M. 2001, Geophys. J. Int., 145, 631 [CrossRef] [Google Scholar]
Appendix A: Keller approximation: Timeindependent random medium
Starting from a timeindependent 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 timedomain, 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
where
The unperturbed equation, assuming a constant background sound speed, is
The corresponding Green’s function G_{0}, solution of 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
with
We assume that the coherent wave field also satisfies a wave equation with a complex wave number k so that
In this case,
Therefore,
where
On the other hand,
Using Eqs. (A.13) and (A.15) in Eq. (A.9), the perturbed wave equation for the coherent wave field is
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 righthand term, to get finally
We note that it is possible to keep k in the righthand 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
where ζ = x′−x. In this case ∂_{ζ}f_{1}(ζ) = − sign(ζ)f_{1}(ζ)/a and , so that
Thus
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: Timedependent random medium
Here, we extend the previous analysis to a timedependent random medium c(x, t). We rewrite the problem as follows:
where
The associated Green’s function, solution of L_{0}G_{0}(t, t′,x, x′) = δ(t − t′)δ(x − x′), is
where Θ is the Heaviside step function. With these new operators, writing the wave field as
it follows that
and
The calculations are similar to those for the timeindependent random medium. Replacing again k by k_{0} in the O(ϵ^{2}) terms, one gets for medium 1
where
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: Secondorder Born approximation
Another theory is the secondorder Born approximation, which we derive here for a timeindependent 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 the unperturbed wave field and the correction such that , the 1storder Born approximation reads
Taking the average, one gets . This means that we have to go down to the second order:
which, averaged, gives
We can compute and easily because these are mostly Eqs. (A.10) and (A.15) replacing by . One finally needs to apply 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
where , so that
The approximate solution in [X, X + L] is
which, since ϵ ≪ 1, can be written (omitting a phase term) in the form 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
which is solution of
Multiplying the equation by ∂_{t}ψ and integrating over space, then applying an integration by parts, we find that
where
is an expression for the energy. Since it is invariant, we are certain that the homogenization expansion converges.
The medium is assumed to vary on length scales much shorter than the wave (for solar granulation the length scale a is at least shorter than the wave length of acoustic waves). We moreover assume the periodicity of the medium: c(x) = c(x + a). We separate the spatial variable x into y_{0}, a slowvarying spatial scale, and y_{1} = y_{0}/η, a fastvarying spatial scale, where η = k_{0}(ω_{0})a ≪ 1 (e.g., Hanasoge et al. 2013). Then
and
We also expand the solution
where ψ_{i} = ψ_{i}(y_{0}, y_{1}, t) = ψ_{i}(y_{0}, y_{1} + a, t). We can now proceed to solving the equation order by order. Order η^{−2} gives
Multiplying by ψ_{0}, integrating over y_{1} and using the argument of periodicity, one gets
meaning that ψ_{0} does not depend on y_{1}. Order η^{−1} then gives
meaning that ψ_{1} does not depend on y_{1} either. Finally, at order η^{0},
Integrating over the fastvarying coordinate y_{1}, invoking periodicity, one finds the following homogenized equation for ψ:
where is a spatial average. The homogenization method, used here for a periodic medium, has been generalized to a statistically homogeneous and ergodic random medium, by making the period tend to ∞ (e.g., Papanicolaou et al. 1982). The spatial average identifies then with the statistical average. The homogenized sound speed c_{h} of the medium is therefore equal to ⟨c^{−2}⟩^{−1/2}. Knowing that c = c_{0} + δc, and . Hence:
We note that the spatial homogenization technique does not make an attenuation arise.
Appendix E: Ray approximation
The geometrical optics theory, or ray theory, is an infinitefrequency approximation. In practice the applicability conditions are (Rytov et al. 1989b):
Under these conditions, the wave travel time inside the random medium starting at x = X is computed as an integral of the slowness over the ray path:
Assuming ergodicity of the random medium, the spatial average identifies with the statistical average and
Appendix F: Comparing theories with numerical simulations in the limit ϵ→ 0
Figure F.1 summarizes the accuracy of the (frozen) Keller theory, the Born secondorder approximation, the spatial homogenization and the ray theory in the smallperturbation regime (ϵ = 0.01). For each simulation, ten sets of 10^{4} realizations were generated to get the error bars. For such a small perturbation, we are in the regime of validity of the Born and Keller theories and the results are in agreement with the numerical simulations for the attenuation and the effective wave speed. The attenuation for medium 2 resulting from the timedomain simulation differs from the attenuation from the frequencydomain one, likely because of numerical diffusion. As k_{0}a is of order unity in our setup, we are not a priori in the regime of validity of the homogenization or the geometrical optic theories. However, the geometrical optics is in good agreement with the numerical simulations for medium 1, despite the fact that the condition is not verified in our simulations. Medium 2 exhibits, just like for ϵ = 0.1, a transition from the homogenization regime at small frequencies (< 1 mHz) to the geometrical optics regime at high frequencies (> 5 mHz).
Fig. F.1.
Comparison of theories with simulations for the average wave field (ϵ = 0.01). Top: attenuation. Bottom: effective wave speed. The triangles are for the simulations in frequency domain (τ → ∞), the squares for those in time domain (τ = 1 day). The two dasheddotted blue lines are the Born solutions for media 1 and 2, while the yellow and orange dashed lines are the Keller solutions. 1 − σ error bars are shown. 
Movies
Movie of Fig. 3 (Access here)
Movie of Fig. 9 (Access here)
All Tables
Theories used in this paper for the effective wave speed and attenuation in a frozen medium as ϵ → 0, and their range of validity.
All Figures
Fig. 1.
Schematics of the problem. 

In the text 
Fig. 2.
Panel a: power spectrum as a function of the adimensional wave number k R_{⊙}. Panel b: power spectrum as a function of frequency. Panel c: spatial autocorrelation. Panel d: temporal autocorrelation. The vertical dotted lines are drawn at the values of the correlation parameters chosen for medium 1, namely τ = 400 s and a = 1 Mm. 

In the text 
Fig. 3.
Top: wave packet propagation through a realization of a random medium (medium 2, with ϵ = 0.1 and τ = 400 s) located between the vertical dashed lines at different time steps. Middle: average over 10 000 realizations. Bottom: square root of the variance of the wave field. See the movie online. 

In the text 
Fig. 4.
Measuring the attenuation with the temporal code. Top: Coherent wave packet at x_{1} = 126 Mm (blue) and x_{2} = 144 Mm (red). Bottom: Natural log of the spectrum of the coherent wave packet at different frequencies. 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. 

In the text 
Fig. 5.
Measuring the effective wave speed. Top: Coherent wave packet (red curve) experiencing a travel time shift compared to the unperturbed wave packet (blue curve). The vertical dashed lines represent the arrival times at x_{2}. Bottom: Fit of a decaying cosine to the real part of the temporal Fourier transform in the random region, at ω/2π = 1 mHz. 

In the text 
Fig. 6.
Theoretical solutions compared with the coherent wave field inside the random medium from the numerical simulation, for a frozen medium 1, at 3 mHz. Keller: blue. Born: red. Homogenization: orange. Geometrical optics: dashed green. The black crosses are the numerical solution. 

In the text 
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. 

In the text 
Fig. 8.
Speed of the coherent wave packet vs frequency for media 1 and 2, after propagation through a perturbed medium. The 1D theory adapted from Keller (dashed lines) as well as the prediction from the homogenization theory and the geometrical optics are shown. 1 − σ error bars are shown. 

In the text 
Fig. 9.
Square root of the variance of the wave field as a function of position at a given time t = 8500 s. Top: medium 1. Bottom: medium 2. See the movie online. 

In the text 
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 timedependent Keller theory. 

In the text 
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. 

In the text 
Fig. F.1.
Comparison of theories with simulations for the average wave field (ϵ = 0.01). Top: attenuation. Bottom: effective wave speed. The triangles are for the simulations in frequency domain (τ → ∞), the squares for those in time domain (τ = 1 day). The two dasheddotted blue lines are the Born solutions for media 1 and 2, while the yellow and orange dashed lines are the Keller solutions. 1 − σ error bars are shown. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.