and started routine observations in May 1996. MDI produces 10242 images

mostly with a cadence of one per minute. The MDI instrument can operate

in one of several modes generating various combinations of Doppler images,

intensity images and magnetograms in full disk and higher resolution. In the

full-disk mode the whole visible disk of the Sun is observed while a (¬xed)

area is imaged with a magni¬cation of 3.2 in the high-resolution mode.

While being located above the disturbing in¬‚uence of the atmosphere

at the Earth-Sun L1 Lagrange point means that the Sun can be observed

continuously and with uniform quality, the telemetry is limited and only a

subset of the images can be transmitted. For two to three months per year

and roughly 8 hours per day a high telemetry rate is available and, depending

on the observing program, Dopplergrams are available each minute.

In order to provide continuous data for the whole year a number of lower-

resolution data products are produced, stored on board the spacecraft and

downlinked several times per day. Most important for helioseismology is

the so-called Medium-l program, in which the velocity images are convolved

with a Gaussian, cropped at 90% of the solar radius and sampled at every

¬fth pixel to provided roughly 2002 pixel images. The temporal coverage

Schou

250

ranges from 95% to well over 99% on time scales of months. The time series

from the Medium-l program are produced up to l = 300.

17.3.3 Other projects

A number of other projects have also been observing solar oscillations over

the years and several are still operating, including ECHO, TON, BISON

and IRIS all of which are ground-based networks and VIRGO and GOLF on-

board the SOHO spacecraft. The observations from several of these projects

have been important for the analysis of, in particular, the solar core.

17.4 Normal mode analysis

The analysis of the basic image data is generally done in a standard series

of steps. First the images are calibrated, a spherical-harmonic transform

(SHT) is performed, the resulting time series are Fourier transformed in time

(and sometimes turned into power spectra) and the spectra are analyzed to

estimate the various mode parameters. Each of these steps is described

below followed by a discussion of some of the problems encountered in the

analysis. For more details see Schou (1998) and Schou et al. (2002).

After this basic analysis an inverse procedure is typically applied to infer

the properties of the solar interior, see Sekii (2003).

17.4.1 Time series generation

The ¬rst step in the analysis is to calibrate the basic images, typically Dopp-

lergrams. Also basic image parameters such as image radius and center need

to be determined.

The next step is to apodize and remap the images to a uniform grid in

longitude and sin θ. The apodization, which is most often a cosine taper

in fractional radius, is done to avoid sharp edges and the resulting ringing

in the subsequent transform. The remapping takes into account the image

center and apparent solar radius, the P- and B-angles, the ¬nite observer-

Sun distance and the non-uniform angular speed of the observer as seen

from the Sun. (The P-angle is the position angle of the solar pole and the

B-angle is the latitude of the sub-observer point.) The non-uniform angular

speed of the observer is important since assuming a constant speed would

lead to an apparent annual variation in the observed rotation rate.

The SHT is generally done by calculating inner products of the individ-

ual remapped images and the spherical harmonics of the target modes. In

Helioseismic data analysis 251

practice this is performed (Brown 1985) by using a Fourier transform in

longitude and multiplying by the Plm ™s in latitude. This multiplication is

quite time consuming and while faster algorithms (Mohlenkamp 1999) do

exist, they are not used, mainly due to their complexity and the fact that

the subsequent peakbagging, rather than the SHT, is the limiting factor in

the overall data analysis.

The Fourier transform of the time series and the construction of power

spectra is in principle straightforward. However, the time series typically

have to be detrended and gap-¬lled in order to remove long-term drifts

and gaps. In the MDI case the detrending is done by subtracting low-

order polynomials, while GONG uses ¬rst di¬erencing. GONG has also

investigated using multi-tapers to construct the power spectra (Komm et

al. 1999).

17.4.2 Peakbagging

The ¬nal and by far the most complex part of the data analysis is what

has become known as peakbagging and ridge ¬tting where the mode param-

eters such as frequencies, linewidths and amplitudes are determined from

the Fourier transforms or power spectra. This complexity is caused by a

combination of properties of the oscillations and the way they are observed.

First of all the fact that we can only observe less than half of the solar

surface means that, although the spherical harmonics form an orthonormal

basis on a sphere, the SHT is not able to separate completely the individual

modes. This leads to what is known as leaks in the time series and power

spectra where, in addition to the target mode, several modes with nearby l

and m appear. To be speci¬c let x be the time series of the oscillations on

the Sun and y the observed time series, then

ylm = clml m xl ,m ,

l ,m

where clml m is the so-called leakage matrix (see Schou & Brown 1994).

Furthermore, the power spectrum of a single damped and stochastically

excited oscillator (such as a single mode) is given by the product of a Lorent-

zian and a random function with an exponential distribution. The width of

the Lorentzian is given by the inverse damping time while the integral under

the pro¬le is given by the average power in the mode. For typical modes in

the middle of the p-mode band the linewidths are of the order 1 µHz.

The combination of the above two issues and the fact that the frequency

spacings between adjacent modes are often comparable to the linewidths

Schou

252

Fig. 17.1. Examples of m-averaged spectra for l = 150 using 72 days of MDI data.

From left to right is n = 0, n = 4 and n = 8. The individual spectra were shifted

by the measured splittings before being averaged. Notice the di¬erent scales on the

frequency axis. On the left panel the peak around ∆ν = 0 is the target mode, the

pairs around ±4 µHz are for l = ±1 (the splitting is caused by m ’ m = ±1) and

the one at 8 µHz is from l = 2 (the split components are m ’ m = ’2, 0 and 2.

means that the ¬tting methods need to include assumptions about the fre-

quency variations, simultaneous ¬ts of di¬erent modes, and/or elaborate

modeling of the e¬ects of the atmosphere, instrument and the SHT. These

problems are illustrated in Figure 17.1, which shows power spectra in the

three main regimes: low frequency where the m leaks (those from the same l

but di¬erent m) can be separated, moderate l and ν where the l leaks (those

from di¬erent l) can be separated but the m leaks blend together and the

ridge ¬tting regime where both l and m leaks blend together.

In the low-frequency regime simple algorithms should work well were it

not for the fact that the signal-to-noise ratio (S/N) is quite low. Also the

linewidth is sometimes less than the frequency resolution (i.e., the lifetime

is longer than the observation time) and thus any given mode may or may

not be excited.

Most of the normal analysis has so far been in the regime where the

m leaks blend together but where the l leaks can still be separated. The

methods developed here have also been used in the low-frequency regime.

At present two main algorithms are being used: the MDI algorithm (Schou

1992) and the GONG algorithm (Anderson et al. 1990; Hill et al. 1996).

These algorithms and some of their associated problems have been described

by Schou et al. (2002).

17.4.2.1 The MDI algorithm

In the MDI algorithm all m values at a given (l, n) are ¬tted simultaneously

and a leakage matrix calculated from our knowledge of the Sun, the instru-

ment and the data analysis is used together with a maximum likelihood

Helioseismic data analysis 253

minimization which takes into account the correlation between the modes.

In order to take the leaks and the associated correlations into account the

Fourier transforms rather than the power spectra are ¬tted. The frequency

interval used is generally quite small (typically 5 linewidths) in order to

speed up the algorithm. The l leaks are calculated from the leakage matrix