PROPAGATION OF SEISMIC WAVES THROUGH A SPATIO-TEMPORALLY FLUCTUATING MEDIUM: HOMOGENIZATION

The Sun evolves continuously over a broad range of spatio-temporal scales. Small-amplitude waves, which are stochastically excited by the action of vigorous near-surface turbulence, propagate through the interior of the Sun and re-emerge at its surface. The Doppler shifting of spectral absorption lines, whose height of formation is altered by these wave motions, provides a direct measurement of the seismic wave field. Because we measure "noise," i.e., a superposition of a multitude of randomly excited waves, we compute and average auto- and cross-correlations of the wave field (between records at different spatial locations). The wave field of the Sun is described as an ergodic and temporally stationary random process, whose statistics are observed to be Gaussian (e.g., Woodard 1997).

A critical goal of helioseismology is to image the properties of turbulence in the interior, the sub-surface structure of magnetic fields, and internal circulations of plasma. However, these phenomena that we want to image are undergoing constant evolution. The fundamental question of how the properties of a temporally changing medium are imprinted in seismic measurements arises. The manner in which long-wavelength waves couple with small-scale, temporally evolving granulation, and whether this represents a strong scattering regime and can therefore explain helioseismic wave damping (e.g., Duvall et al. 1998), may be addressed using these techniques. Small-scale strong granulation flows are thought to distort normal-mode eigenfunctions (e.g., Brown 1984; Murawski & Roberts 1993; Zhugzhda & Stix 1994; Baldner & Schou 2012) and this may be investigated by deriving effective media and studying their anisotropy.

Homogenization, a mathematical technique based on a two-scale asymptotic analysis, is a treatment of wave propagation in the long-wavelength limit. It gives us a means for deriving an effective medium, replacing a rapidly varying medium by a smoother equivalent that to first order is able to model wave propagation on large wavelengths. An analogy may be drawn to Floquet–Bloch descriptions (Floquet 1883; Bloch 1928), which address the energy eigenfunction of an electron within a periodic potential. The theory of homogenization, easily realized for wave propagation in deterministic periodic media, was extended by Kozlov (1979) and Papanicolaou & Varadhan (1982) to ergodic random media by allowing for the periodicity length scale to become infinite. In other words, Kozlov (1979) and Papanicolaou & Varadhan (1982) realized that an ergodic random medium is essentially a periodic medium, but repeating on infinitely long length scales. Thus the method assists us in interpreting the (mostly) horizontally stationary and ergodic random process by which helioseismic waves are described (e.g., Woodard 1997).

In this article, we choose to study the problem of wave propagation amid spatio-temporally periodically fluctuating array speed fluctuations in the context of a much simplified wave equation. For the theory of homogenization, we refer the reader to, e.g., Bensoussan et al. (1978). Homogenization theory, originally developed for periodic structures, also applies in a variety of random media, under the assumption that their statistics are translationally invariant and ergodic in an appropriate sense (see Jikov et al. 1994). For now it is instructive to deal with the periodic case.

The propagation of linear small-amplitude waves in a spatio-temporally evolving medium with flows and magnetic fields is fairly complicated. A variational treatment of the governing equations and the derivation of the formidable set of full equations may be found in Webb et al. (2005). We merely reproduce them here.

Defining the material derivative of the wave displacement by $\dot{{\boldsymbol{\xi }}}$,

$$\begin{equation} \dot{{\boldsymbol{\xi }}}= \frac{\partial {\boldsymbol{\xi }}}{\partial t} + {\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}{\boldsymbol{\xi }}, \end{equation} \tag{ 1 }$$

the oscillations are described by

$$\begin{eqnarray} \rho \partial _t{\dot{{\boldsymbol{\xi }}}} &=& {-}\rho {\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}\dot{{\boldsymbol{\xi }}}+{\boldsymbol{\nabla }}\cdot [ (\rho c^2 - p){\boldsymbol{\nabla }}\cdot {\boldsymbol{\xi }}\,\,{\boldsymbol{I}}\nonumber\\ &&+ (p + B^2/2)({\boldsymbol{\nabla }}{\boldsymbol{\xi }})^t\nonumber \\ &&-({\boldsymbol{B}}{\boldsymbol{B}}:{\boldsymbol{\nabla }}{\boldsymbol{\xi }}- B^2/2{\boldsymbol{\nabla }}\cdot {\boldsymbol{\xi }}){\boldsymbol{I}}- {\boldsymbol{B}}({\boldsymbol{B}}{\boldsymbol{\nabla }}\cdot {\boldsymbol{\xi }}\nonumber\\ &&- {\boldsymbol{B}}\cdot {\boldsymbol{\nabla }}{\boldsymbol{\xi }}) ] - \rho {\boldsymbol{\xi }}\cdot {\boldsymbol{\nabla }}{\boldsymbol{\nabla }}\phi, \end{eqnarray} \tag{ 2 }$$

where t is time, x is space, ∇ is the covariant spatial derivative, I is the dyadic identity tensor, p(x, t) is the background pressure, ρ(x, t) is the background density, c(x, t) is the sound speed, ${\boldsymbol{B}}({\bf x},t)$ is the background magnetic field, ${\boldsymbol{u}}({\bf x},t)$ is the background flow, gravity ${\bf g} = -{\boldsymbol{\nabla }}\phi$, and S(x, t) is the background entropy. The following tensor notation applies $[{\boldsymbol{\nabla }}{\boldsymbol{\xi }}]_{ij} = \partial _i\xi _j$ and $[({\boldsymbol{\nabla }}{\boldsymbol{\xi }})^t]_{ij} = \partial _j\xi _i$, where ∂_i ≡ ∂/∂xⁱ denotes the covariant spatial derivative.

We consider a medium where the sound speed and density are periodically fluctuating in time on a timescale much shorter than the period of the waves of interest. For simplicity's sake, there will be no background flows or magnetic fields, i.e., Equation (2) with ${\boldsymbol{B}}= {\boldsymbol{0}}$, ${\boldsymbol{u}}= {\boldsymbol{0}}$; no entropy waves, ΔS = 0; and we assume constant background pressure, i.e., ${\boldsymbol{\nabla }}p = {\boldsymbol{0}}$. We consider a periodically varying sound speed c(x, t) such that c(x, t + T) = c(x, t), where T is the periodicity timescale, much smaller than the dominant wave period, and t is time. The differential equation of interest is

$$\begin{equation} \rho \partial _t^2 {\xi }- (\rho \partial _t\ln c)\, \partial _t{\xi }- \rho c^2\nabla ^2 {\xi }= 0, \end{equation} \tag{ 3 }$$

where ξ is the scalar wave field displacement. Temporally fluctuating coefficients can pump energy into the wave system. However, this particular form of the wave equation conserves energy. To demonstrate that this is the case, we divide Equation (3) by ρc², multiply it by ∂_tξ, and integrate over volume to obtain

$$\begin{equation} \int d{\bf x} \, \left(\frac{1}{c^2}\partial _t^2{\xi }- \frac{1}{c^3}\partial _t c\, \partial _t{\xi }- \nabla ^2{\xi }\right) \partial _t{\xi }= 0, \end{equation} \tag{ 4 }$$

which may be manipulated further,

$$\begin{equation} \int d{\bf x} \,\left[ \partial _t\left(\frac{1}{c^2}\frac{\dot{{\xi }}^2}{2}\right) + \partial _t \frac{||{\boldsymbol{\nabla }}{\xi }||^2}{2} \right] = 0, \end{equation} \tag{ 5 }$$

where $\dot{{\xi }}= \partial _t{\xi }$ and we have assumed, for the sake of simplicity, that the boundaries are periodic. This allows us to drop boundary-related integrals when using Gauss's theorem to transform the spatial gradient term into the form of Equation (5). Defining the wave energy as

$$\begin{equation} {\mathcal {E}} = \int d{\bf x} \,\left\lbrace \frac{1}{c^2}\frac{\dot{{\xi }}^2}{2}+ \frac{||{\boldsymbol{\nabla }}{\xi }||^2}{2}\right\rbrace\!{,} \end{equation} \tag{ 6 }$$

Equation (5) indicates that it is an invariant, i.e.,

$$\begin{equation} \partial _t {\mathcal {E}} =0. \end{equation} \tag{ 7 }$$

Given the existence of the energy invariant, we are assured that classical techniques of homogenization are applicable. To facilitate the means of retrieving the homogenized solution, we rewrite the wave equation (Equation (3)) thus:

$$\begin{equation} \partial _t\left(\frac{\partial _t{\xi }}{c}\right) - c\nabla ^2{\xi }= 0. \end{equation} \tag{ 8 }$$

Introducing two temporal scales, t₀ and t₁, the fast and slow scales, respectively, and a small parameter ε = ν₀T ≪ 1, where ν₀ is the dominant wave frequency, we note that the first and second time derivatives are

$$\begin{eqnarray} &&\partial _t = {\partial _{t_0}}+ \frac{1}{\varepsilon }{\partial _{t_1}} \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} &&\partial ^2_t = {\partial _{t_0}}^2 + \frac{2}{\varepsilon }{\partial _{t_1}}{\partial _{t_0}}+ \frac{1}{\varepsilon ^2}{\partial _{t_1}}^2. \end{eqnarray} \tag{ 10 }$$

We expand ξ into the series ξ = ξ₀ + εξ₁ + ε²ξ₂ + O(ε³), where ξ_i = ξ_i(x, t₀, t₁) and ξ_i(x, t₀, t₁ + T) = ξ_i(x, t₀, t₁). At order ε⁻², the equation reads

$$\begin{equation} {\partial _{t_1}}\left(\frac{{\partial _{t_1}}{\xi }_0}{c}\right) = 0. \end{equation} \tag{ 11 }$$

Multiplying by ξ₀, integrating over the fast temporal scale, integrating by parts, and invoking periodicity, we obtain

$$\begin{equation} \int dt_1 {\xi }_0\,{\partial _{t_1}}\left(\frac{{\partial _{t_1}}{\xi }_0}{c}\right) = -\int dt_1\, \frac{1}{c}({\partial _{t_1}}{\xi }_0)^2=0, \end{equation} \tag{ 12 }$$

which is negative definite integral unless ${\partial _{t_1}}{\xi }_0 \equiv 0$, forcing ξ₀ to be only a function of t₀, i.e., ξ₀ = ξ₀(x, t₀). At order ε⁻¹, the expansion provides

$$\begin{equation} {\partial _{t_0}}\left(\frac{{\partial _{t_1}}{\xi }_0}{c}\right) + {\partial _{t_1}}\left(\frac{{\partial _{t_0}}{\xi }_0 + {\partial _{t_1}}{\xi }_1}{c}\right) = 0, \end{equation} \tag{ 13 }$$

or

$$\begin{equation} {\partial _{t_1}}\left(\frac{{\partial _{t_0}}{\xi }_0 + {\partial _{t_1}}{\xi }_1}{c}\right) = 0. \end{equation} \tag{ 14 }$$

Invoking the ansatz ${\partial _{t_1}}{\xi }_1 = h(t_1){\partial _{t_0}}{\xi }_0$, together with the constraint arising from periodicity that

$$\begin{equation} \int _0^T dt_1\, h(t_1) = 0, \end{equation} \tag{ 15 }$$

we have

$$\begin{equation} h(t_1) = \frac{c}{\langle c \rangle } - 1, \end{equation} \tag{ 16 }$$

where

$$\begin{equation} \langle c\rangle = \frac{1}{T}\int _0^T dt_1\, c. \end{equation} \tag{ 17 }$$

Thus the following relationship holds:

$$\begin{equation} \frac{{\partial _{t_0}}{\xi }_0 + {\partial _{t_1}}{\xi }_1}{c} = \frac{1}{\langle c\rangle }{\partial _{t_0}}{\xi }_0, \end{equation} \tag{ 18 }$$

and defining t' = t mod T, the corrector ξ₁ is given by

$$\begin{equation} {\xi }_1 = \left(\frac{1}{\langle c \rangle }\int _0^{t^\prime} dt_1\, c - t^\prime\right) {\partial _{t_0}}{\xi }_0. \end{equation} \tag{ 19 }$$

At order ε⁰, we obtain

$$\begin{equation} {\partial _{t_0}}\left(\frac{{\partial _{t_0}}{\xi }_0 + {\partial _{t_1}}{\xi }_1}{c}\right) + {\partial _{t_1}}\left(\frac{{\partial _{t_0}}{\xi }_1 + {\partial _{t_1}}{\xi }_2}{c}\right) - c\nabla ^2 {\xi }_0 = 0, \end{equation} \tag{ 20 }$$

and integrating over the fast temporal scale, invoking periodicity, and using relationship in Equation (18), we obtain the homogenized differential equation

$$\begin{equation} {\partial _{t_0}}^2{\xi }_0 - \langle c \rangle ^2\nabla ^2{\xi }_0 = 0. \end{equation} \tag{ 21 }$$

As we will see in subsequent sections, the bulk modulus ρc² is typically replaced by a tensorial wave speed that varies as a function of the propagation direction. The effective equation in general acquires greater complexity than the original form. It is also important to note that this equation does not always produce a stable solution, and that oscillating coefficients can destabilize the wave equation (see, e.g., Colombini & Spagnolo 1984). This occurs because oscillating coefficients pump energy into the system and in this particular scenario, we control this by introducing a damping term.

4.1. Sound-speed Perturbation

We study Equation (3) numerically in order to characterize the bounds and effectiveness of homogenization at addressing wave speed perturbations. We solve the equation using a pseudo-spectral solver, computing horizontal derivatives using a fast Fourier transform and evolving it in time through the repeated application of an optimized five-stage second-order Runge–Kutta scheme (Hu et al. 1996). The horizontal boundaries are periodic. We set off a one-way Gaussian wave packet at t = 0, of central wavelength 3.33 Mm, central frequency ω₀/(2π) = 3 mHz, and with a nominal sound speed of c₀ = 10 km s⁻¹. The FWHM of this Gaussian wave packet is 10 Mm. The background medium contains a Gaussian ball shaped sound-speed perturbation, whose FWHM is one wavelength, i.e., 3.33 Mm. The amplitude of the perturbation oscillates in time, varying from 0% to 200% in sound speed. The form of this perturbation is described by

$$\begin{equation} c = c_0\left[1 + A\,\exp \left(- \frac{x^2 + y^2}{\sigma ^2}\right)\,\sin ^2\left(\omega t\right) \right], \end{equation} \tag{ 22 }$$

where c₀ = 10 km s⁻¹, A = 2, and σ = 3.33 Mm. It is seen from Equation (17) that, when ω ≫ ω₀, the homogenized, effective sound-speed squared is given by

$$\begin{equation} \langle c\rangle = c_0\left[1 + \frac{A}{2}\,\exp \left(- \frac{x^2 + y^2}{\sigma ^2}\right)\,\right]. \end{equation} \tag{ 23 }$$

In Figure 1, we show the solution at four different instants for four different cases, where r = ω/ω₀ = [0.5, 1.0, 2.0], and the homogenized, time-stationary sound-speed solution. Figure 2 displays a cut of the wave field along the centerline of the x-axis. The homogenized solution together with the corrector is compared with the full solution, and it is seen that the performance is worst for the r = 0.5, 1 cases, where some form of temporal scattering resonance may be occurring.

Zoom In Zoom Out Reset image size

We now consider a medium where the sound speed and density are temporally stationary but spatially fluctuate periodically on a length scale much shorter than the wavelengths of interest. We consider a periodically varying sound speed c(x, x/ε) and density ρ(x, x/ε), where box ${\mathcal {B}}= [0, L_1)\times [0,L_2)$ describes the unit periodic box (much smaller than the dominant wavelength) that is used to tile the entire domain, x is the "slow" coordinate, and y = x/ε the fast coordinate. The differential equation of interest is

$$\begin{equation} \rho \partial _t^2 {\xi }- {\boldsymbol{\nabla }}\cdot (\rho c^2{\boldsymbol{\nabla }}{\xi }) = 0, \end{equation} \tag{ 24 }$$

where ξ is the wave field and ∇ is the covariant spatial derivative. It may be verified that this differential equation possesses an energy invariant, given by

$$\begin{equation} {\mathcal {E}} = \int d{\bf x} \,\left\lbrace \frac{\rho \dot{{\xi }}^2}{2} + \rho c^2\frac{||{\boldsymbol{\nabla }}{\xi }||^2}{2}\right\rbrace, \end{equation} \tag{ 25 }$$

i.e., $\partial _t{\mathcal {E}}=0$. The existence of an energy invariant guarantees the convergence of a homogenization expansion. The two-scale representation of the spatial derivative is given by

$$\begin{equation} {\boldsymbol{\nabla }}= {\boldsymbol{\nabla }}_{\bf x} + \varepsilon ^{-1}{\boldsymbol{\nabla }}_{\bf y}. \end{equation} \tag{ 26 }$$

Introducing the ansatz ξ = ξ₀(x, x/ε, t) + εξ₁(x, x/ε, t) + ε²ξ₂(x, x/ε, t) + ⋅⋅⋅ Collecting terms of O(ε⁻²), we have

$$\begin{equation} {\boldsymbol{\nabla }}_{\bf y} (\rho c^2 {\boldsymbol{\nabla }}{\xi }_0) = 0, \end{equation} \tag{ 27 }$$

which, invoking the periodicity over L, may be manipulated as follows:

$$\begin{equation} \int _{{\mathcal {B}}} d{\bf y} \, {\xi }_0{\boldsymbol{\nabla }}_{\bf y} \cdot (\rho c^2 {\boldsymbol{\nabla }}{\xi }_0) = -\int _{\mathcal {B}}d{\bf y} \, \rho c^2 ||{\boldsymbol{\nabla }}{\xi }_0||^2 = 0, \end{equation} \tag{ 28 }$$

implying that ${\boldsymbol{\nabla }}_{\bf y} {\xi }_0 = {\bf 0}$ or ξ₀ ≡ ξ₀(x, t). The following definition applies:

$$\begin{equation} \int _{\mathcal {B}}d{\bf y} = \int _0^{L_1} dy_1 \int _0^{L_2} dy_2, \end{equation} \tag{ 29 }$$

where y = (y₁, y₂). This is in line with expectations, since, at leading order, the solution is presumably dominated by effects on the scale of the wavelength. At order ε⁻¹, we obtain

$$\begin{equation} {\boldsymbol{\nabla }}_{\bf y} \,{\cdot}\, (\rho c^2 {\boldsymbol{\nabla }}_{\bf y} {\xi }_1) \,{+}\, {\boldsymbol{\nabla }}_{\bf x} \,{\cdot}\, (\rho c^2 {\boldsymbol{\nabla }}_{\bf y} {\xi }_0) \,{+}\, {\boldsymbol{\nabla }}_{\bf y} \,{\cdot}\, (\rho c^2 {\boldsymbol{\nabla }}_{\bf x} {\xi }_0) = 0. \quad \end{equation} \tag{ 30 }$$

Invoking the result that ξ₀ = ξ₀(x, t), this simplifies to

$$\begin{equation} {\boldsymbol{\nabla }}_{\bf y} \cdot [ \rho c^2 ({\boldsymbol{\nabla }}_{\bf y} {\xi }_1 + {\boldsymbol{\nabla }}_{\bf x} {\xi }_0)] = 0. \end{equation} \tag{ 31 }$$

We seek solutions of the form (e.g., Bensoussan et al. 1978) ${\xi }_1 = {\bf h}({\bf y})\cdot {\boldsymbol{\nabla }}_{\bf x} {\xi }_0$, which implies

$$\begin{equation} {\boldsymbol{\nabla }}_{\bf y} \cdot [\rho c^2({\boldsymbol{\nabla }}_{\bf y} {\bf h} + \,{\bf I})]\cdot {\boldsymbol{\nabla }}_{\bf x} {\xi }_0 = 0, \end{equation} \tag{ 32 }$$

where  I is the identity tensor. This produces the following elliptic equation:

$$\begin{equation} {\boldsymbol{\nabla }}_{\bf y} \cdot [\rho c^2({\boldsymbol{\nabla }}_{\bf y} {\bf h} + \,{\bf I})] = {\bf 0}, \end{equation} \tag{ 33 }$$

whose solution gives us the corrector ${\xi }_1 = {\bf h}\cdot {\boldsymbol{\nabla }}{\xi }_0$. This is the classical cell problem in homogenization. Note that this is an implicit equation that does not, in general, possess a closed form or explicit solution. At order ε⁰, we obtain

$$\begin{equation} \rho \partial _t^2 {\xi }_0 \,{-}\, {\boldsymbol{\nabla }}_{\bf x} \,{\cdot}\, [\rho c^2({\boldsymbol{\nabla }}_{\bf x} {\xi }_0 \,{+}\, {\boldsymbol{\nabla }}_{\bf y} {\xi }_1)] \,{-}\, {\boldsymbol{\nabla }}_{\bf y} \,{\cdot}\, (\rho c^2{\boldsymbol{\nabla }}_{\bf y} {\xi }_2) = 0, \end{equation} \tag{ 34 }$$

or

$$\begin{equation} \rho \partial _t^2 {\xi }_0 \,{-}\, {\boldsymbol{\nabla }}_{\bf x} \,{\cdot}\, [\rho c^2({\boldsymbol{\nabla }}_{\bf y} {\bf h} \,{+} \,{\bf I})\,{\cdot}\, {\boldsymbol{\nabla }}_{\bf x} {\xi }_0] \,{-}\,{\boldsymbol{\nabla }}_{\bf y} \,{\cdot}\, (\rho c^2{\boldsymbol{\nabla }}_{\bf y} {\xi }_2) = 0. \end{equation} \tag{ 35 }$$

Integrating over the fast variable and invoking periodicity,

$$\begin{equation} \rho _*\partial _t^2 {\xi }_0 \,{-}\, {\boldsymbol{\nabla }}_{\bf x} \,{\cdot}\, [{\bf C} _*\,{\cdot}\, {\boldsymbol{\nabla }}_{\bf x} {\xi }_0] = 0, \end{equation} \tag{ 36 }$$

where the following definitions hold

$$\begin{eqnarray} \rho _* &=& \frac{1}{L_1 L_2}\int _{{\mathcal {B}}} d{\bf y} \, \rho ({\bf x},{\bf y}), \end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray} {\bf C} _* &=& \frac{1}{L_1 L_2}\int _{{\mathcal {B}}} d{\bf y} \, \rho c^2({\boldsymbol{\nabla }}_{\bf y} {\bf h} + \,{\bf I}). \end{eqnarray} \tag{ 38 }$$

As a consequence of spatial homogenization, the simple scalar wave speed ρc² has now been transformed to wave speed tensor C whose description is obtained by solving the implicit partial differential equation (Equation (33)). Evidently, the case of spatial homogenization is substantially more complicated than the temporal analog.

Zoom In Zoom Out Reset image size

More complicated yet, and in fact much less understood mathematically, is the practically interesting case of fluctuations in both time and space. It is known for some specific choices of the coefficients ρ and c² that the wave equation may not have bounded solutions (e.g., Colombini & Spagnolo 1984). Again, however, with the existence of an energy invariant, the multi-scale expansion is convergent, in which case a dual spatio-temporal homogenization procedure may be applied. See Bensoussan et al. (1978) for additional details on the homogenization of the wave equation in periodic media.

5.1. Numerical Example

Granulation in the Sun is a process that is spatially horizontally "periodic" and substantially smaller than typical acoustic wavelengths. In this simplistic model of wave propagation through granules, we consider a spatially periodic grid of sound-speed perturbations, shown in Figure 3. The elliptic partial differential equation (Equation (33)) is solved by relaxing it to a steady state. The two components of that vector equation are

$$\begin{eqnarray} {\boldsymbol{\nabla }}\cdot [\rho c^2({\boldsymbol{\nabla }}h_x + {\bf e}_x)] &=& 0,\nonumber \\ {\boldsymbol{\nabla }}\cdot [\rho c^2({\boldsymbol{\nabla }}h_y + {\bf e}_y)] &=& 0, \end{eqnarray} \tag{ 39 }$$

where h_x and h_y are the x and y components of vector h and e_x, e_y are the unit vectors along the x and y axes. The relaxation equation for the h_x equation is given by

$$\begin{equation} \partial _t h_x = {\boldsymbol{\nabla }}\cdot [\rho c^2({\boldsymbol{\nabla }}h_x + {\bf e}_x)], \end{equation} \tag{ 40 }$$

where h_x is now a function of time and space. This diffusion equation eventually relaxes to the steady-state ${\boldsymbol{\nabla }}\cdot [\rho c^2({\boldsymbol{\nabla }}h_x + {\bf e}_x)] = 0$. We temporally evolve Equation (40) until a steady state is achieved. The code is tested against an analytical solution that is known for perturbations that are only functions of one coordinate, e.g., where c² = c²(x) only. In such a case, the wave speed tensor C_* has diagonal components, given by the harmonic and simple means.

Zoom In Zoom Out Reset image size

Once we obtain the vector h, we compute the wave speed tensor using Equation (38). The cell problem studied here has cylindrically symmetry, i.e., each of the sound-speed perturbations is azimuthally symmetric around its axis. This implies that the wave speed tensor will contain identical components in the xx and yy planes, providing us another way to test the code.

Once we obtain the homogenized wave speed tensor, which is essentially a constant sound speed in the xx and yy planes, we can compare the true and homogenized solutions. We set off a source at the center of the computational domain and compare the homogenized and true solutions in Figure 4. The corrector, given by $\xi _1 = {\bf h}\cdot {\boldsymbol{\nabla }}\xi _0$, contains the higher-order scattering terms not fully captured by the homogenized solution. We show the corrector also in Figure 5, resulting in the decoherence of the input wave packet, thereby contributing to observed wave attenuation. The wave fields are practically identical, indicating that the infinite-wavelength limit works very accurately at modeling these "granules."

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the cases we have considered thus far, the media consist of periodic arrays of scatterers. Papanicolaou & Varadhan (1982) showed for random media where the time-independent perturbations have short correlation length scales and are drawn from stationary and ergodic distributions, the expectation value of the wave field is given by the ensemble average of Equation (38). The result rests on the argument that random media are the limiting case of a period medium with an infinite periodicity length scale. See Jikov et al. (1994) for additional details on the theory of homogenization in random media. Although mathematically more difficult to establish than homogenization theory in periodic media, the main conclusions drawn in periodic media typically also hold in random media. Indeed, seeing random media as a limit of periodic media with increasing cell size, we may obtain the homogenized coefficients in random media ρ_* and C_* as the limits given in Equations (37) and (38) as the sizes L₁ and L₂ tend to ∞.

In this section, we show a simple case where we tile the two-dimensional computational domain with a randomly generated square of sound-speed perturbations. The perturbations are drawn from a zero-mean uniform distribution with an amplitude of 8 km s⁻¹, where the nominal sound speed is 10 km s⁻¹. This square is successively increased in size until it is the size of the entire domain. In other words, for a computational domain of 512 × 512 we choose tiles of sizes 32 × 32, 64 × 64, 128 × 128, 256 × 256, and 512 × 512. The full medium is then filtered to remove the top third highest spatial frequencies (up to the spatial Nyquist) in order to prevent aliasing from corrupting the numerical simulation (Orszag's two-thirds rule; Orszag 1971). The dominant power in the spectrum of the fluctuations is on length scales smaller than the peak wavelength of the wave, 3.33 Mm. To leading order, we show that the wave field is the same in all the cases in Figure 6. A more thorough investigation may be performed, where for a given size of the cell, the cell problem (Equation (33)) is solved for a number of realizations and the variance of the homogenized coefficient (Equation (38)) is estimated. It can be shown that the variance of the homogenized coefficient falls as N^−d, where N = L₁ = L₂ is the size of the tile and d describes some rate.

Zoom In Zoom Out Reset image size

Optimal rates of convergence for random coefficients with short-range correlations have been obtained (for a slightly modified problem) by Gloria & Otto (2011). These analyses are difficult and not known for large classes of processes.

The WKBJ approximation represents the zero-wavelength limit of wave propagation, where the scale over which the structure changes is substantially larger than the wavelength. However, the Sun displays structure over a broad range of scales and consequently, the asymptotic infinite-wavelength limit is also very important to understand.

It is believed that small-scale granulation likely plays a critical role in scattering waves and distorting eigenfunctions of normal modes (Brown 1984; Baldner & Schou 2012). From the numerical experiments we have performed here, we find that when spatio-temporal scales are separated, scattering will be very weak and that strong scattering happens only in cases where the scales overlap. This may have important implications for granular scattering of waves, where the spatial scale separation is significant but the temporal scales of its evolution are similar. Not surprisingly, asymptotic methods break down in this strong scattering regime. Modeling these effects is an important step toward interpreting seismic measurements appropriately. Most interestingly, the analysis reveals that the effective medium possesses a tensorial wave speed and can potentially induce anisotropy in wave propagation. We compute the corrector, which represents higher-order scattering that contributes to the overall decoherence of the wave field.

A powerful extension of the periodic case (that we have studied here) is to ergodic random media where the probability density function describing the randomness is translationally invariant. Granulation and supergranulation fall into this regime, both being described by translationally horizontally invariant quasi-random processes.

S.M.H. is funded by NASA grant NNX11AB63G. This work is an effort to understand cross-correlations in helioseismology in the context of DFG SFB 963 "Astrophysical Flow Instabilities and Turbulence" (Project A1). S.M.H. and L.G. thank the Courant Institute, New York University for their hospitality.