Open Access
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/0004-6361/202039201
Published online 20 November 2020

© P.-L. Poulier et al. 2020

Licence Creative CommonsOpen 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 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 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 late-arrival fluctuations due to multiply-scattered (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 one-dimensional background with sound speed c0 = 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 time-dependent random fluctuation:

(1)

This is shown in Fig. 1, where the filled circles symbolize the fluctuation. In Eq. (1), δc has a zero mean so that ⟨c⟩ = c0 where angle brackets denote an expectation value.

thumbnail Fig. 1.

Schematics of the problem.

The sound-speed perturbation is specified through the autocorrelation

(2)

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

(3)

In time, we choose an exponential profile

(4)

where τ is the e-folding lifetime. For granulation, we have τ ≈ 400 s (e.g., Title et al. 1989). The temporal power spectrum is Lorentzian,

(5)

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:

(6)

For a granulation-like medium, it is reasonable to choose a = 1 Mm. In Fourier space,

(7)

The second choice (hereafter Medium 2) is a spatial power spectrum of the form (e.g., Baran 2013)

(8)

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−4R where R = 696 Mm, such that the spatial power spectrum peaks at k R = 1500. In real space, for α = 1, we have

(9)

The two power spectra and their corresponding autocorrelation functions are shown in Fig. 2.

thumbnail 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:

(10)

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 (Lynden-Bell & Ostriker 1967)

(11)

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

(12)

In this paper, we assume that the density is constant and consider the following 1D acoustic wave equation

(13)

We implement two numerical codes. The first code is a time-domain code to study the propagation of a wave packet through a time-dependent random medium, based on Eq. (13). As initial condition, we inject at location x0 <  X a wave packet of central frequency ω0 and frequency width σ:

(14)

where

(15)

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 back-scattered 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 = xmax.

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

(16)

with Dirichlet boundary condition at x = 0,

(17)

and the Sommerfeld outgoing radiation condition

(18)

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 c0, depending on frequency, called the effective wave speed.

An approximate representation of the coherent wave field inside the perturbed medium is therefore

(19)

where the effective wave number is

(20)

and the spatial attenuation is

(21)

The effective wave speed is defined by

(22)

Similarly, we define k0(ω), the wave number for an unperturbed wave field, such that

(23)

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.

3. Numerical methods

3.1. Numerical scheme to solve for ϕ(x,t)

In order to solve numerically Eq. (13), we use an explicit finite-difference 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 x0 = 100 Mm, while xmax = 200 Mm and tmax = 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 c0Δtx = 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.

thumbnail 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 second-order discretization scheme with a spatial resolution δx = 4 km. A tridiagonal system is inverted with the tridiagonal matrix algorithm (Thomas algorithm). c0, xmax, X and L are the same as for the time-domain 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 x1 and x2 (x2 >  x1) in an attenuating medium, a plane wave is damped by a factor eki(x2 − x1). The spatial attenuation ki 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 x1 = 126 Mm and x2 = 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 ∈ [x1, x2]. The power of the signal has been attenuated during the propagation from x1 to x2. 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.

thumbnail Fig. 4.

Measuring the attenuation with the temporal code. Top: Coherent wave packet at x1 = 126 Mm (blue) and x2 = 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 [x1, x2]. The vertical dashed lines delimit the location of the perturbation in time (for a wave packet propagating at c0) 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

(24)

where A, xs, and ceff are the free parameters, with xs being a phase shift. A similar fit is done on ϕ0 to take numerical dispersion into account. This is illustrated in Fig. 5.

thumbnail 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 x2. 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 k0a and ϵ, different theories can be used to compute the effective parameters ceff and ki. 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 ceff = c0(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 ceff = cray = c0(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 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.

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.

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 Sect. 5.1.

If k0a ≪ 1 or k0a ≫ 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 second-order theories are consistent and predict that ki is proportional to while the effective wave speed difference is essentially independent of k0 provided that k0a ⪆ 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 k0a ≪ 1 and to for k0a ≫ 1 while the behavior of the effective wave speed is in qualitative agreement. In particular, the effective wave speed is always less than c0. 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 k0a, 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 ceff − c0 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 104 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.

thumbnail 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 time-dependent Keller theory.

thumbnail 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 ki/k0 is a linear function of frequency, meaning that ki 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 ki/k0 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 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 c0. 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 cray = ⟨c−1−1 <  c0.

thumbnail 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 c0 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 ch at 1 mHz to the geometric velocity cray 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 (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 backward-propagating 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).

thumbnail 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 tm = 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:

(25)

(26)

(27)

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%.

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 ki and ceff depend on the correlation time of the medium.

Figure 10 shows the relative errors in the attenuation, eki, and in the effective wave speed difference, ec, between a given correlation time and the τ → ∞ case at 2, 3 and 4 mHz, for medium 2:

(28)

(29)

thumbnail 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.

eki 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, ec 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.

Table 2.

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

(30)

(31)

(32)

As shown in Fig. 11, 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.

thumbnail 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 frozen-medium 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%.

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 second-order 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 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 ki seemingly arises from the presence of a resonance of the function ki(τ) 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 ceff − c0 is similar to that of ki 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%.

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 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 cross-covariance 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

  1. Aki, K., & Richards, P. G. 2002, Quantitative Seismology, 2nd edn. (University Science Books) [Google Scholar]
  2. Aki, K., & Wu, R. S. 1988, Scattering and Attenuation of Seismic Waves, Part I (Springer) [CrossRef] [Google Scholar]
  3. Ball, W. H., Beeck, B., Cameron, R. H., & Gizon, L. 2016, A&A, 592, A159 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Baran, O. A. 2013, Adv. Astron. Space Phys., 3, 89 [Google Scholar]
  5. Bhattacharya, J., Hanasoge, S., & Antia, H. M. 2015, ApJ, 806, 246 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bourret, R. C. 1963, Appl. Sci. Res. Sect. A, 12, 223 [CrossRef] [Google Scholar]
  7. Brown, T. M. 1984, Science, 226, 687 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  8. Delache, P., & Fossat, E. 1988, in Seismology of the Sun and Sun-Like Stars, ed. E. J. Rolfe, ESA SP, 286, 671 [Google Scholar]
  9. Duvall, T. L., Jr, Kosovichev, A. G., & Murawski, K. 1998, ApJ, 505, L55 [NASA ADS] [CrossRef] [Google Scholar]
  10. Gizon, L., & Birch, A. C. 2004, ApJ, 614, 472 [NASA ADS] [CrossRef] [Google Scholar]
  11. Hanasoge, S. M., Gizon, L., & Bal, G. 2013, ApJ, 773, 101 [NASA ADS] [CrossRef] [Google Scholar]
  12. Houdek, G., Trampedach, R., Aarslev, M. J., & Christensen-Dalsgaard, J. 2017, MNRAS, 464, L124 [NASA ADS] [CrossRef] [Google Scholar]
  13. Howe, M. S. 1971, J. Fluid Mech., 45, 785 [NASA ADS] [CrossRef] [Google Scholar]
  14. Keller, J. B. 1964, Proceedings of Symposia in Applied Mathematics (Providence, R.I.: American Mathematical Society) [Google Scholar]
  15. Legendre, G. 2003, Ph.D. Thesis, Universite paris VI [Google Scholar]
  16. Lynden-Bell, D., & Ostriker, J. P. 1967, MNRAS, 136, 293 [Google Scholar]
  17. Müller, G., Roth, M., & Korn, M. 1992, Geophys. J. Int., 110, 29 [CrossRef] [Google Scholar]
  18. Murawski, K., & Roberts, B. 1993a, A&A, 272, 595 [NASA ADS] [Google Scholar]
  19. Murawski, K., & Roberts, B. 1993b, A&A, 272, 601 [NASA ADS] [Google Scholar]
  20. 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]
  21. Rosenthal, C. S., Christensen-Dalsgaard, J., Nordlund, Å., Stein, R. F., & Trampedach, R. 1999, A&A, 351, 689 [NASA ADS] [Google Scholar]
  22. Rytov, S. M., Kravtsov, Y. A., & Tatarskii, V. I. 1989a, Priniciples of Statistical Radiophysics, 3. Elements of Random Fields (Springer) [Google Scholar]
  23. Rytov, S. M., Kravtsov, Y. A., & Tatarskii, V. I. 1989b, Principles of Statistical Radiophysics, 4. Wave Propagation Through Random Media (Springer) [Google Scholar]
  24. Sato, H., Fehler, M. C., & Maeda, T. 2012, Seismic Wave Propagation and Scattering in the Heterogeneous Earth, 2nd edn. (Springer) [CrossRef] [Google Scholar]
  25. Schou, J., & Birch, A. C. 2020, A&A, 638, A51 [CrossRef] [EDP Sciences] [Google Scholar]
  26. Title, A. M., Tarbell, T. D., Topka, K. P., et al. 1989, ApJ, 336, 475 [NASA ADS] [CrossRef] [Google Scholar]
  27. van der Baan, M. 2001, Geophys. J. Int., 145, 631 [CrossRef] [Google Scholar]

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:

(A.1)

The autocorrelation written in Eq. (2) can be simplified to

(A.2)

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

(A.3)

where

(A.4)

(A.5)

(A.6)

The unperturbed equation, assuming a constant background sound speed, is

(A.7)

The corresponding Green’s function G0, solution of where δ is the Dirac delta function, is

(A.8)

where k0 = ω/c0. Keller has shown that one can find a new wave equation for the coherent wave field under the form

(A.9)

with

(A.10)

We assume that the coherent wave field also satisfies a wave equation with a complex wave number k so that

(A.11)

In this case,

(A.12)

Therefore,

(A.13)

where

(A.14)

On the other hand,

(A.15)

Using Eqs. (A.13) and (A.15) in Eq. (A.9), the perturbed wave equation for the coherent wave field is

(A.16)

We can define the complex wave number by

(A.17)

Since the autocorrelation function of the perturbation depends here only on the difference x′−x, I(x) = I. In the small-perturbation approximation, one can also replace k by k0 in the right-hand term, to get finally

(A.18)

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 ∂ζf1(ζ) = − sign(ζ)f1(ζ)/a and , so that

(A.20)

Thus

(A.21)

This formula gives the damping Im(k) = ki of the coherent wave and the effective wave speed ω/Re(k) = ceff 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:

(B.1)

where

(B.2)

(B.3)

(B.4)

The associated Green’s function, solution of L0G0(t, t′,x, x′) = δ(t − t′)δ(x − x′), is

(B.5)

where Θ is the Heaviside step function. With these new operators, writing the wave field as

(B.6)

it follows that

(B.7)

and

(B.8)

The calculations are similar to those for the time-independent random medium. Replacing again k by k0 in the O(ϵ2) terms, one gets for medium 1

(B.9)

where

(B.10)

(B.11)

(B.12)

(B.13)

(B.14)

We have demonstrated here the possibility to develop a time-dependent 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 the unperturbed wave field and the correction such that , the 1st-order Born approximation reads

(C.1)

Taking the average, one gets . This means that we have to go down to the second order:

(C.2)

which, averaged, gives

(C.3)

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

(C.4)

where , so that

(C.5)

The approximate solution in [X, X + L] is

(C.6)

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 time-independent random medium c(x), we consider the variable

(D.1)

which is solution of

(D.2)

Multiplying the equation by ∂tψ and integrating over space, then applying an integration by parts, we find that

(D.3)

where

(D.4)

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 y0, a slow-varying spatial scale, and y1 = y0/η, a fast-varying spatial scale, where η = k0(ω0)a ≪ 1 (e.g., Hanasoge et al. 2013). Then

(D.5)

and

(D.6)

(D.7)

We also expand the solution

(D.8)

where ψi = ψi(y0, y1, t) = ψi(y0, y1 + a, t). We can now proceed to solving the equation order by order. Order η−2 gives

(D.9)

Multiplying by ψ0, integrating over y1 and using the argument of periodicity, one gets

(D.10)

meaning that ψ0 does not depend on y1. Order η−1 then gives

(D.11)

meaning that ψ1 does not depend on y1 either. Finally, at order η0,

(D.12)

Integrating over the fast-varying coordinate y1, invoking periodicity, one finds the following homogenized equation for ψ:

(D.13)

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 ch of the medium is therefore equal to ⟨c−2−1/2. Knowing that c = c0 + δc, and . Hence:

(D.14)

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 infinite-frequency approximation. In practice the applicability conditions are (Rytov et al. 1989b):

(E.1)

(E.2)

(E.3)

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:

(E.4)

Assuming ergodicity of the random medium, the spatial average identifies with the statistical average and

(E.5)

Appendix F: Comparing theories with numerical simulations in the limit ϵ→ 0

Figure F.1 summarizes the accuracy of the (frozen) Keller theory, the Born second-order approximation, the spatial homogenization and the ray theory in the small-perturbation regime (ϵ = 0.01). For each simulation, ten sets of 104 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 time-domain simulation differs from the attenuation from the frequency-domain one, likely because of numerical diffusion. As k0a 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).

thumbnail 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 dashed-dotted 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

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.

Table 2.

Relative error at τ = 400 s for the measured quantities for media 1 and 2.

All Figures

thumbnail Fig. 1.

Schematics of the problem.

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

Measuring the attenuation with the temporal code. Top: Coherent wave packet at x1 = 126 Mm (blue) and x2 = 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 [x1, x2]. The vertical dashed lines delimit the location of the perturbation in time (for a wave packet propagating at c0) 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
thumbnail 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 x2. 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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.

In the text
thumbnail 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
thumbnail 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 dashed-dotted 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 (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.