A&A 412, 257-265 (2003)
DOI: 10.1051/0004-6361:20031361

Sensitivity kernels for time-distance inversion based on the Rytov approximation

J. M. Jensen^1,2 - F. P. Pijpers¹

1 - Theoretical Astrophysics Center, Institute for Physics and Astronomy, Aarhus University, 8000 Aarhus C, Denmark
2 - Department of Earth Sciences, University of Aarhus, Denmark

Received 22 December 2000 / Accepted 19 August 2003

Abstract
The study of the solar convection zone has been given a new impetus with what promises to be a very powerful new tool in time-distance helioseismology (Duvall et al. 1993). Using data obtained with this technique, it has become possible to do a tomographic reconstruction of structures in the convection zone. Most of the work done on inversion of these data has relied on the ray approximation to deliver sensitivity kernels for the inversion procedures (Kosovichev 1996). Inversions using non ray-theoretical kernels were shown by Jensen et al. (2001). In this paper we go beyond the ray approximation and derive sensitivity kernels based on the Rytov approximation and Green's function theory. We derive an expression for the sensitivity kernel using ray-based Green's functions and show how further approximations can be done to increase the computational efficiency.

Key words: Sun: helioseismology - Sun: general

1 Introduction

For more than a decade helioseismology has been unraveling the secrets of the solar interior by inferring the solar properties from measurements of eigenmode frequencies. From these measurements it has been possible to obtain very detailed knowledge about the global structure of the sun through the use of inversion techniques. Since this method uses the eigenmodes of the entire sun it is difficult to obtain information about very localized structures. A further problem is the inability of the global helioseismic method to resolve any azimuthal variations. Using the newly developed technique of time-distance helioseismology (Duvall et al. 1993) it has become possible to measure travel times of acoustic waves between different points on the solar surface. From these measurements it is possible to obtain information about very localized structures such as sun spots (Kosovichev 1996), convection cells (Kosovichev & Duvall 1997) and meridional flows (Giles et al. 1997). This information is obtained using inversion techniques which are very similar to those used in terrestrial seismology. To use these inversion techniques it is necessary to obtain sensitivity kernels that describe the forward problem. Previous inversion studies of helioseismic time-distance data used sensitivity kernels obtained in the ray approximation (Kosovichev 1996; Kosovichev & Duvall 1997; Giles et al. 1997). This high-frequency approximation has also been widely used in terrestrial seismology with great success even though its validity can be questioned in many cases. The approximation is even more questionable in the solar case as the ratio of wavelength to ray length is larger than in the terrestrial case. The fact that the ray approximation is not strictly applicable for the Sun was pointed out by Bogdan (1997), who showed that a wave packet propagating in a solar-like atmosphere is not confined to the ray path but samples a large volume of the surrounding medium. Spetzler & Snieder (2001) examined the regimes of validity for ray theory and scattering theory. They found that the scattering theory is superior when the inhomogenieties are smaller than the widths of the first Fresnel zone. This is often the case in time-distance helioseismology. The problem with the ray approximation is well known in geophysics and several suggestions about how to calculate sensitivity kernels for sound-speed perturbations without the ray approximation have been made (e.g. Woodward 1992; Snieder & Lomax; 1996; Marquering et al. 1999). Such sensitivity kernels have also been calculated for the helioseismic case (Birch & Kosovichev 2000; Jensen et al. 2000b). Jensen et al. (2000a) applied these to synthetic data to test the inversion procedure while Jensen et al. (2001) used such kernels on solar data.

This paper aims to develop these techniques from geophysics to the helioseismic case and to show the theoretical foundation for the approximate sensitivity kernels which were proposed in Jensen et al. (2000b) and applied to data by Jensen et al. (2001) to study sun-spot structures.

2 Wave equation for acoustic waves in the convection zone

We wish to study the propagation and scattering of sound waves in the solar convection zone using diffraction tomography. An early analysis of a similar problem is found in Slaney et al. (1984) for the case of micro-wave tomography.

In Appendix A we derive a wave equation for acoustic waves traveling in the solar convection zone in the presence of flows. The wave equation can be written in Fourier domain as

$$\left ( \nabla ^2 + \frac{\omega^2 - \omega_{\rm ac}^2}{c^2} \right ) \Psi = 0$$ (1)

where $\Psi$ is the wavefield, $\omega$ is the frequency and $\omega_{\rm ac}$ is the acoustic cut-off frequency. The wave field is given as $\Psi=\sqrt{\rho}c^2 \nabla \cdot \vec{\xi}$ where $\vec{\xi}$ is the displacement, $\rho$ is the density and c is the sound speed. This is a linearized adiabatic wave equation obtained by using the following approximations: the Cowling approximation, nonmoving time-independent background (on the time-scale of the oscillatory motion), no magnetic fields and neglecting internal gravity waves and the effect of gravity on acoustic waves.

3 Scattering

We now wish to write the wave-propagation problem as a scattering problem where sound-speed variations are considered small deviations from a reference model. The sound-speed perturbations are given as $\delta c = c - c_0$, where c₀ is the reference sound-speed model. The wave equation can be rewritten as

$$\left ( \nabla ^2 + k_0^2 \right ) \Psi = - 2 k_0^2 \frac{\delta c}{c_0} \Psi$$ (2)

where k₀ is given as

$$k_0^2 =\frac{\omega^2 - \omega_{\rm ac}^2}{c_0^2} \cdot$$ (3)

In Appendix D we show that k₀ is the wave number in the unperturbed reference media. In writing the wave equation in this way we have neglected the effect of the sound-speed perturbation on the acoustic cut-off frequency. The right hand side of Eq. (2) acts as a source term where the waves are scattered of the perturbations.

4 The Rytov approximation

The wave equation given by Eq. (2) is non linear and in order to solve it one has to apply some approximation. The most widely used approximation is the Born approximation which Birch & Kosovichev (2002) and Gizon & Birch 2002 have used in helioseismology. In the Born approximation one obtains a linear relation between the slowness perturbations of the medium and the scattered wavefield, but it is not straightforward to separate the effects of the perturbation on the phase and amplitude of the scattered waves. This separation occurs naturally in the Rytov approximation which we will use in this paper (for further details on the Rytov approximation see e.g. Ishimaru 1978; Slaney et al. 1984 or Snieder & Lomax 1996). In the helioseismic literature the Rytov approximation has previously been used by Brüggen 2000 to investigate the effect of solar convection of the eigen-mode frequencies.

The Rytov approximation is derived from the wavefield expressed as

$$\Psi = {\rm e}^{ \Phi } ,$$ (4)

where $\Phi$ is a complex phase function. The phase function $\Phi$ can be written as a sum of the phase function for the incoming waves $\Phi_0$ and the phase function of the scattered waves $\Delta \Phi$

$$\Phi =\Phi_0 + \Delta \Phi .$$ (5)

This means that the incoming waves are given by

$$\Psi_0 ={\rm e}^{\Phi_0} ,$$ (6)

and the total wavefield is then given by

$$\Psi =\Psi_0 {\rm e} ^{\Delta \Phi} .$$ (7)

Since all the phase functions are complex functions we shall consider their real and imaginary parts as

 

$$\Phi \gamma +i\theta \;$$ =

$$\Phi_0 \gamma_0 +i \theta_0 \; \;$$ =

$$\Delta \Phi \Delta \gamma +i \Delta \theta \;.$$ = (8)

The total wavefield and the incoming wavefield is then

 

$$\Psi {\rm e}^{\gamma +i\theta } = A ~ {\rm e}^{i\theta}\; ,$$ =

$$\Psi_0 {\rm e}^{\gamma_0 +i\theta_0} = A_0 ~ {\rm e}^{i\theta_0} \; ,$$ = (9)

where $A={\rm e}^{\gamma}$ and $A_0={\rm e}^{\gamma_0}$ are the amplitudes of the wavefields. Inserting Eqs. (8) and (9) into Eq. (7) and equating real and imaginary parts yields

  

$$\Delta \gamma \ln \left ( \frac{A}{A_0} \right ) ,$$ = (10)

$$\Delta \theta \theta - \theta_0 \; .$$ = (11)

From this it is clear that the real part of the complex phase perturbation represents fluctuations in the amplitude of the wave field whereas the imaginary part represents phase fluctuations.

Inserting Eq. (4) into the wave equation yields

$$\left ( \nabla \Phi \right )^2 + \nabla^2 \Phi + k_0^2 = - 2 k_0^2 \frac{\delta c}{c_0}\cdot$$ (12)

From this a wave equation for the phase-function perturbation $\Delta \Phi$ can be obtained as (see Appendix B)

$$\left ( \nabla ^2 + k_0^2 \right ) \Psi_0 \Delta \Phi = - \P... ...\Phi \right ) ^2 - 2 k_0^2 \frac{\delta c}{c_0} \right \}\cdot$$ (13)

To simplify this equation we use the Rytov approximation. Here we assume that $\nabla \Delta \Phi$ is of the same order as the perturbations (or more accurately that the dimensionless quantity $\nabla \Delta \Phi /k_0$ is of the same order as $\delta c/c_0$) and then neglect terms which are second order in the perturbed quantities (see e.g. Chernov 1960). Applying the Rytov approximation we obtain the following equation

 

$$\left ( \nabla ^2 + k_0^2 \right ) \Psi_0 \Delta \Phi = - 2 k_0^2 \frac{\delta c}{c_0} \Psi_0 \; .$$ (14)

Using Green's function theory Eq. (14) can be solved as (cf. Morse & Ingard 1968)

$$\Psi_0 ({\vec r}_1) \Delta \Phi({\vec r}_1) \approx \int_V 2 ... ...ec r}_1) \Psi_0 ({\vec r^\prime}) {\rm d} {\vec r^\prime} \; ,$$ (15)

where $\Delta \Phi({\vec r}_1)$ is the wave field observed at ${\vec r}_1$ and G₀ is the Green's function of the reference medium. The Green's function is a solution to the equation

$$\left ( \nabla ^2 + k_0^2 \right ) G_0 ({\vec r}^\prime \vert {\vec r}_1) = -\delta({\vec r}^\prime-{\vec r}_1) \; ,$$ (16)

where $\delta$ is the delta function. Thus $G_0({\vec r}^\prime \vert {\vec r}_1)$ is the spatial response to a point source at ${\vec r}^\prime$with frequency $\omega$ observed at  ${\vec r}_1$. For a homogeneous three-dimensional medium the Green's function is given by

$$G_0 ({\vec r}_1 \vert {\vec r^\prime})= \frac{{\rm e}^{ik\ver... ...^\prime}) }} {4 \pi \vert{\vec r}_1-{\vec r^\prime\vert}} \; ,$$ (17)

where $\tau({\vec r}_1,{\vec r^\prime})$ is the travel time from ${\vec r}^\prime$ to ${\vec r}_1$.

Inserting the expression for the object function into Eq. (15) and introducing the Rytov wave path, $\mathcal{L}_{\rm R}$, we get

  

$$\Delta \Phi({\vec r}_1) \approx \int_V \frac{\delta c( {\vec r^\p... ...athcal{L}_{\rm R} ({\vec r^\prime} \vert {\vec r}_1) {\rm d}{\vec r^\prime} \;,$$ (18)

$$\mathcal{L}_{\rm R} ({\vec r^\prime} \vert {\vec r}_1) = 2 k_0^2 ... ...r^\prime}\vert{\vec r}_1)\Psi_0 ({\vec r^\prime})}{\Psi_0 ({\vec r}_1)}\; \cdot$$ (19)

If the incoming wavefield is generated by a point source at ${\vec r}_2$ then $\Psi_0$ is the Green's function for the background medium and the Rytov wave path becomes

$$\mathcal{L}_{\rm R} ({\vec r^\prime} \vert {\vec r}_1 , {\vec... ...me\vert{\vec r}_2)} {G_0 ({\vec r}_1\vert{\vec r}_2)} \; \cdot$$ (20)

Figure 1 shows an example of a Rytov wave path. Shown is the imaginary part of the wave path which is the sensitivity kernel for phase perturbations. There are two interesting features of this wave path. First there are the bands of alternating sensitivity; these corresponds to the different Fresnel zones and we can see the wide first Fresnel zone. The first Fresnel zone is the volume from which scattered waves will interfere constructivly with the direct arrival at the point of observation. Second we can see that along the ray path is an area of low sensitivity. Actually the sensitivity along the ray path is zero. This counter-intuitive result has previously been noted by Woodward (1992) and Marquering et al. (1999). It is important to realize that there is zero sensitivity on the ray path only for phase perturbations. The amplitude perturbations have non-zero sensitivity on the ray path. To understand why this is so, one can consider a perturbation placed on the ray path. The perturbation will act as a point source and since this source is located on the ray path the wavefield emitted by the source will arrive at the point of observation in phase with the wavefield which was scattered. Thus the scattered waves will not give rise to any phase perturbation, but may change the amplitude.

	Figure 1: Horizontal and vertical cross section of the imaginary part of a Rytov wave path in a homogeneous medium. The reference velocity is 20 km s-1 and the period is 5 min. Light areas corresponds to positive sensitivity and dark areas to negative sensitivity
Open with DEXTER

5 The correlation function

In time-distance helioseismology the wavefield is not directly measured. The (processed) data are usually in the form of a correlation C between the radial velocity measured at two areas on the solar surface (cf. Duvall 1993; D'Silva 1996a):

$$C({\vec r_1}, {\vec r_2}, \tau) = \int_{0}^{T} V_{\rm r}({\vec r_1}, t+\tau) V^*_{\rm r}({\vec r_2}, t) {\rm d} t \; .$$ (21)

Here ${\vec r}_1$ and ${\vec r}_2$ represent locations on the solar surface; and T is the total time of the observations. Integration over a (small) area around the surface locations is suppressed for clarity of notation. The radial velocity $V_{\rm r}$ is measured as a function of time t and the correlation function is constructed as a function of time delay $\tau$ between the two locations. The asterisk stands for taking the complex conjugate. In this case where the measured velocity is a real time series a complex function is generated using a Hilbert transform. Here we work with the radial velocity whereas what is observed is the line of sight velocity. In local area helioseismology one often works in a flat Sun approximation where these are the same otherwise one has to make a geometrical correction to take this into account. As is the case for the wave field $\Psi$, here too it is convenient to work instead in the Fourier domain. The Fourier transform of C is $\Gamma$:

$$\Gamma = v_r(\omega,{\vec r_1}) v_r^*(\omega,{\vec r_2}) \; .$$ (22)

Here v_r is the Fourier transform of the radial component of the velocity. From the expression for the wave field in Appendix A one finds that $v_r =i B \nabla_r \Psi$, where $B= \rho_0^{-1/2} \omega (\omega^2-\omega_{\rm ac}^2)^{-1}$is constant at the solar surface. Making use of Eq. (7) and expanding yields:

$$\Gamma B^2 \nabla_{r_1} \Psi \nabla_{r_2} \Psi^*$$ =

$$B^2 \nabla_{r_1} \left(\Psi_0 {\rm e}^{\Delta \Phi}\right) \nabla_{r_2} \left( \Psi_0^* {\rm e}^{\Delta \Phi^*}\right)$$ =

$$B^2 {\rm e}^{\Delta\Phi({\bf r_1}) + \Delta\Phi({\bf r_2})^* } \left[\nabla_{r_1}\Psi_0 \nabla_{r_2}\Psi_0^* + \right.$$ =

$$+ \Psi_0^*({\vec r_2})\nabla_{r_1}\Psi_0\nabla_{r_2}\Delta\Phi^* + \Psi_0({\vec r_1})\nabla_{r_2}\Psi_0^*\nabla_{r_1}\Delta\Phi$$

$$+\left.\Psi_0({\vec r_1})\Psi_0^*({\vec r_2})\nabla_{r_1} \Delta\Phi\nabla_{r_2} \Delta\Phi^*\right] \; .$$ (23)

The first of the four terms in square brackets is simply the Fourier transform of the correlation function of the unperturbed medium, $\Gamma_0$. The last of the four terms is of second order in $\Delta \Phi$ and is discarded. Making use of the Fourier transform of the correlation function in the unperturbed medium $\Gamma_0$ and of Eq. (6) yields:

$$\Gamma = \Gamma_0 \left[1 + {\nabla_{r_2}\Delta\Phi^*\over\na... ...rm e}^{\Delta\Phi({\vec r_1}) + \Delta\Phi({\vec r_2})^*} \; .$$ (24)

As long as the Rytov approximation is valid the second and third term between the square brackets are also small compared to unity, and thus changes in the prefactor to the exponential term are negligible. This implies that the Fourier transforms of the correlation functions for the perturbed and unperturbed medium behave much as the wave field itself (cf. Eq. (7)).

When evaluating $\Delta \Phi({\vec r}_1)$ and $\Delta \Phi({\vec r_2})$ one has to remember that the source of the waves is the turbulence generated by convection. Thus the source is a distributed source rather than a localized one. For a distributed noise source, $s({\vec r}^\prime,\omega)$, the background wave field is given as

$$\Psi _0({\vec r},\omega) = \int G_0({\vec r}\vert {\vec r}) s({\vec r}^\prime,\omega) {\rm d} {\vec r}^\prime \; .$$ (25)

It is this expression for $\Psi_0$ one has to use in Eq. (19) when calculating the wave paths. These wave paths for a distributed source should then be used to evaluate $\Delta \Phi$ in the correlation function. Woodard (1997) and Gizon & Birch (2002) give a more in depth discussion of distributed sources and their effect on the crosscorrelation function and sensitivity kernels.

In practice time-distance helioseismology has worked with arrival times obtained in a way described in Appendix C, for which one can work with either the wavefield or the correlation function. Since the wavefield and correlation function behave in much the same way the differences between the arrival times extracted from them should be small. The theory in this paper has been applied to the wavefield but could just as well be applied to the correlation function. For full waveform analysis one should use the correlation function as described here.

6 Sensitivity kernels for wave travel times

Since the imaginary part of the complex phase-function perturbation is sensitive to the phase perturbation it is sensitive to the travel-time delays. In Appendix C we show how to obtain information about the travel-time perturbation from the imaginary part of the phase-function perturbation. For a monochromatic wavefield the travel-time perturbation is given by $\delta \tau_i = \Delta \theta_i / \omega$, where the subscript i denotes the measurement configuration. Together with Eq. (18) this gives the travel-time perturbation as

$$\delta \tau _i= \frac{\Im (\Delta \Phi_i )}{\omega}= - \int_V... ...a c}{c_0} \mathcal{K}_{\omega,i} {\rm d}{\vec r^\prime} \; ,$$ (26)

where $\Im$ denotes the imaginary part and

 

$$\mathcal{K}_{\omega,i} ({\vec r^\prime}) - \frac{\Im \left \{ \mathcal{L}_R ({\vec r}^\prime \vert {\vec r}_1, {\vec r}_2)\right \}} {\omega}$$ = (27)

$$- 2 \frac{k_0^2}{\omega} \Im \left \{ \frac{G_0 ({\vec r^\prime}\... ...vec r}^\prime\vert{\vec r}_2)} {G_0 ({\vec r}_1\vert{\vec r}_2)} \right \}\cdot$$ = (28)

If the source is not monochromatic then we show in Appendix C how to measure both the group- and phase travel-time perturbation. For the group travel-time the perturbation is given by

$$\delta \tau _{g,i}= \frac{\displaystyle \int \vert \Psi_0\ver... ...ystyle \int \vert \Psi_0\vert^2 \omega^2 {\rm d} \omega} \; ,$$ (29)

where $\vert \Psi_0\vert^2$ is the power spectrum of the solar oscillations. Using the above to get the travel-time perturbation gives

$$\delta \tau _{g,i}= \int_V \mathcal{K}_{i,c} \frac{\delta c}{c_0} {\rm d}{\vec r^\prime} \; ,$$ (30)

where

$$\mathcal{K}_{i,c}= \frac{\displaystyle \int \vert\Psi_0\vert^... ...isplaystyle \int \vert\Psi_0\vert^2 \omega ^2{\rm d}\omega} \;$$ (31)

is the frequency-integrated sensitivity kernel for sound-speed perturbations.

In the case of the phase travel time the perturbations are shown in Appendix C to be given as

$$\delta \tau_{p,i}= \frac{\displaystyle\int \vert \Psi_0\vert^... ...aystyle\omega_0 \int \vert \Psi_0\vert^2 {\rm d} \omega} \; ,$$ (32)

where $\omega _0$ is the central frequency of the wave packet. From this expressions sensitivity kernels similar to Eqs. (31) for the group travel-time perturbation can be derived. In the following we will primarily work with the definition for group travel-time perturbations but because of the similarities between Eqs. (29) and (32) very similar results can be obtained for the phase travel-time perturbations.

The equations of this section show how to calculate sensitivity kernels from knowledge of the Green's function for the reference medium and the power spectrum. These equations are a general result which is independent of the way the Green's functions are calculated. In the following sections we will show one way of calculating Green's function and examples of sensitivity kernels obtained using these Green's functions. The kernels given here are a generalized form of what Woodward (1992) called the bandlimited ray path, and Snieder & Lomax (1996) proved a similar result for the travel time of narrow-band signals.

7 Constructing Green's functions

7.1 Ray-based Green's functions

We assume the reference medium is sufficiently smooth for the ray approximation to be valid thus we can solve Eq. (16) using the ray approximation and thereby obtain an expression for $G_0 ({\vec r}\vert{\vec r^\prime})$. In the following we will consider a source at ${\vec r}^\prime$. To find the solution we write G₀ as a function of a slowly varying amplitude and a phase as

$$G_0= A_0 {\rm e}^{i \theta_0}.$$ (33)

In Appendix D we show that the Green's function in this case is given as

$$G_0 ({\vec r^\prime}\vert{\vec r})= A_0 ({\vec r^\prime}\vert{\vec r}){\rm e}^{i\omega \tau_{p0} ({\vec r^\prime,r})} \; ,$$ (34)

where $\tau_{p0}$ is the phase travel time which is given as

$$\tau_{p0} ({\vec r^\prime},{\vec r}) = \int_{\Gamma_i} \frac{1}{v_{p0}(l )} {\rm d} l \; ,$$ (35)

where $v_{p0}=\omega / k$ is the phase speed and $\Gamma_i$ is the ray path from ${\vec r}^\prime$ to ${\vec r}$. The ray path is given by

$$\frac{{\rm d} {\vec r}}{{\rm d}t} = \frac{\partial \omega}{\partial {\vec k}} \; ,$$ (36)

where $\partial \omega / \partial {\vec k}$ is the group speed. The amplitude, A₀ is shown to be given for a cylindrically symmetric problem by

$$A_0 ({\vec r}^\prime\vert{\vec r} )= \frac{1}{4 \pi} \left \v... ...{v_{p0}({\vec r^\prime}) J({\vec r})} \right \vert ^{1/2} \; ,$$ (37)

where $\Theta$ is the declination of the take-off angle and Jis the Jacobian. Eqs. (34)-(37) now fully specify the ray-based Green's function. To find the travel times and Jacobian one has to solve the ray equations. D'Silva (1994 and 1996b) gives the full ray equations for the solar case.

In the regions where the acoustic cut-off frequency is small the group- and phase travel times become identical. The travel-time perturbations corresponds to the group travel-time perturbations given by Eq. (29). In this case the Green's functions are given as

$$G_0 (\vec{r}^\prime \vert \vec{r}) = \frac{1}{4 \pi} \left \v... .../2} {\rm e}^{i \omega \tau_0 (\vec{r} , \vec{r^\prime}) } \; ,$$ (38)

where $\tau_0=\int_{\Gamma_i} c_0^{-1} {\rm d}l$. Here the efficient parameterization of the ray equations given by Gibson et al. (1972) can be used and it is possible to calculate the travel times using a software packet such as "punch'' (Hole & Zelt 1995). The kernels presented in the following all have been calculated in this approximation.

7.2 Global mode approach

Another way of obtaining Green's functions is from the global eigenmodes. If we have the complete set of eigenfunctions ($\xi _n$) and the corresponding eigenvalues ( $\lambda _n$) then the Green's function is given by (Zwillinger 1989)

$$G_0 (\vec{r} \vert \vec{r^\prime}) = \sum_n \frac{\xi _n (\ve... ...})}{\lambda _n \int \xi _n^2 (\vec{r}) {\rm d} \vec{r} } \cdot$$ (39)

Inserting this into Eq. (20) gives an expression for the Rytov wave path in terms of the global- mode eigenfunctions from which one can obtain the sensitivity kernels through Eq. (28) together with Eq. (31). This procedure corresponds to the approach used by Birch & Kosovichev (2000a) to calculate sensitivity kernels.

8 Frequency-integrated sensitivity kernel

To obtain a sensitivity kernel for the travel time we insert Eq. (34) into Eq. (20). This gives an expression for the Rytov wave path as

 

$$\mathcal{L}_{R,i}(\vec{r^\prime}) 2 k_0^2(\vec{r^\prime}) \frac{A_0 (\vec{r^\prime}\vert\vec{r_1} )... ...(\vec{r_1}\vert\vec{r_2})} {\rm e}^{i\omega \Delta \tau_i(\vec{r^\prime)}} \; ,$$ (40)

where

$$\Delta \tau _i (\vec{r^\prime})= \tau_{p0} (\vec{r}_1,\vec{r^... ... (\vec{r^\prime , r}_2) - \tau_{p0} (\vec{r}_1,\vec{r_2}) \; ,$$ (41)

is the time delay between the direct arrival at ${\vec r}_1$ of waves originating at ${\vec r}_2$ and the arrival of the waves which have been scattered at ${\vec r}^\prime$. Inserting Eq. (40) into Eq. (28) gives the single-frequency sensitivity kernel as

$$\mathcal{K}_{\omega,i} = - 2 \frac{k_0^2}{\omega} \frac{A_0 (... ... (\vec{r_1}\vert\vec{r_2})} \sin ( \omega \Delta \tau _i) \; .$$ (42)

A kernel obtained from Eq. (42) with the amplitudes given in Eq. (38) is shown in Fig. 2.

The spectrum of the solar oscillation spans a range of frequencies with a peak around periods of five minutes. To obtain a frequency-integrated sensitivity kernel for the solar case we use a Gaussian power spectrum with a central frequency of 3.15 mHz and a FWHM of 1.5 mHz. Figure 3 shows a frequency-integrated sensitivity kernel obtained, using this power spectrum along with Eqs. (31), and (42). Here we see that the sensitivity kernel has the shape of a hollow banana. The shape of the kernel, with the property that the sensitivity is zero along the ray path, gave rise to what is called the banana-doughnut paradox by Marquering et al. (1999). We see in Fig. 3 that the single-frequency sensitivity kernels interfere destructively in the outer Fresnel zones and constructively in the first Fresnel zone. This illustrates the finding of Snieder & Lomax (1996) that the travel-time perturbations can be written as an averaging integral over the first Fresnel zone. Similar effects were also noted by Woodward (1992).

9 Approximate sensitivity kernels

In the preceding section we calculated frequency-integrated sensitivity kernels using the ray-based Green's functions. In this section we show how examples of approximation which can be made to simplify the amplitude terms or the frequency integration. We then investigate the effects of these approximations on the kernels. The point of this section is twofold. First we wish to shown how to efficiently calculte approximate kernels which are an improvement over the Fresnel volumes or "fat rays'' suggested in terrestrial seismology (e.g. Cervený & Soares 1992; Husen & Kissling 2001) Secondly we wish to show how the sensitivity kernels suggested and used by Jensen et al. (2000a, 2000b, 2001) are related to the frequency integrated kernels.

	Figure 2: Cross sections along and perpendicular to the ray path of single frequency sensitivity kernel obtained from ray based Green's functions for a frequency of 3.2 mHz. The kernels are normalized by c02. The solid line in the upper panel shows the ray path.
Open with DEXTER

As shown in Appendix D the amplitude can be approximated in terms of the travel time as $A_0 \propto \tau_{p0}^{-1}$ which corresponds to the amplitude in the homogeneous case. If we use these approximate amplitudes we find that the sensitivity kernel is given by

$$\mathcal{K}_{\omega,i}^{\rm app} = \frac{K(\vec{r}_1,\vec{r}_... ...prime,r}_2) } \sin \left ( \omega \Delta \tau _i \right ) \; ,$$ (43)

where $K(\vec{r_1},\vec{r_2})$ is a normalization constant determined by constraining the integral of $-\mathcal{K}v_{p0}$(the sensitivity kernels for a phase slowness perturbation) to be equal to the length of the ray path. This is a very efficient way of calculating approximate sensitivity kernels since one only needs to calculate the travel times and the length of the ray path to obtain the sensitivity kernel.

In order to emulate the effect of cancellation outside the first Fresnel zone Jensen et al. (2000b) used a Gaussian envelope to restrict the sensitivity to the first few Fresnel zones for the two dimensional problem. In three dimensions an approximate kernel obtained from this approach is given by:

$$\mathcal{K}_i^{\rm Gauss} = \mathcal{K}_{\omega_0,i} \exp{ \... ...left ( \alpha \omega _0 \Delta \tau_i \right)^2 \right )} \; ,$$ (44)

where $\omega _0$ is the central frequency for the power spectrum and $\alpha$ controls the degree of cancellation outside the first Fresnel zone. In practice we found that an $\alpha$ of 0.13 was a good choice to emulate the effect of the power spectrum we are using.

	Figure 3: Cross sections along and perpendicular to the ray path of sensitivity kernel obtained from ray based Green's functions for a Gaussian frequency profile (see text). The kernels are normalized by c02. The solid line in the upper panel shows the ray path.
Open with DEXTER

If we take Eq. (42) as the central frequency kernel and use a Gaussian to emulate the cancellation outside the first Fresnel zone as in Eq. (44) we obtain for the approximate sensitivity kernel for sound-speed perturbations:

$$\begin{array}{l} \mathcal{K}_{i,u}^{\rm Gauss}= \hspace{6.65... ...eft (\alpha \omega_0 \Delta \tau_i \right)^2 }) . } \end{array}$$

If we replace the amplitudes with the approximate amplitudes of Eq. (43) we get a sensitivity kernel given by

$$\begin{array}{l} \mathcal{K}_{i, u}^{\rm app}= \hspace{6.65cm... ...eft (\alpha \omega_0 \Delta \tau_i \right)^2 }) . } \end{array}$$

Here we have neglected the effects of the acoustic cut-off frequency by writing $\tau_p \approx \tau$ and $v_{p0} \approx c_{0}$The kernel suggested by Jensen et al. (2000b) is very similar to this kernel, the main difference being a factor 1/c₀. If this factor is left out the discrepancy between the approximate kernel and the frequency integrated kernel will increase, thus this kernel represents an improvement over the previously suggested kernels.

	Figure 4: Sensitivity variation with depth across the ray path at a distance of 3 Mm (top) and 20 Mm (bottom) from the source. Shown are $\mathcal{K}_{i,u}$ (solid), $\mathcal{K}_{i,u}^{\rm Gauss}$ (dashed) and $\mathcal{K}_{i,u}^{\rm app}$ (dotted). The kernels are normalized by c0.
Open with DEXTER

Figure 4 shows two cross sections through the different kernels. One cross section is taken halfway between ${\vec r}_1$ and ${\vec r}_2$ at a distance of 20 Mm right where the ray has its turning point. The other is taken much closer to ${\vec r}_1$ at a distance of 3 Mm from it. First of all we can see that the use of a Gaussian to emulate the frequency integration works well as the two kernels, $\mathcal{K}_{i,u}^{\rm Gauss}$ and $\mathcal{K}_{i,u}$, are very similar. One difference between the two kernels is that the approximate kernel is a bit wider than the frequency-integrated kernel. This is due to the factor $\omega ^2$ in the expression for $\mathcal{K}_{\frac{\delta u}{u}}$ (Eq. (31)). This factor increases the influence of the high-frequency single-frequency kernels and since these kernels have smaller first Fresnel zones than the kernel at the central frequency we see this shrinking of the frequency-integrated kernel.

The last kernel, $\mathcal{K}_{i,u}^{\rm app}$, has the same overall shape as the other kernels but the amplitude is more off. It is lower than the other kernels in the cross section close to ${\vec r}_1$ and higher in the other cross section. This is because the amplitudes were approximated using just knowledge of the travel times. Thus although $\mathcal{K}_{i,u}^{\rm app}$ is not in as good agreement with $\mathcal{K}_{i,u}$ as $\mathcal{K}_{i,u}^{\rm Gauss}$ it gives a very efficient way of calculating a kernel which is in reasonable good agreement with the frequency-integrated sensitivity kernel.

Jensen et al. (2001) showed the first inversion results obtained using non-ray sensitivity kernels. The kernels they used were similar to $\mathcal{K}_{i,u}^{\rm app}$, except for a factor of 1/c₀which unfortunately was obmitted in the earlier kernels. Inversions using the new kernels show that the same overall structures can be recovered although some of the details change. Use of the more correct kernels in the inversion procedures should increase the reliability of the results.

10 Comparison with global-mode approach

Figure 5 shows a comparison between our frequency-integrated sensitivity kernels, Born kernels obtained from global modes and ray-theoretical kernels. The kernels shown are one-dimensional kernels for a sound-speed perturbation which is spherical or planar depending on the geometry of the kernels. The Born kernels are described in Birch & Kosovichev (2000a) and are obtained in a spherical geometry whereas the frequency-integrated kernels are calculated for a planar geometry. Both kernels use the Gaussian power spectrum described in Sect. 8. Two different ray-theoretical kernels are shown; one obtained in a spherical geometry which includes the effect of the acoustic cut-off frequency and one obtained in a planar geometry neglecting the acoustic cut-off frequency. All kernels are for an offset of 10 heliospheric degrees.

We here see that the shape of the frequency-integrated kernel and the Born kernel are very similar although there are some differences. The amplitude of the frequency-integrated kernel is larger near the ray turning point than the amplitude of the Born kernel. The right panel of Fig. 5 shows a close up of the kernels near the surface. Here the Born kernel shows some oscillations which the frequency-integrated kernel lacks. We can also see that the frequency-integrated kernel and the Born kernel follow the ray path from the same geometry close to the surface. For the ray kernels the differences near the surface are due to the effect of the acoustic cut-off frequency while the differences near the turning point are due to geometrical effects. The differences between the frequency-integrated kernel and the Born kernel can be due to geometrical effects, the neglecting of the acoustic-cut-off frequency, different filtering and the use of different approximations (Rytov and Born approximations). The exact causes of the discrepancies are not yet fully understood.

Birch & Kosovichev (2001) showed a comparison between one-dimensional sensitivity kernels obtained from a global-mode approach in the Born approximation (Birch & Kosovichev 2000a) and ray theoretical sensitivity kernels. Their Fig. 1 shows the sensitivity kernels for an offset of 22.5 degrees. By comparing this figure with our Fig. 5 we can see that since the ratio of wavelength to ray length is smaller for the larger offset the kernel is closer to the ray kernel although there still are larger discrepancies near the upper and lower turning point.

Figure 5: Comparison of 1D sensitivity kernels for an offset of 10 heliospheric degrees. Shown are a frequency-integrated sensitivity kernel (solid) and a ray-approximation kernel (dashed) both obtained in a planar geometry and without inclusion of the acoustic cut-off frequency These are shown together with a sensitivity kernel obtained from global-mode summation (dash-dotted) and a ray-approximation kernel (dotted), both from a spherical geometry with inclusion of the acoustic cut-off frequency. To the right is shown a close up of the kernels near the surface.

Open with DEXTER

11 Extending the theory

The theory presented in this paper show how to calculate sensitivity kernels for a point to point measurement geometry. In time-distance helioseismology normal practice has been to correlate a central pont with an surrounding annulus or quadrants of such an annulus (e.g. Duvall et al. 1996). D'Silva (2001) showed that the travel time obtained from this procedure is the average of the point to point travel times weighted with the wave field amplitude. The frequency-integrated kernels presented in this paper are an weighted average over frequency. It is straight forward to expand this weighted average to also include the spatial averaging.

The Green's functions used to calculate the sensitivity kernels have been calculated without inclusion of the acoustic cut-off frequency. The equations for the frequency-integrated kernels presented in Sect. 7 are independent of the approximations made in calculation of the Green's functions and the same expressions for the sensitivity kernels can be used to calculate kernels taking the acoustic cut-off frequency into account The same holds true for the expressions for approximate sensitivity kernels except where noted in the text.

In time-distance helioseismology measurements are made for waves traveling in opposite directions. The differences and averages of these measurements are then taken in order to isolate the effects of sound-speed perturbations from flows (see e.g. Kosovichev & Duvall 1997). In this paper only kernels concerning sound-speed perturbations have been calculated and a future goal is to extend the theory to include flows. Another goal for future work is to derive a wave theoretical treatment of wave propagation in the presence of magnetic structures in the solar convection zone. This will be necessary in order to fully investigate sunspots and other magnetic features in the Sun.

12 Concluding remarks

In this paper we show how to compute sensitivity kernels for helioseismic time-distance tomography using the Rytov approximation and the Green's function approach.

The use of non ray-theoretical kernels in the inversion of helioseismic time-distance data increases the accuracy of the inversion results. Preliminary tests using similar two-dimensional kernels (Jensen et al. 2000a) demonstrate this. These kernels along with rapid inversion techniques such as the MultiChannelDeconvolution technique (Jacobsen et al. 1999) provides a powerful tool which would allow reliable fast inversion of helioseismic time-distance data.

Expressions for both the group- and phase travel-time perturbations are given in this paper. Using these expressions the differences between the group- and phase travel-time perturbations could be used to improve the inversions.

It is also shown that the full correlation function can be directly related to the slowness perturbations through a linear inverse equation for which the kernels are derived. This opens the way for full waveform inversion on the time-distance data.

Acknowledgements

The authors wish to thank Aaron Birch for providing the one dimensional Born kernels and ray kernels for the spherical geometry and for fruitfull discussions. JMJ especially wishes to acknowledge Aaron Birchs help regarding sensitivity kernels for flows. This research was supported by the Danish National Research Foundation through the establishment of the Theoretical Astrophysics Center. The authors also wish to thank an anonymous referee for helping improve the original manuscript.

References

 
• Birch, A. C., & Kosovichev, A. G. 2000, Sol. Phys., 192, 193 In the text NASA ADS  
• Birch, A. C., & Kosovichev, A. G. 2001, in SOHO10/GONG2000: Helio- and Astereoseismology at the Dawn of the millennium, ESA SP-464, 187 In the text  
• Bogdan, T. J. 1997, ApJ, 477, 475 In the text NASA ADS  
• Brüggen, M. 2000, MNRAS, 312, 887 In the text NASA ADS  
• Cerverný, V., & Soares, J. E. P. 1992, Geophysics, 57, 902 In the text  
• Chernov, L. A. 1960, Wave Propagation in Random Media (New York: McGraw-Hill) In the text  
• Christensen-Dalsgaard, J., Däppen, W., Ajukov, S. V., et al. 1996, Science, 272, 1286 NASA ADS  
• Christensen-Dalsgaard, J. 1997, Lecture Notes on Stellar Oscillations, 4th ed. (D.f.I. Print, Aarhus) In the text  
• D'Silva, S. 1994, ApJ, 435, 881 NASA ADS  
• D'Silva, S. 1996a, ApJ, 469, 964 In the text NASA ADS  
• D'Silva, S. 1996b, ApJ, 462, 519 NASA ADS  
• D'Silva, S. 2001, ApJ, 547, L81 NASA ADS  
• Duvall, T. L., Jefferies, S. M., Harvey, J. W., & Pomerantz, M. A. 1993, Nature, 362, 430 In the text NASA ADS  
• Gibson, B. S., Odegard, M. E., & Sutton, G. H. 1979, Geophysics, 44, 185  
• Giles, P. M., Duvall, Jr., T. L., & Scherrer, P. H., 1997, Nature, 390, 52 In the text NASA ADS  
• Gizon, L., & Birch, A. C. 2002, ApJ, 571, 966 In the text NASA ADS  
• Gough, D. 1986, in Hydrodynamics and magneto hydrodynamical problems in the Sun and stars, ed. Y. Osaki (Uni. of Tokyo Press), 117  
• Gough, D. O. 1993, Course 7, Linear adiabatic stellar pulsation, in Astrophysical fluid dynamics, Les Houches Session XLVII, ed. J.-P. Zahn, & J. Zinn-Justin (Amsterdam: Elsevier), 399 In the text  
• Hole, A. J., & Zelt, C. 1995, Geophys. J. Int., 121, 427 In the text  
• Husen, S. & Kissling, E. 2001, Physics of the Earth and planetary interiors, 123, 127 In the text  
• Ishimaru, A. 1978, Wave Propagation and Scattering in Random Media. (New York: Academic Press) In the text  
• Jacobsen, B. H., Møller, I., Jensen, J. M., & Effersø, F. 1999, Phys. Chem. Earth, 24, 215 In the text  
• Jensen, J. M., Jacobsen, B. H., & Christensen-Dalsgaard, J. 2000a, In Inverse Methods II, Lecture Notes in Earth Sciences (Berlin-Heidelberg: Springer), 187 In the text  
• Jensen, J. M., Jacobsen, B. H., & Christensen-Dalsgaard, J. 2000b, Sol. Phys., 192, 231 In the text NASA ADS  
• Jensen, J. M., Duvall, T. L., Jacobsen, B. H., & Christensen-Dalsgaard, J. 2001, ApJ, 553, L193 In the text NASA ADS  
• Jensen, F. B., Kuperman, W. A., Porter, M. B., & Schmidt, H. 1994, Computational Ocean Acoustics, AIP series in modern acoustics and signal processing, (New York: AIP Press) In the text  
• Kosovichev, A. G., 1996, ApJ, 461, L55 In the text NASA ADS  
• Kosovichev, A. G., & Duvall, Jr. T. L. 1997, in SCORe'96: Solar Convection and Oscillations and their Relationship, ed. F. P. Pijpers, J. Christensen-Dalsgaard, & C. S. Rosenthal (Dordrecht: Kluwer), 241 In the text  
• Marquering, H., Dahlen, F. A., & Nolet, G. 1999, Geophys. J. Int., 137, 805 In the text NASA ADS  
• Morse, P. M., & Ingard, K. U. 1968, Theoretical Acoustics (New York: McGraw-Hill)  
• Munk, W., Worchester P., & Wunsch C. 1995, Ocean Acoustic Tomography (Cambridge: Cambridge University Press)  
• Priest, E. R. 1982, Solar Magnetohydrodynamics (Dordrecht: D. Reidel Publ.) In the text  
• Slaney, M., Kak, A. C., & Larsen, L. E. 1984, IEEE Trans. Microwave Th. Tech., MTT-32, 860-873 In the text  
• Snieder, R., & Lomax, A. 1996, Geophys. J. Int., 125, 796 In the text  
• Spetzler, J., & Snieder, R. 2001, Geophys. J. Int., 145, 786 In the text NASA ADS  
• Woodward, M. J. 1992, Geophysics, 57, 15 In the text  
• Zhao, J., Kosovichev, A. G., & Duvall, T. L. 2001, ApJ, 557, 384 NASA ADS  
• Zwillinger, D. 1989, Handbook of differential equations (Boston: Academic Press) In the text

  
Online Material

Appendix A: Derivation of wave equation

In this Appendix we derive a wave equation for acoustic waves traveling in the solar convection zone in the presence of subsurface flows and magnetic fields. The derivation for stratified media follows work done by Gough (1993), while Christensen-Dalsgaard (1997) treat flows. Here we present a treatment which include both effects and derive a wave equation for this case.

Governing equations

The governing ideal MHD equations are (see e.g. Priest 1982)

$$\frac{\displaystyle{\rm d}\rho}{\displaystyle{\rm d}t} + \rho \nabla \cdot {\vec V} =0$$ (A.1)

$$\rho \frac{\displaystyle{\rm d}{\vec V}}{\displaystyle{\rm d}t} = - \nabla p + \rho {\vec g}$$ (A.2)

$$\frac{\displaystyle{\rm d}p}{\displaystyle{\rm d}t} = c^2 \frac{\displaystyle{\rm d}\rho}{\displaystyle{\rm d}t}$$ (A.3)

where $\rho$ is the density, p is the pressure, ${\vec V}$ is the fluid velocity, ${\vec g}$ is the gravitational acceleration and c is the adiabatic sound speed. Furthermore ${\rm d}/{\rm d}t$ is the time derivative following a moving fluid element and given as

$$\frac{\displaystyle{\rm d}}{\displaystyle{\rm d}t} = \frac{\d... ...splaystyle\partial t} + \left ( {\vec V} \cdot \nabla \right).$$ (A.4)

We write the involved quantities as the equilibrium values and addition small perturbations associated with the oscillations.

$${\vec V}={\vec U} + {\vec v}$$ (A.5)

$$\rho = \rho_0 + \rho^\prime$$ (A.6)

$$p = p_0 + p^\prime$$ (A.7)

where ${\vec U}$ is the flow velocity and ${\vec v}$ is the fluid velocity associated with the oscillations. In the following we will derive a waveequation in the case where there are no flows present and thus neglect ${\vec U}$. The subscript 0 denotes the equilibrium values and the prime denotes the perturbations. Furthermore we write the time-dependence of the solution as

$$\Psi({\vec r},t)=\Psi({\vec r})\exp(-i\omega t) \; .$$ (A.8)

Stratified media

We now consider the effect of density stratification. Here the perturbed equations are

$${\frac{\displaystyle\partial \rho^\prime}{\displaystyle\partial t} + \nabla \cdot (\rho_0 {\vec v}) = 0 \; ,}$$ (A.9)

$${\rho_0\frac{\displaystyle\partial {\vec v}}{\displaystyle\partial t} = -\nabla p ^\prime +\rho ^\prime {\vec g}}$$ (A.10)

$${\frac{\displaystyle\partial p^\prime}{\displaystyle\partial t} +... ...^\prime}{\displaystyle\partial t} + ({\vec v} \cdot \nabla) \rho_0 \right ).}.$$ (A.11)

From the perturbed equations a wave equation can be obtained for the displacement $\vec{\xi}$ as

$$\rho_0 \frac{\partial^2 \vec{\xi}}{\partial t^2} = \nabla \le... ...cdot {\vec g}) - {\vec g} \nabla \cdot (\rho_0 \vec{\xi}) \; .$$ (A.12)

Here we have also used that the reference medium is in hydrostatic equilibrium and thus we have that $\nabla p_0=\rho_0 {\vec g}$. The last terms on the right hand side, containing ${\vec g}$, give rise to internal gravity waves, the so-called g-modes (e.g. Gough 1993). These terms only affect p-modes propagating deep below the convection zone. Since we are concerned with p-modes in the convection zone we neglect this term. This gives the wave equation as

$$\frac{\partial^2 \vec{\xi}}{\partial t^2} = \rho_0^{-1}\nabla \left ( c^2 \rho_0 \nabla \cdot \vec{\xi} \right ).$$ (A.13)

Following Gough (1993) we now change variable to $\chi=\nabla \cdot \vec{\xi}$ and the wave equation becomes

$$\frac{\partial^2 \chi}{\partial t^2} = \nabla^2 \left ( c^2 \... ... + \nabla \left ( \rho_0^{-1} c^2 \chi \nabla \rho_0 \right ).$$ (A.14)

Making another change of variable to $\Psi$= $\sqrt{\rho_0}c^2\chi$gives the wave equation as

$$\left ( \nabla ^2 - \frac{\omega_{\rm ac}^2}{c^2} - \frac{1}{c^2} \frac{\partial^2}{\partial t^2} \right ) \Psi = 0$$ (A.15)

where

$$\omega_{\rm ac} =\frac{c}{2H}\sqrt{1+2\frac{{\rm d} H}{{\rm d}z}}$$ (A.16)

is the acoustic cut-off frequency and H is the density scale height given as

$$H^{-1}= \frac{{\rm d} \ln \rho}{{\rm d} z} \cdot$$ (A.17)

If we insert the time dependence, $\exp( -i \omega t)$, into the wave equation it reduces to

$$\left ( \nabla ^2 + \frac{\omega^2-\omega_{\rm ac}^2}{c^2} \right ) \Psi =0 \; .$$ (A.18)

Appendix B: Rytov wave equation

We wish to obtain a wave equation for the complex phase perturbation $\Delta \Phi$. One way of doing this was demonstrated by Slaney et al. (1984) and we follow their approach here. Inserting Eq. (4) into the wave equation and using that

$$\nabla ^2 \Psi = \left ( \left ( \nabla \Phi \right )^2 + \nabla^2 \Phi \right ) e^{ \Phi } \; ,$$ (B.1)

yields

$$\left ( \nabla \Phi \right )^2 + \nabla^2 \Phi + k_0^2 = - 2 k_0^2 \frac{\delta c}{c}$$ (B.2)

Inserting Eq. (5) into Eq. (B.2) gives

$$\left ( \nabla \Phi_0 \right )^2 + 2 \nabla \Phi_0 \cdot \nab... ...bla^2 \Delta \Phi + k^2_0 = - 2 k_0^2 \frac{\delta c}{c}\cdot$$ (B.3)

Since $\Phi_0$ is the complex phase function for the incoming waves it satisfies

$$\left ( \nabla \Phi_0 \right )^2 + \nabla^2 \Phi_0 + k_0^2 =0 .$$ (B.4)

Hence Eq. (B.3) reduces to

$$2 \nabla \Phi_0 \cdot \nabla \Delta \Phi + \nabla^2 \Delta \... ...abla \Delta \Phi \right ) ^2 - 2 k_0^2 \frac{\delta c}{c}\cdot$$ (B.5)

To progress further we have to use that

$$\nabla^2 \left ( \Psi_0 \Delta \Phi \right ) = \Delta \Phi \... ..._0 \cdot \nabla \Delta \Phi + \Psi_0 \nabla ^2 \Delta \Phi ~ .$$ (B.6)

As $\Psi_0$ is a solution to the homogeneous Helmholtz equation $\nabla^2 \Psi_0= -k_0^2 \Psi_0$ and from Eq. (6) we find that $\nabla \Psi_0 = \Psi_0 \nabla \Phi_0$. Using this Eq. (B.6) can be rearranged to

$$\Psi_0 \left \{ 2 \nabla \Phi_0 \cdot \nabla \Delta \Phi + \n... ... \} = \left (\nabla^2 + k_0^2 \right ) \Psi_0 \Delta \Phi \; .$$ (B.7)

Inserting Eq. (B.5) into Eq. (B.7) and rearranging gives

$$\left ( \nabla ^2 + k_0^2 \right ) \Psi_0 \Delta \Phi = - \P... ...a \Phi \right ) ^2 - 2 k_0^2 \frac{\delta c}{c} \right \}\cdot$$ (B.8)

This is the wave equation we were looking for.

Appendix C: Determination of the travel time

We wish to determine travel-time differences by fitting the perturbed wavefield to the unperturbed wavefield. The method described here is a generalisation of the method employed by Snieder & Lomax (1996). We assume that the wave field (or correlation function) can be described by a Gabor wavelet which is a cosine modulated with a Gaussian envelope. D'Silva (1996a) and Kosovichev & Duvall (1997) both showed that the correlation function can be written in this way using a ray theoretical approach and a global mode approach respectively. This gives the wave field or correlation function as

$$\Psi_0(t) = A \exp \left [- \frac{\delta \omega ^2}{4} (t - \... ... \right ] \exp \left [ i \omega _0 ( t -\tau_{\rm p}) \right ]$$ (C.1)

where $\tau_{\rm g}$ is the group travel time associated with the peak of the envelope while $\tau_{\rm p}$ is the phase travel time associated with the zero crossings. The amplitude A is a constant, $\omega _0$ is the central frequency of the power spectrum and $\delta \omega$ controls the width of the envelope. Kosovichev & Duvall (1997) showed that the cross correlation function could be expresed in this way. The Fourier transform of this wave field is given as

$$\Psi_0(\omega) = A \exp \left [- \frac{(\omega - \omega_0)^2}... ... \left [ - i \omega _0 ( \tau_{\rm g} -\tau_{\rm p}) \right ].$$ (C.2)

If we look for perturbations to the phase- and group travel times, $\delta \tau_{\rm p}$ and $\delta \tau_{\rm g}$, we find the perturbed wave field as

$$\Psi(\omega) = a \Psi_0 \exp \left ( i \omega \delta \tau_{\r... ...mega _0 ( \delta \tau_{\rm g} - \delta \tau_{\rm p}) \right ].$$ (C.3)

We wish to find the travel time perturbations by fitting this expression for the wave field to the perturbed wave field obtained from the Rytov approximation $\Psi = \Psi_0 {\rm e}^{\Delta \Phi}$. The misfit is given by the function M as

$$\int \left \vert \Psi - a \Psi_0 {\rm e}^{ i \left [ \omega \delt... ...a \tau_{\rm g} - \delta \tau_{\rm p}) \right ] } \right \vert ^2 {\rm d} \omega$$                                    M = (C.4)

$$\int \left \vert \Psi_0 \right \vert^2 \left \vert {\rm e}^{\Delt... ...a \tau_{\rm g} - \delta \tau_{\rm p}) \right ] } \right \vert^2 {\rm d} \omega.$$ = (C.5)

To find the the travel time perturbations we set the derivative of misfit with respect to the perturbation to zero as

$$\frac{\partial M}{\partial \delta \tau_i} =0$$ (C.6)

where $\delta \tau_i$ can be either the phase- or group travel-time perturbation. These derivatives can be written as

$$\frac{\partial M}{\partial \delta \tau_i} = \int \vert \Psi_0... ...rtial \mathcal{A}}{\partial \delta \tau_i} {\rm d} \omega \; ,$$ (C.7)

where

$$\mathcal{A} = \omega \delta \tau_{\rm g} - \Delta \theta - \omega _0 ( \delta \tau_{\rm g} - \delta \tau_{\rm p}) \; .$$ (C.8)

Group travel-time perturbation

Assuming that the travel-time perturbation for the group and phase speed will tend to cancel we get that $\omega _0 ( \delta \tau_{\rm g} - \delta \tau_{\rm p}) \ll \omega \delta \tau_{\rm g}$. Whether this approximation is justified over the entire frequency range of interest remains to be investigated. Using this we find the derivative of M to be

$$\frac{\partial M}{\partial \delta \tau_{\rm g}} = \int \vert ... ...mega \delta \tau_{\rm g} - \Delta \theta ) {\rm d} \omega \; .$$ (C.9)

Expanding the sine term to first order we obtain

$$\int \vert \Psi_0\vert^2 \omega {\rm e}^{\Delta \gamma} (\omega \delta \tau_{\rm g} - \Delta \theta) {\rm d} \omega = 0\; .$$ (C.10)

From this we find that the group speed perturbation is given as

$$\delta \tau_{\rm g}= \frac{\displaystyle\int \vert \Psi_0\ver... ...ystyle\int \vert \Psi_0\vert^2 \omega^2 {\rm d} \omega} \cdot$$ (C.11)

Here we have assumed the amplitude perturbation to be small so that to first order ${\rm e}^{\Delta \gamma} \approx 1$.

We note that if one used only one travel-time parameter the equation obtained would be identical to this.

Phase travel-time perturbation

Close to the central frequency where the power spectrum has its maximum the term $(\omega - \omega _0) \delta \tau_{\rm g}$ should be small compared to $\omega_0 \delta \tau_{\rm p}$. We therefore neglect the first term and find that for the phase travel-time perturbation we get

$$\frac{\partial M}{\partial \delta \tau_{\rm p}} = -\int \vert... ..._0 \delta \tau_{\rm p} - \Delta \theta ) % {\rm d} \omega \; .$$ (C.12)

Expanding the sine term to first order gives

$$\int \vert \Psi_0\vert^2 {\rm e}^{\Delta \gamma} (\omega_0 \delta \tau_{\rm p} - \Delta \theta) {\rm d} \omega = 0\; .$$ (C.13)

From this we find that the phase travel-time perturbation is given as

$$\delta \tau_p= \frac{\displaystyle\int \vert \Psi_0\vert^2 \... ...yle\omega_0 \int \vert \Psi_0\vert^2 {\rm d} \omega} \; \cdot$$ (C.14)

Appendix D: Ray-based Green's functions

In this appendix we wish to solve the wave equation using the ray approximations. This will allow us to derive ray-based Green's functions for the reference medium which can be used to calculate sensitivity kernels for the perturbed medium which takes diffraction into account. All quantities in this appendix are related to the reference model thus the subscript 0 has been obmitted for ease of notation. For a more detailed treatment see e.g. Gough (1993). First we start by writing the wave field in terms of a slowly varying amplitude and a phase component as

$$\Psi= A {\rm e}^{i\theta} \; .$$ (D.1)

Inserting this into the wave equation gives the two equations as

$$\nabla \theta \cdot \nabla \theta = \frac{\dot{\theta}^2 - \omega_{\rm ac}^2}{c^2} + \frac{\nabla^2 A}{A}$$ (D.2)

$$A \nabla ^2 \theta + 2 \nabla A \nabla \theta + \frac{\partial}{\partial t} \left ( A^2 \dot{\theta} \right )= 0.$$ (D.3)

The first of these equations is known as the eikonal equation while the second is the transport equation. We then define the local frequency and wave number in terms of the phase as

$$\frac{\partial \theta}{\partial t}=- \omega \; , \quad \nabla \theta = {\vec k} \; .$$ (D.4)

Since the amplitude is assumed slowly varying the amplitude term in the eikonal equation can be neglected. Inserting Eq. (D.4) into the eikonal equation gives the dispersion relation between frequency and wave number

$$k^2 = \frac{\omega^2 - \omega_{\rm ac}^2}{c^2}$$ (D.5)

For a stationary medium the frequency $\omega$ is constant and the phase of the wave field is given as

$$\theta = \int_\Gamma {\vec k} \cdot {\rm d}{\vec r} - \omega t$$ (D.6)

where the integral is taken along the ray path, $\Gamma$, which is given by

$$\frac{{\rm d} {\vec r}}{{\rm d}t} = \frac{\partial \omega}{\partial {\vec k}} \cdot$$ (D.7)

Here $\partial \omega / \partial {\vec k}$ is the group speed. This is the speed at which the energy of the wave field propagates. The spatial varying part of the wave field can be written as

$$\Psi = A {\rm e}^{i \omega \tau _{\rm p}}\; ,$$ (D.8)

where $\tau_{\rm p}$ is the phase travel time which is given as

$$\tau _{\rm p} = \int _\Gamma \frac{{\vec k}}{\omega} \cdot {\rm d}{\vec r}\; .$$ (D.9)

Here $\omega / {\vec k}= {\vec v}_{\rm p}$ is the phase velocity.

The travel-time perturbation can be calculated using Fermats Principle which states that to first order a perturbation to ${\vec k}$does not perturb the ray path (see e.g. Gough 1993). The the travel-time perturbation is then given as

$$\delta \tau _{\rm p} = \frac{1}{\omega} \int _{\Gamma_0} \delta {\vec k} \cdot {\rm d}{\vec r}\;$$ (D.10)

where $\delta {\vec k} = {\vec k} - {\vec k}_0$ is the wave number perturbation, ${\vec k}_0$ is the wave number in the reference medium and the ray path $\Gamma_0$ is calculated in the reference medium. Equation (D.10) is the base for doing inversion of time-distance data in the ray approximation.

In a stationary medium the amplitude is also independent of time, thus the last term in the transport equation is zero. The remaining terms can with the help of Eq. (D.4) be rewritten as

$$\nabla \cdot \left (A^2 {\vec k} \right ) = 0.$$ (D.11)

This is a conservation equation for the amplitude. Using Gauss theorem one can see that this equation implies that for a ray tube

$$A^2k{\rm d}S = \mbox{constant}$$ (D.12)

where ${\rm d}S$ is the cross section of the ray tube. Thus we find that

$$A^2({\vec r_2}) = A^2({\vec r_1})\frac{v_{\rm p}({\vec r_2})}{v_{\rm p}({\vec r_1})} \frac{{\rm d}S_1}{{\rm d}S_2} \; ,$$ (D.13)

where ${\vec r}_1$ and ${\vec r}_2$ denotes the position of the two cross sections and v_p is the phase speed. Equation (D.13) describes the change in amplitudes in terms of the geometrical development of the ray tube. If we approximate the development of an angular area element ${\rm d}S$ as ${\rm d}S = r^2 {\rm d} \theta \approx (v_{\rm p} \tau_{\rm p}) ^2 {\rm d} \theta$, where ${\rm d} \theta$ is the angular area on the unit circle and r is the length along the ray path, we find that

$$A({\vec r}) = A({\vec r^{\prime}}) \frac{\tau_{\rm p}({\vec r^{\prime}})}{\tau_{\rm p}({\vec r})} \; ,$$ (D.14)

where we have taken ${\vec r}_1$ to be the source position ${\vec r}^\prime$. One complication arises because the amplitude close to the source should correspond to the homogeneous case and according to Eq. (17) the amplitude term A is divergent at the source. It turns out that $A \tau_{\rm p}$ is a bounded quantity at the source and thus we find that

$$A({\vec r}) \propto \frac{1} {\tau_{\rm p}({\vec r})} \; \cdot$$ (D.15)

This means that the amplitude is proportional to the reciprocal of the travel time. The same approximation is obtained if one assumes the sound speed to be constant. Using this approximation it is possible to obtain information about the behavior of the amplitude from the travel time without further calculations.

To describe the geometrical behavior of a ray tube in a more quantitative way one can use the Jacobian of the problem (Jensen et al. 1994). In three dimensions we have the Jacobian given by:

$$J = \left \vert \frac{\partial {\vec r} }{\partial (l,\Theta,... ...\frac{\partial z}{\partial \psi} \end{array} \right \vert \; ,$$ (D.16)

where $\Theta$ and $\psi$ are the declination and azimuthal of the take-off angle. Equation (D.13) can be expressed in terms of the Jacobian as

$$A ({\vec r})= A ({\vec r^\prime}) \left \vert \frac{v_{\rm p}... ... p}({\vec r^\prime}) J({\vec r})} \right \vert ^{1/2} \; \cdot$$ (D.17)

Here ${\vec r}^\prime$ again denotes the source position. As before the amplitude is divergent at the source thus posing a problem. Jensen et al. (1994) showed that $A^2 ({\vec r^\prime}) J({\vec r^\prime})$ is a bounded quantity and gave the solution for a cylinder symmetric problem as

$$A ({\vec r}\vert{\vec r^\prime} )= \frac{1}{4 \pi} \left \ver... ... p}({\vec r^\prime}) J({\vec r})} \right \vert ^{1/2} \; \cdot$$ (D.18)

Close to the source this expression corresponds to the amplitude for a homogeneous medium. Since the homogeneous Green's function is integrable at the source this ensures that the ray based Green's function also is integrable at the source.