### 4.3 Time-distance helioseismology

The purpose of time-distance helioseismology (Duvall Jr et al., 1993) is to measure and interpret the travel times of solar waves between any two locations on the solar surface. A travel time anomaly contains the seismic signature of buried inhomogeneities within the proximity of the ray path that connects two surface locations. An inverse problem must then be solved to infer the local structure and dynamics of the solar interior (see, e.g., Jensen, 2003, and references therein).

#### 4.3.1 Fourier filtering

The first operation in time-distance helioseismology is to track Doppler images at a constant angular velocity to remove the main component of solar rotation, as is done in ring-diagram analysis (see Section 4.2.1). Postel’s azimuthal equidistant projection is often used in time-distance helioseismology. The resulting data cube is Fourier transformed to obtain .

A filtering procedure is then applied to the data (Duvall Jr et al., 1997). First, frequencies below 1.5 mHz are filtered out in order to remove granulation and supergranulation noise. The data are further filtered to select parts of the wave propagation diagram.

 Table 1: Parameters of Fourier filters used in p-mode time-distance helioseismology. The Gaussian phase-speed filters are defined by . Courtesy of T.L. Duvall and S. Couvidat.
 distance range (Mm) ( km s–1) ( km s–1) 1 3.7 – 8.7 12.8 2.6 2 6.2 – 11.2 14.9 2.6 3 8.7 – 14.5 17.5 2.6 4 14.5 – 19.4 25.8 3.9 5 19.4 – 29.3 35.5 5.3 6 26.0 – 35.1 39.7 3.0 7 31.8 – 41.7 43.3 3.2 8 38.4 – 47.5 47.7 3.6 9 44.2 – 54.1 52.3 4.5 10 50.8 – 59.9 57.2 3.8 11 56.6 – 66.7 61.1 3.4

In the case of p modes, a phase speed filter is applied to the data, and the f-mode ridge is filtered out. This choice of filter is based on the fact that acoustic waves with the same horizontal phase speed travel the same horizontal distance (see Bogdan, 1997). Thus, to measure the travel time for acoustic waves propagating between two surface locations, it is appropriate to consider only those waves with the same phase speed. The choice of the phase speed depends on the travel distance. A list of Gaussian phase-speed filters,

is provided in Table 1 for various distances. Here, is the mean phase speed, is the dispersion, and is the wavenumber. The filtered signal is given by . The larger the phase speed, the deeper p-mode wavepackets penetrate inside the Sun.

A separate filtering procedure is applied for surface gravity waves, which are used to probe the near surface layer. In this case the filter function is 1 if and 0 otherwise. The parameter controls the width of the region around the f-mode dispersion relation . A reasonable choice is . This value allows for large Doppler frequency shifts introduced by flows, and does not let the ridge through.

There is some freedom in the selection of the Fourier filters. For example one may construct filters that depend explicitly on the direction of the wavevector , and not simply on the wavenumber (see, e.g., Giles, 1999).

#### 4.3.2 Cross-covariance functions

The basic computation in time-distance helioseismology is the temporal cross-covariance between filtered signals at two points and on the solar surface,

where is the temporal sampling and is the time duration of the observation. Times are evaluated at discrete values . Multiplication by the temporal window function is already included in the definition of , and the sum over is a short notation to mean the sum over all discrete times in the interval . The normalization factor is a correction that becomes significant only when is small. We note that in practice it is faster to compute the time convolution in the Fourier domain. The positive time lags () give information about waves moving from to , and the negative time lags about waves moving in the opposite direction.

The cross-correlation is useful as it is a phase-coherent average of inherently random oscillations. It can be seen as a solar seismogram, providing information about travel times, amplitudes, and the shape of the wave packets traveling between any two points on the solar surface. Figure 13 shows a theoretical cross-covariance as a function of distance between and , and time lag . This time-distance diagram was computed for a spherically symmetric solar model (see Sekii and Shibahashi, 2003). The first ridge corresponds to acoustic waves propagating between the two points without additional reflection from the solar surface. The next ridge corresponds to waves which arrive after one reflection from the surface, and the ridges at greater time delays result from waves arriving after multiple bounces. The backward branch associated with the second ridge corresponds to waves reflected on the far-side of the Sun. In most applications, only the direct (first-bounce) travel times are measured. Figure 14 shows the p-mode cross-covariance function at short distances, computed with the Fourier filters given in Table 1. In this example, cross-covariance functions were averaged over pairs of points at fixed distance to reduce noise. An example cross-covariance function for the f-mode ridge is plotted in Figure 15.

The travel times of the wave packets are measured from the (first-bounce) cross-covariance function. Local inhomogeneities in the Sun will affect travel times differently depending on the type of perturbation. For example, temperature perturbations and flow perturbations have very different signatures. Given two points and on the solar surface the travel time perturbation due to a temperature anomaly is, in general, independent of the direction of propagation between and . However, a flow with a component directed along the direction will break the symmetry in travel time for waves propagating in opposite directions: Waves move faster along the flow than against it. Magnetic fields introduce a wave speed anisotropy and will have yet another travel time signature (this has not been detected yet).

#### 4.3.3 Travel time measurements

Cross-covariance functions can have a large amount of realization noise due to the schochastic nature of solar oscillations, and it has proved difficult to measure travel times between two individual pixels on the solar surface. Some averaging is often required. The time over which the cross-covariances are computed is usually (this puts a limit on the temporal scale of the solar phenomena that can be studied). In order to further enhance the signal-to-noise ratio, Duvall Jr et al. (1993) and Duvall Jr et al. (1997) suggested to average over points that belong to an annulus or quadrants centered at . For instance, to measure flows in the east-west direction in the neighborhood of point , cross-covariances are averaged over two quadrant arcs and that include points a distance from :

Figure 16 shows the geometry of the averaging procedure. An example cross-covariance is shown in Figure 17 in the f-mode case. Correlations at positive times () correspond to waves that propagate westward, and correlations at to waves that propagate eastward. As shown in Figure 17, cross-covariances can be further averaged along lines of constant phase within a range of distances. Likewise an average cross-correlation is constructed from south-north quadrants to measure flows in the meridional direction. Another average is obtained by averaging over the whole annulus,
where is an annulus of radius centered at . This average can be used to measure separately the waves that propagate inward and outward from the central point.

At fixed and , the cross-covariance function oscillates around two characteristic (first-bounce) times . Center-to-annulus or center-to-quadrants travel times are often measured by fitting Gaussian wavelets (Duvall Jr et al., 1997Kosovichev and Duvall Jr, 1997). This procedure distinguishes between group and phase travel times by allowing both the envelope and the phase of the wavelet to vary independently. The positive-time part of the cross-correlation is fitted with a function of the form

where all parameters are free, and the negative-time part of the cross-correlation is fitted separately with
The times and are the so-called phase travel times. The basic observations in time-distance helioseismology are the travel time maps and , measured for each of the three averaged cross-correlations , , and . The travel time differences are mostly sensitive to flows while the mean travel times are sensitive to wave-speed perturbations. Maps of measured travel times are shown in Figure 18: Most of the signal in these maps is due to supergranular flows (15 – 30 Mm length scales). We note that an alternative definition of travel time can be obtained by fitting a model cross-covariance function to the data (Gizon and Birch, 2002). This last definition is often used in geoseismology (see, e.g., Zhao and Jordan, 1998).

In order to maximize the potential resolution of time-distance helioseismology it is desirable to obtain travel times from cross-covariances measured with shorter and with as little spatial averaging as possible. However conventional fitting methods will fail when the cross-covariance is too noisy. A new robust definition of travel time was introduced by Gizon and Birch (2004) to measure travel times between individual pixels and as short as . According to this definition, the point-to-point travel time for waves going from to , denoted by , and the travel time for waves going from to , denoted by , are given by

with the weight functions given by
In this expression is a (smooth) reference cross-covariance function computed from a solar model and . The window function is a one-sided function (zero for negative) used to separately examine the positive- and negative-time parts of the cross-correlation. The window is used to measure , and is used to measure . A standard choice is a window that isolates the first-skip branch of the cross-covariance.

This definition has a number of useful properties. First, it is very robust with respect to noise. The fit reduces to a simple sum that can always be evaluated whatever the level of noise. Second, it is linear in the cross-covariance. As a consequence, averaging various travel time measurements is equivalent to measuring a travel time on the average cross-covariance. This is unlike previous definitions of travel time that involve non-linear fitting procedures. Third, the probability density functions of and are unimodal Gaussian distributions. This means, in particular, that it makes sense to associate an error to a travel time measurement. The Born sensitivity kernels discussed in Section 4.3.5 were derived according to this definition of travel times.

#### 4.3.4 Noise estimation

In global helioseismology, it is well understood that the precision of the measurement of the pulsation frequencies is affected by realization noise resulting from the stochastic nature of the excitation of solar oscillations (see, e.g., Woodard, 1984Duvall Jr and Harvey, 1986Libbrecht, 1992Schou, 1992). It is important to study these properties since the presence of noise affects the interpretation of travel time data. In particular, the correlations in the travel times must be taken into account in the inversion procedure.

An interesting approach, pioneered by Jensen et al. (2003a), consists of estimating the noise directly from the data by measuring the rms travel time within a quiet Sun region. The underlying assumptions are that the fluctuations in the travel times are dominated by noise, not by ’real’ solar signals, and that the travel times measured at different locations can be seen as different realizations of the same random process. By ’real’ solar signals we mean travel time perturbations due to inhomogeneities in the solar interior that are slowly varying over the time of the observations. Jensen et al. (2003a) studied the correlation between the center-to-annulus travel times as a function of the distance between the central points, at fixed annulus radius.

Gizon and Birch (2004) derived a simple model for the full covariance matrix of the travel time measurements. This model depends only on the expectation value of the filtered power spectrum and assumes that solar oscillations are stationary and homogeneous on the solar surface. The validity of the model is confirmed through comparison with MDI measurements in a quiet Sun region. Gizon and Birch (2004) showed that the correlation length of the noise in the travel times is about half the dominant wavelength of the filtered power spectrum. The signal-to-noise ratio in quiet-Sun travel time maps increases roughly like the square root of the observation time and is maximum for a distance near half the length scale of supergranulation.

#### 4.3.5 Travel time sensitivity kernels

In this section we show example computations of travel time sensitivity kernels. Kernels are the functions that connect small changes in the solar model with small changes in travel times

The sum over is taken over all possible relevant types of perturbations to the model, for example local changes in sound speed, density, flows, magnetic field, or source properties. For each type of perturbation there is a corresponding kernel that depends on , the observation points and that the travel time is measured between, and a spatial location that ranges over the entirety of the solar interior. In this section we will describe the various types of approximate calculations of the that have been done.

Notice that kernels can be computed for quantities other than travel times. For example, in geophysics Dahlen and Baig (2002) computed the first-order sensitivity of the amplitude of the observed waveform to a small local change in sound speed. In particular it might be helpful to compute kernels for the mean frequency and amplitude of the cross-correlations, with the aim of using these quantities to help constrain inversions.

##### 4.3.5.1 Ray approximation.
The first efforts at computing the sensitivity of travel times to changes in the solar model (see, e.g., Kosovichev, 1996D’Silva et al., 1996D’Silva, 1996Kosovichev and Duvall Jr, 1997) were based on the ray approximation. In the ray approximation the travel time perturbation is approximated as an integral along the ray path. Fermat’s principle (see, e.g., Gough, 1993) says that in order to compute the first order change in the travel time we do not need to compute the first order change in the ray path, i.e., we can simply take the integral over the unperturbed ray path. In particular, we have (Kosovichev and Duvall Jr, 1997):

where the ray path connecting and is denoted by . The increment of path length is and is a unit vector directed along the path, in the sense of going from to . The other quantities in Equation (64) are the change in the sound speed , the mass flow , the Alfvén velocity , the acoustic cut-off frequency , and the phase speed . Ray theory is based on the assumption that the perturbations to the model are smooth and that the wavepacket frequency bandwith is very large. Bogdan (1997) showed that the energy density of a realistic wavepacket was substantial away from the ray path. This result strongly suggested that perturbations located away from the ray path could have substantial effects on travel times. It is now well known that ray theory fails when applied to perturbations that are smaller that the first Fresnel zone (see, e.g., Hung et al., 2001Birch et al., 2001).

##### 4.3.5.2 Finite-wavelength kernels in the single-source approximation.
The three approximations that have been used to treat small scale perturbations for time-distance helioseismology are the Born approximation, the Rytov approximation, and the Fresnel zone approximation. Here we describe some results that have been obtained using these methods in the single-source approximation. The single-source approximation, which was motived by the “Claerbout Conjecture” (Rickett and Claerbout, 1999) and by analogy with previous ray-theory based work, says that the one-way travel time perturbation can be found by looking at the time-shift of a wavepacket created at one observation location and then observed at the other. Gizon and Birch (2002) showed that the single-source approximation ignores a potentially important scattering process; this issue is discussed in detail in the paragraph on distributed source models. Motivated by the use of the Born approximation to compute travel time sensitivities in the context of geophysics (see, e.g., Marquering et al., 1998), Birch and Kosovichev (2000) used the Born approximation to compute travel time kernels for time-distance helioseismology. Figure 19 shows an example result. The shape of the kernel is the characteristic “banana-doughnut” shape first described by Marquering et al. (1998). The kernel is hollow along the ray path, travel times are not sensitive to small-scale sound-speed perturbations located along the ray path. This has been explained in terms of wavefront healing (Nolet and Dahlen, 2000). The ringing in the kernel with distance from the ray path is a result of the finite bandwidth of the wavefield.

Another approach to producing finite-wavelength kernels is the Fresnel zone approximation (see, e.g., Sneider and Lomax, 1996Jensen et al., 2000). The Fresnel zone approach makes a simple parametrized approximation to the kernels. In particular, in the Fresnel zone approximation, the travel time kernels, for fixed and , are written as (Jensen et al., 2001)

where is the ray-theoretical travel time from to , and . The parameter is the central frequency of the wavepacket, is the corresponding period, and is related to the frequency bandwidth of the wavepacket. The amplitude of the kernel is determined by demanding that the total integral of the kernel be equal to , which is the ray theory value for the total integral of the kernel. Notice that the value of the kernel goes to infinity as tends to or . Jensen et al. (2000) showed, by comparison with direct simulation, that a two dimensional analogue of the Fresnel zone kernel (Equation (65)) was a reasonable approximation to the actual linear sensitivity of travel time to changes in the sound speed.

The Rytov approximation gives a first order correction to the phase of the wavefield, instead of a first order correction to the wavefield itself as one obtains from the Born approximation. Jensen and Pijpers (2003) used the Rytov approximation to write travel time kernels for the effects of sound-speed perturbations. The results were qualitatively similar to the Born approximation results of Birch and Kosovichev (2000). No comparison has been made, in the helioseismology literature, of the ranges of validity of the Rytov and Born approximations.

##### 4.3.5.3 Distributed source models.
Gizon and Birch (2002), motived by Woodard (1997), gave the first comprehensive treatment of the linear forward problem for time-distance helioseismology. Such a treatment must include a physical model of wave excitation which takes into account the fact that waves are stochastically excited by sources (granules) that are distributed over the solar surface (see Section 3.2). An example calculation (Figure 20) demonstrates that the single-source approximation is qualitatively incorrect, as it ignores a scattering process that can be dominant for some perturbations. Figure 20 compares an f-mode travel time kernel for local changes in the wave damping rate, computed in the single-source and distributed-source models. The hyperbolic feature in the distributed-source kernel is not present in the single-source model. As discussed by Gizon and Birch (2002), the single-source approximation ignores the scattered wave that is seen at the observation point that is treated as a source in the single-source approximation. Birch et al. (2004) used the recipe derived by Gizon and Birch (2002) to obtain the linear sensitivity of p-mode travel times to changes in sound speed. The results showed that the sensitivity depends on the filtering that is done to the wavefield before the travel times are measured. Figure 21 shows sound-speed kernels computed for three different cases: a model with no filters, a model with an MTF (see Section 3.4) roughly like the one for MDI full-disk data, and a model with an MTF and also a phase-speed filter of the type employed in routine time-distance analysis. The three kernels are all strikingly different. In the unfiltered case, the kernels looks much like the traditional “banana-doughnut” kernel (see Figure 19). With the introduction of the MTF, the travel time becomes more sensitive to deeper sound-speed perturbations and less sensitive to the near-surface perturbations. This is because the MTF preferentially removes high wavenumbers, and thus modes with low phase speeds and shallow turning points. With the introduction of the phase-speed filter the kernel again changes drastically. The phase-speed filter removes all of the waves with lower turning point much below the turning point of the target ray. As a result, the sensitivity much below the lower turning point of the ray is greatly reduced. This example emphasizes the importance of including the filtering in models of travel time sensitivity kernels. The strong dependence of the shape of the kernels on the filtering suggests that the kernels could to be tailored by adjusting the filters that are used in the data analysis procedure.

#### 4.3.6 Inversions of travel times

The inverse problem is to determine the perturbations to the solar model that are consistent with a particular set of observed travel times . For the inverse problem, the kernel functions and the noise covariance are in general assumed to be known.

Unlike in the ring-analysis case, where a set of one dimensional inversions are performed (see Section 4.2.3), the time-distance inversion procedure is three-dimensional. So far, only RLS inversions have been implemented in time-distance helioseismology. In the RLS approach, the minimization problem reads

where we have assumed, for simplicity of notation, a diagonal covariance matrix. The are the noise estimates for each of the travel times . The regularization parameter is and is the regularization operator. Notice that we have not assumed any particular parametrization of the perturbations to the solar model . The ideal choice of regularization operator and the optimal method for choosing the regularization parameter are currently open questions.

The minimization of , Equation (66), has been done using either the LSQR algorithm (Paige and Saunders, 1982), the multi-channel deconvolution (MCD; Jacobsen et al., 1999), or via singular value decomposition (Hughes and Thompson, 2003). The inputs to either of these methods are the travel time kernels, which result from the linear forward problem (see Section 3.5), the regularization operator, the choice of regularization parameter, and in general also the covariance matrix of the noise.

The central question about inversions is the degree to which any particular inversion method is able to retrieve the true Sun, e.g., the actual sound-speed variations and flows in the interior, from a set of observed travel times. The degree to which any inversion method succeeds will presumably depend on a number of things: (i) the accuracy of the forward model, (ii) the noise level in the travel times, (iii) the depths and spatial scales of the real variations in the Sun, (iv) the number and type of travel times that are available as input to the inversion, and (v) the accuracy of the travel time noise covariance matrix. The two main approaches that have been taken to study inversion methods are either to invert artificial data or to invert real data for perturbations that are relatively well known from global helioseismology or from surface measurements. It is also useful to compare the results of different inversion methods applied to the same set of observed travel times.

The approach of inverting real data is appealing as it is by definition “realistic”. This approach for time-distance was pioneered by Giles (1999) who did ray theory based inversions to measure differential rotation. The results were in approximate agreement with the results of global model inversions, which validated time-distance inversions applied to slowly-varying large scale flows. Gizon et al. (2000) inverted f-mode travel times for horizontal flows in the near-surface. Figure 22 shows a comparison between the projection of the inferred flows onto the line of sight and the directly observed line-of-sight Doppler velocity. The correlation coefficient between the two is 0.7 (Gizon et al., 2000), which shows that their inversion method, an iterative deconvolution, is able to approximately retrieve the horizontal flow velocities from the travel times.

There have been a number of efforts to validate inversions by the generation and inversion of artificial data. These tests are useful as they provide an intuitive understanding of various regularization schemes, the choice of regularization parameter, the effects of limited sets of input travel times, the effects of incorrect assumptions regarding the noise covariance, the potential resolution of any particular inversion, and the cross-talk between different model parameters. A drawback to testing inversions with artificial data is that one always wonders about the realism of the artificial data. Hopefully, in the not distant future quite realistic data will be available for the purpose of testing inversion schemes.

Couvidat et al. (2004) compared inversions done using ray-approximation kernels with inversions done with Fresnel-zone kernels. The results were quite similar for depths between the surface and the lower turning point of the deepest rays used in the inversion. Inversions done with the ray kernels cannot retrieve sound speed perturbations below the lower turning point of the deepest ray, while the Fresnel-zone kernels extend below the depth of the deepest ray. The results of Couvidat et al. (2004), as well as those of Giles (1999), suggest that even though the ray approximation overestimates travel times for small scale perturbations this effect does not seriously corrupt the large scale features found by inversion.

Another test of inversion methods was done by Zhao and Kosovichev (2003a). The two main results of this study, regarding inversions methods, were that there is cross-talk between upflows and convergent flows, and between downflows and divergent flows, and that the MCD and LSQR methods give essentially the same results when applied to quiet sun MDI data. Figure 23 shows an example of a test of the LSQR method. The initial model is shown in the top panel. Travel times were generated using the ray approximation. The travel times were then inverted using five iterations of LSQR to obtain the flow field shown in the second panel of Figure 23. Notice that there are small vertical flows near the surface at the centers and boundaries of the supergranulations cells. After 100 iterations of LSQR, these small vertical velocity features are removed (bottom panel of Figure 23).

The first end-to-end test of time distance helioseismology using artificial data was done by Jensen et al. (2003b). In this important study, artificial data were generated by numerical wave propagation through a three-dimensional stratified and horizontally inhomogeneous model. Waves were generated by stochastic sources distributed over a layer 50 km below the upper boundary of the surface of the model. The artificial wave field was “observed” at the spatial resolution and temporal cadence of MDI high-resolution observations. This artificial data set was then subjected to standard time-distance analysis, and the inversion was performed with the Fresnel-zone kernels described by Jensen and Pijpers (2003). Figure 24 shows a comparison of the original sound-speed model and the result of the inversion of the artificial data. Notice that the inversion recovers much of the original structure. The inversion result contains noise, which is to be expected from travel times computed from finite time series. The only noise source in this example is realization noise, noise resulting from the stochastic nature of the wave excitation.