Up: The art of fitting
Subsections
It is well known that p modes are stochastically excited oscillators
(Kumar et al. 1988). The source of excitation lies in the many
granules covering the Sun. The modes are assumed to be independently excited
provided that their spatial scale is larger than the granule size
(Chang 1996). From the equation of an oscillator, the
statistics of the p-mode profile can be derived as:
| |
(9) |
where t is the time, x is the displacement, is the damping
term or the linewidth,
is the frequency of the mode and F(t) is the forcing
function. From this equation the Fourier transform of x can be
written as:
| |
(10) |
where and are the Fourier transform
of x(t) and F(t). From the large number of granules, it can be derived that the forcing function
is normally distributed. Therefore the 2 components (the real and imaginary
parts) of the Fourier transform of the forcing function are also normally distributed. For the p
modes, each component of the Fourier transform is normally distributed
with a mean of zero, and a variance given by:
| |
(11) |
The square of the modulus of , or
power spectrum, has a with 2 degree of freedom statistics
and its mean is twice that of Eq. (11).
This is the p-mode profile that is usually approximated by a
Lorentzian profile. Similarly other effects such as asymmetry can be
introduced in the profile of Eq. (11).
Instruments integrating over the solar surface the velocity or the
intensity signal observe a superposition of various modes of
different degrees. They are mainly sensitive to the low-degree
modes (). For a
given l, they detect a mixing of azimuthal order m for which a
visibility is prescribed (Toutain & Gouttebroze 1994;
Christensen-Dalsgaard & Gough 1982). Most often they can only
detect modes for which l+m is even. Since the Fourier components of the
observed time series have a normal
distribution, and since the different m are uncorrelated, the statistics of the power
spectra of unresolved observation is a with 2 degrees of
freedom. Toutain & Appourchaux (1994) gave an analysis of
the problem associated with these observations; we will not repeat it
here.
When the solar image is resolved many more degrees can be detected
making the data analysis somewhat more complicated. In order to
extract a single l,m mode from resolved observations, one has to apply a
specific spatial filter or weight to the velocity or intensity
images. Most often these weights are such that imperfect
isolation of the l,m mode is achieved; especially because the most
commonly used filters (spherical harmonics) are not orthogonal over
half a sphere. This leads to the existence of
other modes in the Fourier spectrum generated for a given l,m filter.
Therefore, the observed
Fourier spectrum is a linear combination of the modes
to be detected. This linear combination of the modes can be understood
as modes leaking into each other spectrum: this is represented by the so-called
leakage matrix. These leakages will produce correlations
between the different Fourier spectra. These correlations will modify
the statistics of the Fourier spectra, such that their power spectra
cannot be described as a with 2 degrees of freedom.
Therefore the statistics of the 2l+1 power spectra of a given l cannot be derived
from the product of 2l+1 with 2 degrees of freedom as in
Appourchaux et al. (1995).
Nevertheless, the real and imaginary parts of the Fourier
spectra will still be normally distributed; in other words, the
Fourier spectra have a multi-normal distribution defined by a
covariance matrix. This
fact will be used to derive the statistics of the observation. The
covariance matrix is the sum of the noise and mode covariance
matrices, which are not necessarily the same. Last but not least,
the theoretical probability distribution has to be
generated using the previous covariance matrix.
In summary, to understand the statistics of resolved observation, one has
to follow four steps:
- Compute leakage matrices,
- Compute mode covariance matrices (related to the leakage),
- Compute noise covariance matrices,
- Generate the likelihood from the theoretical probability distribution.
Each step is described in detail hereafter.
Due to the spherical symmetry of the Sun, the most likely
weights to be used to isolate the modes are the spherical harmonics Yl,m. Here we
generalize the approach to any weight Wl,m. The result is the
observation of Fourier spectra that are related to what we
want to detect, i.e. the Fourier spectra of the individual Fourier
spectrum , by the so-called leakage matrix (Schou
1992; Schou & Brown 1994). The following expression
can be derived for as many different degrees as needed; for simplicity we
wrote it for 2 different degrees l,l' as:
| |
(12) |
where and are 2 complex vectors made
each of 2l+2l'+2 component: 2l+1
components for l, 2l'+1 components for l' and
is the leakage matrix of both l and
l'. The dimension of is . The
coefficient of the leakage matrix can be computed as:
| |
(13) |
with:
| |
(14) |
where , , * denotes the complex conjugate, are the angles
in a spherical coordinate system, is the integration domain,
is the generalized velocity or intensity perturbation
of the mode (l',m'). The are not necessarily the usual spherical harmonics Yl',m'. For instance, the horizontal component of the velocity perturbation, and the intensity perturbation due the distortion of the the surface by the modes are both expressed as derivative of spherical harmonics.
is an apodization function, nl,m is a sensitivity correction factor
associated with Wl,m. The ratio ensures that .The apodization function A is the product of 3 different function as:
| |
(15) |
is the natural apodization function due to the way the images
are obtained: for intensity this is the limb darkening (), and
for velocity the projection factor ().
is the data analysis apodization: for data re-mapped on the
Sun's surface it is unity; for
no re-mapping, it is the projection factor ().
is the artificial apodization that can take into account the
non-linear velocity (or intensity) response of the instrument over the
solar disk, or can help to reduce limb effects. Here we must point out
that the leakage matrix has a useful property such as:
| |
(16) |
It shows that is in general not hermitian nor symmetrical.
Nevertheless, when , it is possible with a
proper sensitivity factor correction of Wl,m to have
such a property. In this case the sensitivity correction is given by:
| |
(17) |
which is the "natural" normalization factor of the perturbation
. Of course in this latter case, we have:
| |
(18) |
Unfortunately, the leakage matrix does not always have
such a nice property, especially because .
This was the case for the ground-based Luminosity Oscillations Imager
(LOI) (Appourchaux et al. 1994) and for the GONG instrument
(Hill 1997, private communication). In both cases, this is not
produced by the observation techniques but by the data analysis techniques.
If the
weight functions Wl,m and the observed perturbations have the
same symmetry properties as the spherical harmonics Yl,m (or if
), the leakage
matrix is real as shown by Schou (1992). In addition the
leakage elements of are zero if l+m+l'+m'
is odd; this is the case when the Sun is not tilted with respect
to the observer's axis of reference (P = 0, B = 0). If the axes of reference
of Wl,m differ
from that of the Yl,m these 2 properties can be lost. For
instance, an incorrect orientation of the Sun axis with respect to
the detector axis could lead to a complex leakage matrix; or a Sun
seen at an angle give a real leakage matrix with non-zero elements with l+m+l'+m'
odd. This latter property has been used by Gizon et al. (1997)
to infer the inclination of the Sun's core.
Equation (13) is valid when the size of the pixel is small compared with the
spatial scale of the degree. When the pixels are larger, one should
write the following:
| |
(19) |
where the are given by:
| |
(20) |
where is the area defined by the i-th pixel and
wi(l,m) is the weight applied to the i-th pixel to extract
the l,m mode. Equation
(19) is the more general form used for the LOI
(Appourchaux & Andersen 1990). As a starting point, the
wi(l,m) can also be taken as the .
To compute the covariance of the complex vector as a real number we form
the vector defined as:
In absence of noise,
the covariance matrix of the vector can be generated using a
complex notation:
| |
(21) |
is a super matrix where and are the real and imaginary
parts of a complex matrix which elements are given by:
| |
(22) |
where is the variance of the l'',m'' mode which
profile is given by
Eq. (11), in which is a function of m. The real and
imaginary parts of Eq. (22)
will give respectively the covariance of the
real (or imaginary) part of , and the covariance between
the real and imaginary part of . It is
obvious from Eq. (22) that is hermitian.
Schou (1992)
gave an equation similar to Eq. (22) for a real leakage matrix
and for a single degree. Here we add a subtlety
to the formulation of Schou (1992), the matrix
can be decomposed as follows:
| |
(23) |
where is the transpose of a matrix. The elements of
and are given by:
| |
(24) |
| |
(25) |
We will see later on that this decomposition is of prime importance
for understanding the statistics of the observation.
Unfortunately, the observed vector include a noise
contribution. Due to the way the data are combined, the noises between the
different 2l+2l'+2 components of this vector are also correlated.
Schou (1992) gave the
correlation matrix when the filter used are spherical harmonics Yl,m.
A more general formulation can be written as:
| |
(26) |
is a super matrix where and are the
real and imaginary parts of the complex matrix
. The dimension of is
. Its elements are given by:
| |
(27) |
with
| |
(28) |
where a is an apodization function which characterizes through
how the noise varies over the solar
image, assuming that the noise is uncorrelated between different
points on the Sun; are defined in Eq. (15).
When the instrumental noise is low, a is derived from the characteristics
of the solar noise. The evaluation of is less
straightforward than that of because we need to
know a model of the solar noise. An easier way to understand the noise correlation
is to built the
ratio covariance matrix or "pseudo" noise leakage matrix
as:
| |
(29) |
Here we can see the similarity between and . In velocity,
the granulation noise is rather low at the center of the
disk and then increases towards the limb; the meso- and
super-granulation exhibits somewhat different or complementary
center-to-limb variations. In intensity, the
granulation noise is a function of the number of granules; the noise
is larger at the center of the disk and decreases slowly towards the
limb. In addition the solar noise in intensity has no contribution
from mesogranulation (Fröhlich et al. 1997), making the
spatial dependence of the noise almost independent of frequency across the
p-mode range. This is not the case in velocity where mesogranulation still
contributes to the noise in the p-mode range. Therefore in intensity the
apodization a is closer to A than in velocity, making the ratio
covariance matrix very close to the leakage
matrix . Although is
not mathematically useful, it is a matrix easy to visualize and understand
(See Part II). The ratio matrix has some properties of the leakage matrix
like being not necessarily hermitian. This is not the case of
which is hermitian by definition.
Again, when the size of the pixel is large compared with the
spatial scale of the degree, Eq. (27) is rewritten as
follows:
| |
(30) |
where bi is the variance of the noise of pixel i. Equation
(30) is the more general form used for the LOI.
The statistical distribution of the Fourier spectra or of the
vector is a multi-normal distribution. The
probability density is given by:
| |
(31) |
where d is the number of elements of , is a short notation for the following matrix: ; this is the matrix given
by the sum of the p-mode and noise covariance matrix; the p modes and
the noises are assumed to be independent of each other. The matrix
can also be built from sub-matrices as:
;
as a result is also hermitian. Equation
(31) is the most general formulation for any multi-normal
distribution with a given covariance matrix (Kendall &
Stuart 1967).
Using Eq. (31), we can write the likelihood L of an observation
of at N different frequencies
as given by:
| |
(32) |
We assumed that the frequency bins are independent of each other.
This is the case when the data have no gaps. For unresolved
observation having gap,
the expression of the likelihood becomes extremely complicated as shown by
Gabriel (1994). For resolved observation
having gaps, as for the LOWL data of Tomczyk et al. (1995), it
is impracticable to use the full formulation of the likelihood: Tomczyk
et al. (1995) used Eq. (32) as an approximation for fitting the
LOWL data.
In principle, given the observed vector , it is always possible in
the absence of noise to recover the vector . Due to the presence
of noise only a solution close to the ideal one can be found that
will minimize the correlation between the components. Provided that
the leakage matrix can be inverted, we have by analogy to Eq.
(12):
| |
(33) |
where .
Then we can write a similar equation for and
as:
| |
(34) |
where is defined as:
| |
(35) |
is a super matrix where and are the
real and imaginary parts of the complex matrix
. Using Eq.
(34)
to replace by in Eq. (32) we can rewrite
this latter as:
| |
(36) |
with given by:
| |
(37) |
We recognize in Eq. (36) the probability density of the vector
to a
constant (i.e. ). As a matter of fact, it is well
known that using a linear transformation similar to that of Eq.
(34) will produce the new covariance matrix of
as written in Eq. (37) (Davenport &
Root 1958). It can be easily shown using Eqs. (23) and
(37) that the matrix is diagonal and its element are
given by:
| |
(38) |
where l''=l or l' and . Therefore Eq. (37) is the sum of a diagonal matrix representing
the correlation between the p modes; and of a new noise
covariance matrix representing the correlation of the
components of the vector after the transformation of
Eq. (33). It means that has no
correlation due to the p modes as we could expect from Eq. (33):
the leakage matrix of is the identity matrix.
In summary, there is no gain in
fitting data for which the leakage matrix is the identity matrix: the 2 approaches are
identical . The main problem is really to know the leakage matrices,
not only theoretically but also experimentally: this is the subject of
the Part II.
It can be derived from Eq. (37) that it is also possible to
remove correlation due to the noise by replacing by a
proper matrix associated with . The derivation of
this matrix is given in Appendix A.
When a single degree is observed, it is quite simple to maximize the
likelihood of Eq. (32) using , or using
as in Eq. (36). For low degree and low frequency modes, this is possible for l=0,
2, 3. As soon as the mode linewidth increases, at high
frequencies, the assumption of a single degree is not valid anymore.
For example, l=0 and l=1 overlap with l=2 and l=3, respectively.
At high frequencies, the effect of the aliasing degree should be
taken into account.
For the other low degree modes, the likelihood
becomes somewhat more complicated. It is well known, that in the diagramme of l=1, there are
leaks coming from other degrees. The l=6 and l=9 modes overlap with the l=1, while the l=3 modes overlap only at higher frequencies when the linewidth is larger than about 5 Hz. In the diagramme of l=4, there are
leaks of l=7 and vice versa (Appourchaux et al. 1997). The
leaks have severe effects on determination of the p-mode parameters of the
l=1. When many degrees are overlapping, one should use Eq. (32)
using the covariance matrix for l and l'. Nevertheless, we do not
advice to do so for fitting the p modes; it has some severe computer speed
penalty. Instead we advice to clean the data by inverting the full leakage
matrix taking into account the effects of the various degrees on each
other, in a similar way to Eq. (33). For example for l=1, one
should compute the leakage matrix on a sub-space of degrees namely for l=1, 6 and 9. These 3 degrees interact strongly in the Fourier spectra. For l=4 and l=5 one should compute the leakage matrix on sub-spaces for l=4 and 7, and for l=5 and 8. The advantage of this direct cleaning is to replace the original aliasing degrees by new aliasing degrees which are further away, in frequency, from the modes to be fitted. This technique has been applied to the LOI and GONG data, and is developed in Part II.
Last but not least, when the signal-to-noise
ratio is high (i.e. we neglect in Eq.
(37)), the elements of the vector
are all independent of each other, leading to a statistical distribution which is
a product of with 2 degree of freedom. This is an
approximation which is useful and less incorrect that using this
statistics for the GONG data for the vector as in
Hill et al. (1996).
Up: The art of fitting
Copyright The European Southern Observatory (ESO)