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
are filtered out in order to remove granulation and supergranulation noise. The data are further
filtered to select parts of the wave propagation diagram.
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., 1997
; Kosovichev 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 (
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, 1984; Duvall and Harvey, 1986; Libbrecht, 1992; Schou, 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, 1996
; D’Silva et al., 1996; D’Silva, 1996; Kosovichev 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., 2001; Birch 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, 1996; Jensen 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
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.