next previous
Up: The art of fitting


3 The statistics of p-mode spectra

3.1 Single mode

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:

 d}t^{2}}+2\pi\gamma\frac{{\rm d}x}{{\rm
 \end{displaymath} (9)
where t is the time, x is the displacement, $\gamma$ is the damping term or the linewidth, $\nu_{0}$ 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:

 \end{displaymath} (10)
where $\tilde{x}(\nu)$ and $\tilde{F}(\nu)$ 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:  
 \end{displaymath} (11)
The square of the modulus of $\tilde{x}(\nu)$, or power spectrum, has a $\chi^{2}$ 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).

3.2 Unresolved observations

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 ($l \le 4$). 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 $\chi^{2}$ with 2 degrees of freedom. Toutain & Appourchaux (1994) gave an analysis of the problem associated with these observations; we will not repeat it here.

3.3 Resolved observations

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 $\chi^{2}$ 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 $\chi^{2}$ 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:

Each step is described in detail hereafter.

3.3.1 Leakage matrices

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 $y_{m}^{l}(\nu)$ that are related to what we want to detect, i.e. the Fourier spectra of the individual Fourier spectrum $x_{m'}^{l'}(\nu)$, 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:  
 \end{displaymath} (12)
where $\vec{x}(\nu)$ and $\vec{y}(\nu)$ are 2 complex vectors made each of 2l+2l'+2 component: 2l+1 components for l, 2l'+1 components for l' and $\tens{\cal{C}}^{(l,l')}$ is the leakage matrix of both l and l'. The dimension of $\tens{\cal{C}}^{(l,l')}$ is $(2l+2l'+2) 
\times (2l+2l'+2)$. The coefficient of the leakage matrix can be computed as:  
 \end{displaymath} (13)
 ...theta,\phi)A(\theta,\phi) \sin \theta {\rm
 d}\theta {\rm d}\phi
 \end{eqnarray} (14)
where $m=-l,\ldots,l$, $m'=-l',\ldots,l'$, * denotes the complex conjugate, $\theta,\phi$ are the angles in a spherical coordinate system, $\cal{D}$ is the integration domain, $\tilde{Y}_{l',m'}$ is the generalized velocity or intensity perturbation of the mode (l',m'). The $\tilde{Y}_{l',m'}$ 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. $A(\theta, \phi)$ is an apodization function, nl,m is a sensitivity correction factor associated with Wl,m. The ratio ensures that $\tens{\cal{C}}^{(l,l)}_{m,m}=1$.The apodization function A is the product of 3 different function as:  
A(\theta,\phi)=A_{\rm n}(\theta,\phi)A_{\rm
d}(\theta,\phi)A_{\rm a}(\theta,\phi)\end{displaymath} (15)
$A_{\rm n}$ is the natural apodization function due to the way the images are obtained: for intensity this is the limb darkening ($I(\mu)$), and for velocity the projection factor ($\mu=\sin \theta \cos \phi$). $A_{\rm d}$ 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 ($\mu=\sin \theta \cos \phi$). $A_{\rm a}$ 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:
 \end{eqnarray} (16)
It shows that $\tens{\cal{C}}^{(l,l')}$ is in general not hermitian nor symmetrical. Nevertheless, when $W_{l,m}=\tilde{Y}_{l,m}$, 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:
 ...eta,\phi)A(\theta,\phi) \sin \theta {\rm d}\theta
 {\rm d}\phi}\end{displaymath} (17)
which is the "natural" normalization factor of the perturbation $\tilde{Y}_{l,m}$. Of course in this latter case, we have:
 \end{eqnarray} (18)
Unfortunately, the leakage matrix does not always have such a nice property, especially because $W_{l,m}\ne\tilde{Y}_{l,m}$. 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 $\tilde{Y}_{l,m}$ have the same symmetry properties as the spherical harmonics Yl,m (or if $W_{l,m}=\tilde{Y}_{l,m}=Y_{l,m}$), the leakage matrix is real as shown by Schou (1992). In addition the leakage elements of $\tens{\cal{C}}_{m,m'}^{(l,l')}$ 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 $B \neq 0$ 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:  
 ...(l',m')} }
 \end{displaymath} (19)
where the $\tilde{y}_{i}$ are given by:  
 A(\theta,\phi) \sin \theta {\rm d}\theta {\rm d}\phi
 \end{displaymath} (20)
where ${\cal{D}}_{i}$ 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 $\tilde{y}_{i}^{(l,m)}$.

3.3.2 p-mode covariance matrix

To compute the covariance of the complex vector $\vec{y}(\nu)$ as a real number we form the vector $\vec{z}_{\vec{y}}(\nu)$ defined as:
\vec{z}_{\vec{y}}^{\rm T}(\nu)=(\mbox{Re}(\vec{y}^{\rm
 T}),\mbox{Im}(\vec{y}^{\rm T})).
In absence of noise, the covariance matrix $\tens{M}(\nu)$ of the vector $\vec{z}_{\vec{y}}(\nu)$ can be generated using a complex notation:
\tens{\cal{M...{M}}_{i}(\nu) & \tens{\cal{M}}_{r}(\nu)\\ \end{array}\right).\end{displaymath} (21)
$\tens{M}^{(l,l')}$ is a super matrix where $\tens{\cal{M}}_{r}(\nu)$ and $\tens{\cal{M}}_{i}(\nu)$ are the real and imaginary parts of a complex matrix $ \tens{\cal{M}}^{(l,l')}$ which elements are given by:  
 \end{displaymath} (22)
where $f_{m''}^{l''}(\nu)$ is the variance of the l'',m'' mode which profile is given by Eq. (11), in which $\nu_{0}$ 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 $\vec{y}$, and the covariance between the real and imaginary part of $\vec{y}$. It is obvious from Eq. (22) that $ \tens{\cal{M}}^{(l,l')}$ 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 $\tens{M}^{(l,l')}(\nu)$ can be decomposed as follows:  
 ...{w}^{\rm T}(\nu) & \tens{v}^{\rm T}(\nu)\\ \end{array}\right)
 \end{displaymath} (23)
where $\rm T$ is the transpose of a matrix. The elements of $\tens{v}$ and $\tens{w}$ are given by:  
 \end{displaymath} (24)
 \end{displaymath} (25)
We will see later on that this decomposition is of prime importance for understanding the statistics of the observation.

3.3.3 Noise covariance matrix

Unfortunately, the observed vector $\vec{y}(\nu)$ 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:  
\tens{\cal{B...{B}}_{i}(\nu) & \tens{\cal{B}}_{r}(\nu)\\ \end{array}\right)\end{displaymath} (26)
$\tens{B}^{(l,l')}$ is a super matrix where $\tens{\cal{B}}_{r}(\nu)$ and $\tens{\cal{B}}_{i}(\nu)$ are the real and imaginary parts of the complex matrix $\tens{\cal{B}}^{(l,l')}$. The dimension of $\tens{\cal{B}}^{(l,l')}$ is $(2l+2l'+2) 
\times (2l+2l'+2)$. Its elements are given by:
 a(\theta,\phi) \sin \theta {\rm d}\theta {\rm d}\phi
 \end{eqnarray} (27)
 \end{eqnarray} (28)
where a is an apodization function which characterizes through $\sigma^{2}_{\odot}(\theta,\phi,\nu)$ how the noise varies over the solar image, assuming that the noise is uncorrelated between different points on the Sun; $A_{\rm a}, A_{\rm d}$ are defined in Eq. (15). When the instrumental noise is low, a is derived from the characteristics of the solar noise. The evaluation of $\tens{\cal{B}}^{(l,l')}$ is less straightforward than that of $\tens{\cal{C}}^{(l,l')}$ 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 $\tens{\cal{R}}$ as:  
 \end{displaymath} (29)
Here we can see the similarity between $\tens{\cal{R}}$ and $\tens{\cal{C}}$. 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 $\tens{\cal{R}}^{(l,l')}$ very close to the leakage matrix $\tens{\cal{C}}^{(l,l')}$. Although $\tens{\cal{R}}^{(l,l')}$ 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 $\tens{\cal{B}}^{(l,l')}$ 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:  
 \end{displaymath} (30)
where bi is the variance of the noise of pixel i. Equation (30) is the more general form used for the LOI.

3.3.4 Probability density of the observation and likelihood

The statistical distribution of the Fourier spectra or of the vector $\vec{z}_{\vec{y}}$ is a multi-normal distribution. The probability density is given by:  
 ...{\vec {y}}(\nu)}}{(2\pi)^{d/2}\sqrt{\vert\tens{V}(\nu)\vert}}
 \end{displaymath} (31)
where d is the number of elements of $\vec{z}_{\vec{y}}$, $\tens{V}$ is a short notation for the following matrix: $\tens{V}^{(l,l')}(\nu)=
\tens{M}^{(l,l')}(\nu)+\tens{B}^{(l,l')}(\nu)$; 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 $\tens{V}^{(l,l')}(\nu)$ can also be built from sub-matrices as: $\tens{\cal{V}}^{(l,l')}=\tens{\cal{M}}^{(l,l')}+\tens{\cal{B}}^{(l,l')}$; as a result $\tens{\cal{V}}^{(l,l')}$ is also hermitian. Equation (31) is the most general formulation for any multi-normal distribution with a given covariance matrix $\tens{V}$ (Kendall & Stuart 1967).

Using Eq. (31), we can write the likelihood L of an observation of $\vec{z}_{\vec{y}}(\nu_{i})$ at N different frequencies $\nu_{i}$ as given by:  
 \end{displaymath} (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 $\vec{y}$, it is always possible in the absence of noise to recover the vector $\vec{x}$. 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):  
 \end{displaymath} (33)
where $\tens{\cal{C}}=\tens{\cal{C}}^{(l,l')}$. Then we can write a similar equation for $\vec{z}_{\vec{y}}$ and $\vec{z}_{\vec{\tilde{x}}}$ as:  
 \end{displaymath} (34)
where $\tens{C}$ is defined as:  
\tens{\cal{C...{C}}_{i}(\nu) & \tens{\cal{C}}_{r}(\nu)\\ \end{array}\right)\end{displaymath} (35)
$\tens{C}^{(l,l')}$ is a super matrix where $\tens{\cal{C}}_{r}$ and $\tens{\cal{C}}_{i}$ are the real and imaginary parts of the complex matrix $\tens{\cal{C}}^{(l,l')}$. Using Eq. (34) to replace $\vec{z}_{\vec{y}}$ by $\vec{z}_{\vec{\tilde{x}}}$ in Eq. (32) we can rewrite this latter as:  
 \end{displaymath} (36)
with $\tens{V'}$ given by:  
\tens{V'}=\tens{C}^{-1} \tens{V} {\tens{C}^{\rm T}}^{-1}=
\tens{B}^{(l,l')} {\tens{C}^{\rm T}}^{-1}.\end{displaymath} (37)
We recognize in Eq. (36) the probability density of the vector $\vec{z}_{\vec{\tilde{x}}}(\nu)$ to a constant (i.e. $\vert\tens{C}\vert^{-N}$). 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 $\tens{V'}$ of $\vec{z}_{\vec{\tilde{x}}}$ as written in Eq. (37) (Davenport & Root 1958). It can be easily shown using Eqs. (23) and (37) that the matrix $\tens{D}(\nu)=\tens{C}^{-1}
\tens{M}^{(l,l')} {\tens{C}^{\rm T}}^{-1}$ is diagonal and its element are given by:  
 \end{displaymath} (38)
where l''=l or l' and $m''=-l'',\ldots,l''$. 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 $\vec{\tilde{x}}$ after the transformation of Eq. (33). It means that $\vec{\tilde{x}}$ has no correlation due to the p modes as we could expect from Eq. (33): the leakage matrix of $\vec{\tilde{x}}$ 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 $\tens{C}$ by a proper matrix associated with $\tens{B}^{(l,l')}$. The derivation of this matrix is given in Appendix A.

3.3.5 The use of the likelihood in practice

When a single degree is observed, it is quite simple to maximize the likelihood of Eq. (32) using $\vec{y}$, or using $\vec{\tilde{x}}$ 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 $(m,\nu)$ 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 $\mu$Hz. In the $(m,\nu)$ 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 $\tens{B}^{(l,l')}$ in Eq. (37)), the elements of the vector $\vec{z}_{\vec{\tilde{x}}}(\nu)$ are all independent of each other, leading to a statistical distribution which is a product of $\chi^{2}$ 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 $\vec{z}_{\vec{y}}(\nu)$ as in Hill et al. (1996).

next previous
Up: The art of fitting

Copyright The European Southern Observatory (ESO)