MODELING SOLAR OSCILLATION POWER SPECTRA. II. PARAMETRIC MODEL OF SPECTRAL LINES OBSERVED IN DOPPLER-VELOCITY MEASUREMENTS

Helioseismic measurements of the deep solar interior are based on the inversions of the solar oscillation frequencies and frequency splittings inferred from observations of the velocity or intensity field on the visible part of the solar surface. With the exceptional quality and duration of the data now available from dedicated ground-based (GONG, BiSON) and space projects (SOHO Michelson Doppler Imager (MDI), GOLF, SDO Helioseismic and Magnetic Imager (HMI)), it has become obvious that progress in improving the accuracy and resolution of helioseismic measurements is limited not by the solar noise and finite mode lifetimes, but by the current procedures for measuring the p-mode frequencies and frequency splittings from the observed power spectra. The most damaging aspect of current procedures is the systematic errors they introduce; these are clearly seen both when comparing the results obtained with different data analysis pipelines (Schou et al. 2002; Basu et al. 2003) and in the helioseismic inversions where they manifest themselves as internal inconsistencies in the input data sets (Vorontsov 2002; Vorontsov et al. 2002). Systematic errors cannot be reduced by extending the duration of the observations and thus mask any gains in signal-to-noise ratio that come from the longer data set.

One source of the problems that limits the accuracy of the frequency measurements is uncertainties in the spatial response function (also known as the leakage matrix) for the observations. The spatial response function depends on the area of the Sun observed and on the details of the instrumentation. The latter can change with time and if the instrument is not accessible for comprehensive regular calibration (one issue with the Solar and Heliospheric Observatory (SOHO) MDI instrument is a gradual drift of the focal distance; see, e.g., Korzennik et al. 2004), this can be particularly damaging. Fortunately, this problem can be mitigated by using an adaptive response function (Vorontsov & Jefferies 2005; Paper I).

Another source of the problems that we address here is related to the modeling of the intrinsic line profiles. Current (published) frequency measurements are based on fitting individual lines with Lorentzian profiles. This is motivated by the analogy of solar oscillations with a simple damped single harmonic oscillator. However, many of the spectral peaks in the solar oscillation power spectrum are asymmetric. Although this asymmetry has typically been accommodated by adding an extra asymmetry parameter (e.g., Korzennik 2005) to the Lorentzian model, there are downsides to this approach: "asymmetric Lorentzian" line profiles are only applicable in the immediate vicinity of a sharp resonance, and the addition of the extra parameter can reduce the stability of the frequency measurements.

Ideally, theoretical simulation of observational power spectra would incorporate physical modeling of wave excitation and damping. Unfortunately, this approach is not feasible due to our poor knowledge of the turbulent convection and its interaction with waves in the outer solar layers. This lack of knowledge, however, does not mean that there is not a way to provide a model suitable for helioseismic inversions that is based on some physically motivated parameterizations. There is an analogy here with a similar difficulty that arises when the p-mode frequencies are used to study the solar internal structure. In this case, we have poor knowledge of the structure of the near-surface layers and of the physics of the wave behavior there. However, because the wave propagation in these layers is very nearly vertical, it is reasonable to expect that the near-surface effects change slowly with frequency and depend on frequency only (at least when the degree l is not too high). This a priori knowledge allows us to separate the uncertain effects into a frequency-dependent "surface phase shift," which in itself can be used later as a separate diagnostic of the near-surface layers. Without implementing this knowledge, which allows us to split the inverse problem, we would not be able to infer anything about the solar interior from the solar oscillation frequencies (more discussion on this topic can be found in Gough & Vorontsov 1995).

In developing our spectral model, we will follow the very same logic and will split the problem in exactly the same manner. We will make no attempt to consider the physics behind the mode excitation and damping. Instead, we will only assume that the two mechanisms operate close to the solar surface (an assumption that is supported by the observations). For a spherically symmetric Sun, both the excitation strength and damping should change slowly with frequency (slow on a scale of frequency separation between modes of consecutive radial order) and should depend on frequency only when the degree l is not too high. At a higher degree, just like in the "surface phase shift," the dependence on l at constant frequency is expected to emerge in proportion to l(l + 1), which governs the inclination of the acoustic ray paths to vertical in the leading order.

This work on spectral modeling is a part of a wider project targeted at streamlining the helioseismic inversion process (by analyzing the observational power spectra directly, without frequency measurements; Jefferies & Vorontsov 2004). However, it can also serve the more limited goal of improving the accuracy of frequency measurements. By allowing enough flexibility in the spectral model, we will substantially reduce systematic errors. In addition, by reducing the total number of fitting parameters (when this is allowed by the data), we will improve the stability of the measurements and reduce random errors.

Obviously, eliminating systematic errors without introducing new ones is only possible if the model is physically relevant and the underlying basic assumptions are sound. There are three tests that can be used to verify this. The first is how well the model fits the observations. The second is whether the fitted parameters meet the original expectations for their behavior (e.g., whether or not the excitation amplitude falls on a single function of frequency over a wide range in degree). The third and final test is helioseismic inversion, which is sensitive to the internal consistency of the input data set (the measured frequencies).

We will develop our model in the temporal frequency domain. We assume the Sun to be spherically symmetric (deviations from spherical symmetry are usually treated as small perturbations) and expand all the variables in our model in terms of spherical harmonics. We also assume that the wave excitation comes from uncorrelated excitation events that are localized somewhere near the solar surface and are random in both time and in the angular coordinates over the solar surface: this makes their contributions to the observational signal uncorrelated, a property that comes from the orthogonality of spherical harmonics. This, in turn, allows us to reduce the excitation problem to one of considering the response of the dynamical system to a single excitation event.

We start by considering a simplified ideal measurement in which the entire solar surface is observed so that a particular spherical-harmonic component in the signal can be perfectly isolated. We also consider the signal, defined over the solar surface, as a scalar field (e.g., an intensity perturbation or radial velocity). Complications arising from the line-of-sight projection of the vector velocity field that occurs in realistic measurements, which are limited by the visible part of the solar surface, will be addressed at the end of this section.

Before we present a mathematical derivation, it is convenient to first discuss the functional form of our spectral model, as well as its basic underlying assumptions.

The general expression for the amplitude of a signal observed in a particular spherical harmonic is

$$\begin{equation} {\cal A}_{\rm obs}= {{\cal A}_1e^{2i\theta }+{\cal A}_2\over 1-Re^{2i\theta }}+{\cal B}_{\rm c}. \end{equation} \tag{ 1 }$$

All the parameters in the right-hand side of this expression are slowly varying functions of frequency. θ is the phase integral of the trapped acoustic wave, a monotonically increasing function of frequency that takes values π(n + 1) at frequencies of acoustic resonances (n is the radial order). Energy dissipation is assumed to be localized in the sub surface layers, and the effects of energy losses are parameterized by the surface acoustic reflectivity of absolute value R. ${\cal A}_1$ and ${\cal A}_2$ are two complex excitation amplitudes. The complex variable ${\cal B}_c$ designates the "coherent component" of the solar noise (or the "direct visibility" of the excitation source), which is thought to be responsible for the different sign of line asymmetries in intensity and velocity measurements (Roxburgh & Vorontsov 1997; Nigam et al. 1998). A random (uncorrelated) component of the solar noise will be added to the analysis later (Equation (19)).

We provide two derivations of Equation (1). The first derivation is based on the analysis of the linear adiabatic oscillation equations. We note that the excitation domain is not limited to the deeper layers, where the adiabatic approximation to the oscillation equations is thought to be valid, and can well extend to the near-surface layers where the physics of the oscillations is poorly understood. However, the general functional form of Equation (1) does not change when the excitation operates partly from the non adiabatic domain (as we will see from the second derivation, which does not rely on the adiabatic approximation).

Separation of the variables in the oscillation equations is performed in a standard way with the displacement field u and the Eulerian pressure perturbation p' written as

$$\begin{equation} {\bf u}({\bf r},t)=[ U(r)\hat{\bf r}Y_{\ell m}(\theta,\phi) +V(r)\nabla _1Y_{\ell m}(\theta,\phi)]\,e^{-i\omega t}, \end{equation} \tag{ 2 }$$

$$\begin{equation} p^{\prime }({\bf r},t)=p_1(r)Y_{\ell m}(\theta,\phi)\,e^{-i\omega t}. \end{equation} \tag{ 3 }$$

Here $\nabla _1=\hat{\bf \theta }(\partial /\partial \theta) +(1/\sin \theta)\hat{\bf \phi }(\partial /\partial \phi)$ is the angular part of the gradient operator, and the hat symbol is used to designate unit vectors. The oscillation equations are implemented in the low-density envelope, where the Cowling approximation (which neglects perturbations of the gravitational field) is locally applicable. In this approximation the oscillation equations can be written as a single second-order equation:

$$\begin{equation} {d^2\psi \over d\tau ^2}+ [\omega ^2-V(\tau)]\psi =0, \end{equation} \tag{ 4 }$$

with

$$\begin{equation} \psi =\rho _0^{-1/2}r \left({1\over c^2}-{\tilde{w}^2\over r^2}\right)^{1/4} \left(1-{N^2\over \omega ^2}\right)^{-1/2}p_1 \end{equation} \tag{ 5 }$$

and the independent variable τ, defined as

$$\begin{equation} {d\tau \over dr}={1\over c} \left(1-\tilde{w}^2{c^2\over r^2}\right)^{1/2}, \end{equation} \tag{ 6 }$$

where $\tilde{w}^2=l(l+1)/\omega ^2$. The acoustic potential V is given by

$$\begin{eqnarray} V&=&N^2+\left[{r^2h\over \omega ^2-N^2}\left({1\over c^2}-{\tilde{w}^2\over r^2}\right)^{1/2}\right]^{-1/2}\nonumber\\ &&\times {d^2\over d\tau ^2} \left[{r^2h\over \omega ^2-N^2}\left({1\over c^2}-{\tilde{w}^2\over r^2}\right)^{1/2}\right]^{1/2}, \end{eqnarray} \tag{ 7 }$$

where N is the Brunt–Väisäla frequency and

$$\begin{equation} h=\exp \int \nolimits _0^r\left({N^2\over g_0}-{g_0\over c^2}\right)dr. \end{equation} \tag{ 8 }$$

It is straightforward to show that the time-averaged energy flux, carried by an arbitrary local solution to the adiabatic oscillation equations in the upward direction, is

$$\begin{equation} {i\over 4}\omega r^2(U^*p_1-Up_1^*)={i\over 4\omega }W(\psi,\psi ^*), \end{equation} \tag{ 9 }$$

where W(ψ₁, ψ₂) is the Wronskian of the two solutions ψ₁ and ψ₂,

$$\begin{equation} W(\psi _1,\psi _2)=\psi _1{d\psi _2\over d\tau }-\psi _2{d\psi _1\over d\tau }, \end{equation} \tag{ 10 }$$

and the star symbol designates complex conjugate.

We specify the external solution, which satisfies outer boundary conditions, by a sum of a wave incident on the surface and a reflected wave:

$$\begin{equation} \psi _e={\cal A}_e(\psi _++R\psi _+^*), \end{equation} \tag{ 11 }$$

where ψ₊ is an upward wave, $\psi _+^*$ is a downward wave, R is a real-valued reflectivity coefficient, and ${\cal A}_e$ is a complex-valued function of frequency. The internal solution, which satisfies the inner boundary conditions, is specified as a standing wave:

$$\begin{equation} \psi _i={\cal A}_i(e^{i\theta }\psi _++e^{-i\theta }\psi _+^*), \end{equation} \tag{ 12 }$$

where θ is the phase integral across the acoustic cavity; when the reflectivity R is set to 1, the two solutions match each other when θ = π(n + 1) (eigensolutions of the adiabatic problem).

To address the wave excitation, we add Dirac's δ-function δ(τ − τ_s) to the right-hand side of Equation (4). A particular solution of the resulting inhomogeneous equation is the Green's function

$$\begin{eqnarray} &&G(\tau,\tau _s)\nonumber\\ &&=\left\lbrace\begin{array}{@{}ll@{}} -W^{-1}(\psi _+,\psi _+^*)[\psi _+(\tau _s)\psi _+^*(\tau)-\psi _+^*(\tau _s)\psi _+(\tau)], & \tau \le \tau _s \\ 0, & \tau \ge \tau _s.\end{array} \right.\nonumber\\ && \end{eqnarray} \tag{ 13 }$$

A general solution that satisfies the outer boundary conditions is G(τ, τ_s) + ψ_e, where ψ_e is given by Equation (11). By matching this solution with the internal solution ψ_i (Equation (12)), we obtain the amplitude of the external solution

$$\begin{equation} {\cal A}_e=-W^{-1}(\psi _+,\psi _+^*) {\psi _+(\tau _s)e^{2i\theta }+\psi _+^*(\tau _s)\over 1-Re^{2i\theta }}. \end{equation} \tag{ 14 }$$

We see that the right-hand side of Equation (14) fits the functional form of the right-hand side of Equation (1). Since an arbitrary excitation source can be considered as a linear superposition of "elementary" excitations described by δ-functions, it will bring an observational signal of the functional form specified by Equation (1).

An alternative derivation comes from considering wave propagation and interference in the acoustic cavity in analogy with light propagation in a Fabry–Perot interferometer (e.g., Duvall et al. 1993). Consider two waves emitted by an excitation source that is localized near the surface: an upward wave and a downward wave. The upward wave is reflected almost immediately from the surface and then undergoes multiple further reflections from the surface as it travels up and down across the acoustic cavity. The signal observed at the surface is evaluated as a sum of a geometric series with the common ratio Re^2iθ, where 2θ is the phase lag between two consecutive reflections from the surface. This signal gives the term ${\cal A}_2/(1-Re^{2i\theta })$ in Equation (1). The other term in Equation (1), ${\cal A}_1e^{2i\theta }/(1-Re^{2i\theta })$, comes from the wave that is emitted downward and accumulates an additional phase lag of 2θ before reaching the surface.

To obtain a convenient expression for the signal power $|{\cal A}_{\rm obs}|^2$, we introduce the variable transformation

$$\begin{equation} \tan \theta ={1-R\over 1+R}\tan \varphi, \end{equation} \tag{ 15 }$$

which gives

$$\begin{equation} {1\over 1-Re^{2i\theta }}={1+Re^{2i\varphi }\over 1-R^2},\quad {e^{2i\theta }\over 1-Re^{2i\theta }}={R+e^{2i\varphi }\over 1-R^2}. \end{equation} \tag{ 16 }$$

Equation (1) can thus be written as

$$\begin{equation} {\cal A}_{\rm obs}={b+ae^{2i\varphi }\over 1-R^2} \end{equation} \tag{ 17 }$$

with

$$\begin{equation} b=R{\cal A}_1+{\cal A}_2+(1-R^2){\cal B}_c,\quad a={\cal A}_1+R{\cal A}_2, \end{equation} \tag{ 18 }$$

and we have

$$\begin{equation} \left|{\cal A}_{\rm obs}\right|^2 =\left[{A\cos \left(\varphi -S\right)\over 1-R^2}\right]^2 +B^2, \end{equation} \tag{ 19 }$$

where

$$\begin{equation} A^2=4\left|ab\right|, \end{equation} \tag{ 20 }$$

$$\begin{equation} B^2=\left(|a|-|b|\over 1-R^2\right)^2+B_u^2, \end{equation} \tag{ 21 }$$

$$\begin{equation} 2S=\arg b -\arg a, \end{equation} \tag{ 22 }$$

and we have added an uncorrelated component of the solar background $B_u^2$ to the right-hand side of Equation (19).

We now designate L_U and L_V as the two components of the response function (leakage matrix) that specify the sensitivity of the spatial filter of the Doppler-velocity measurements to the vertical and horizontal velocity components, respectively. The observational power spectrum is then given by

$$\begin{equation} P=\left|L_U+hL_V\right|^2 \left\lbrace \left[{A\cos \left(\varphi -S\right)\over 1-R^2}\right]^2 +B^2\right\rbrace, \end{equation} \tag{ 23 }$$

where h is the ratio of horizontal and vertical components of the displacement field on the solar surface (V/U, see Equation (2)). A simple theoretical prediction for this quantity is

$$\begin{equation} h={GM\over R^3\omega ^2}. \end{equation} \tag{ 24 }$$

This value comes from the horizontal component of the momentum equation in the Cowling approximation, when the Lagrangian pressure perturbation is zero; see, e.g., Cox (1980). It is thus applicable to any divergence-free flow at any frequency, not only at the frequencies of acoustic resonances. Our practical experience with fitting solar velocity power spectra provides no evidence that this ratio shall be modified, a conclusion that is in agreement with Rhodes et al. (2001). We note, however, that a small variation of h cannot be distinguished from a small distortion of the leakage matrix caused by an error in the plate scale (Paper I). Discussions of other efforts of measuring h can be found in Rhodes et al. (2001) and Korzennik et al. (2004).

Equation (23) represents the spectral model that we use for matching the observed power spectra by way of an iterative maximum-likelihood technique. At each degree l, the parameters A, R, S, and B are expected to change slowly with frequency (we discard their possible variation with azimuthal order m in the results presented below). A is an effective excitation amplitude and R is the acoustic reflectivity of the solar surface. These two parameters replace the mode amplitude and line width that are used in traditional spectral models for the oscillation modes. Since the excitation and damping mechanisms are thought to operate in the near-surface layers, where wave propagation is nearly vertical, we expect A and R to depend on frequency only (not on the degree), at least when the degree l is not too high. The parameter S specifies line asymmetry: its relation with "asymmetric Lorentzian" line profiles is discussed in Appendix A. B is a composite background with contributions from both correlated and uncorrelated components of the solar noise. The functional form of Equation (23) assumes that horizontal and vertical components of the background are in the same ratio as those of the resonant signal. This is hardly the case for solar granulation noise. The inaccuracy of this assumption leads to a distortion of the measured values of B². We will get back to this point later when discussing our numerical results.

For a spherically symmetric configuration, the positions of the spectral resonances (as well as the line profiles at given R and S) are governed by the phase integral θ(ω), which has to be specified in the spectral model for each value of degree l. Ideally, θ(ω) is computed from an equilibrium solar model that fits the seismic data. As such a model is not available (it is itself a result of helioseismic structural inversion of the measured oscillation frequencies), we use an iterative approach where θ(ω) is parameterized by the positions of the spectral resonances ω_n, where θ(ω_n) = π(n + 1). A continuous variation of θ(ω) between the resonances is represented by a spline with prescribed gradients dθ/dω at the resonant frequencies ω = ω_n. The values of these gradients are evaluated using solutions to the adiabatic oscillation equations computed for a proxy solar model. Using the unrestricted variational principle described in Vorontsov et al. (2013), it is straightforward to show that

$$\begin{equation} {1\over \omega ^2}W(\psi _i,\psi _e)=\delta (\omega ^2) \int \nolimits _0^R\rho _or^2[U^2+l(l+1)V^2]dr. \end{equation} \tag{ 25 }$$

In this equation, ψ_i is the solution to the adiabatic oscillation equations that satisfies the inner boundary conditions, and ψ_e is the solution that satisfies the (conservative) outer boundary conditions and $\delta (\omega ^2)=\omega ^2-\omega _n^2$. Using this variational principle, we have

$$\begin{equation} {d\theta \over d\omega }|_{\omega =\omega _n}= {2\omega ^2\over \psi ^2+{1\over \omega ^2}\left({d\psi \over d\tau }\right)^2} \int \nolimits _0^R\rho _or^2[U^2+l(l+1)V^2]dr, \end{equation} \tag{ 26 }$$

with ψ and dψ/dτ evaluated at the depth where V(τ) ≪ ω² (see Equation (4)). We note that due to the high-frequency asymptotic properties of solar p modes, the functions θ(ω) are close to linear functions of frequency ω. A small inaccuracy in the gradients dθ/dω leads only to a small distortion of the measured values of the acoustic reflectivity R (see Equation (15): the line width in the power spectrum depends on both R and dθ/dω).

For the simple scenario where A, R, and B do not depend on frequency, θ changes linearly with ω, and S = 0; the power spectrum predicted by Equation (19) can be represented by a sum of Lorentzian profiles (see Appendix A).

Using Equations (A8) and (A9) we find that the line width at half power, measured in units of angular frequency ω, is

$$\begin{equation} \Lambda =-2\gamma =-{\ln R\over d\theta /d\omega }. \end{equation} \tag{ 27 }$$

This relation shows that even when R is degree-independent, the line width depends explicitly on the degree l, being inversely proportional to the "mode mass" (the integral in Equation (26) is proportional to mode energy).

The oscillation frequencies and frequency-splitting coefficients are measured using an iterative procedure. Our preliminary analysis (with results presented in the next section) was targeted at measuring the parameters A, R, S, and B for individual modes to check if they reveal the expected behavior (a variation with frequency only, at least when the degree l is not too high). For this reason, the adjustments of the frequencies and splittings were alternated with improvements of A, R, S, and B as functions of frequency at given l. A dependence of A, R, S, and B on frequency inside the fitting domains was allowed for by adding a uniform shift to a single (piecewise linear) function of frequency. This function describes the average variation of the corresponding parameter with frequency for lower-degree modes (updated in the outer iteration after the inner iterations converge). In the results presented below, the size of the fitting domain in frequency was taken to be 10 line widths for measuring A and R, an average of the frequency spacings between modes of consecutive order when measuring S and B, and about 6 μHz when measuring the frequencies and frequency splittings. Separate tests where performed to verify that the fitted parameters do not suffer from systematic variations when the sizes of the fitting intervals change.

The leakage matrix was computed as described in Paper I, but with two extensions to the analysis targeted at further improvement of its accuracy. One extension is a full description of the effects due to a non zero solar B angle (in Paper I these effects were only allowed as a perturbation). This extension to the leakage-matrix computation is described in Appendix B. The other extension is a more accurate treatment of the mode coupling by differential rotation, which is now performed as described in Vorontsov (2007). Spatial leaks in the range of ±30 in l and m were included in the spectral modeling. The large number of leaks had to be taken into account for addressing B at frequencies higher than about 2 mHz, where the solar background appears to be buried below the resonant signals of spatial leaks.

We note that the spectral model described by Equation (23) is applicable not only to acoustic resonances, but also to the high-frequency interference pikes, or pseudomodes (R → 0) observed at frequencies higher than the "acoustic cutoff" frequency (about 5 mHz). In addition, adaptation of the spectral model to observations in intensity only requires a modification of the response function (L_U + hL_V in Equation (23)).

A representative example of the one-year SOHO MDI power spectrum, obtained from the first year of observations (at low solar activity), and its model is shown in Figure 1. Here the spectra corresponding to different values of azimuthal order m was shifted in frequency and then averaged for better visual comparison. The observed complexity of the asymmetric spectral lines is due to the contribution of multiple spectral leaks. We note the noise level is difficult to detect as it is below the amplitudes of resonant signals.

Zoom In Zoom Out Reset image size

The parameters of the spectral model are illustrated in Figures 2–5. The results are shown for p modes of radial order n from 1 to 10 and for the f modes of degree l ⩽ 200. The maximum-likelihood solution was obtained by an iterative improvement of both the spectral parameters (A, R, S, and B) and the resonant frequencies and frequency splittings. The resulting agreement between the model and the data is satisfactory, as can be judged from a proper merit function (Anderson et al. 1990) and from the visual inspection of the m-averaged power spectra. However, a small systematic inaccuracy in the predicted amplitudes of the spatial leaks remains. The origin of this mismatch is not yet understood; it may be related, in part, to the asymmetric distortion of the point-spread function of the MDI instrument.

Zoom In Zoom Out Reset image size

The spectral parameters of individual modes do not depend on the degree and collapse to slowly varying functions of frequency only, when the degree l is not too high (less than about 100). The composite background B² (Figure 5) may look to be an exception to this. However, an accurate measurement of the background at frequencies higher than about 2 mHz is a difficult task as the background level appears to be significantly smaller than the resonant signals coming from the spatial leaks. As a result, the measurement of B² can be distorted by small inaccuracies in the leakage matrix. The fitted background is also significantly higher than the average at the lowest values of degree l (e.g., p₁₀ mode of l = 2 in Figure 5). A likely explanation of this excess is instrument-related noise due to variations in the exposure times for the series of observations acquired at different wavelengths by MDI and used to compute the solar velocity signal (J. Schou 2013, private communication). This noise can probably be modeled by adding an l = 0 component to the (otherwise degree-independent) solar background B²(ω).

Interestingly, the analysis of solar f-modes reveals the same values of the excitation amplitudes A and "acoustic reflectivity" R as solar p-modes of similar frequencies (Figures 2 and 3). This indicates that the excitation and damping mechanisms do not distinguish between p and f modes, despite the difference in their physical nature (the f modes are incompressible waves). The composite background of f modes appears to be smaller than that of p modes (Figure 5). Since the fitted background B² at low frequencies is dominated by the granulation noise, one possible explanation is that its contribution to the observational power should be modeled with a smaller (or zero) value of h, which specifies the ratio of horizontal and vertical velocities. When processing the data we used a simple theoretical value (Equation (24)) that corresponds to incompressible adiabatic motion. This approximation is hardly relevant to granulation noise.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 6 shows the rotational splitting coefficients resulting from the measurement process described above, in comparison with published splitting coefficients (Schou 1999). The horn-like structures in the published results, which signify systematic errors, are eliminated. As indicated by a detailed analysis (Vorontsov et al. 2009), the dominant part of the systematic errors came from discarding the effects of mode coupling by differential rotation in the original version of the SOHO MDI data analysis pipeline (which was later improved to include these effects, among others; see Larson & Schou 2008, who also observed a reduction in the horn-like structures).

Zoom In Zoom Out Reset image size

The differences between the centroid frequencies and their published values are shown in Figure 7. The major part of the discrepancies is due to the line asymmetries that were not accounted for in the original version of the MDI data analysis pipeline. Smaller-scale features are due apparently to the combination of the mode-coupling effects with the plate-scale error. The centroid frequencies measured in this work have been used in a recent study targeted at the seismic diagnostics of the equation of state (Vorontsov et al. 2013). Interestingly, it was found that these frequencies allow us to achieve a significantly better agreement with solar models. This is quite an unusual finding: better accuracy of observational data brings better agreement with theoretical models.

Zoom In Zoom Out Reset image size

When tested with high-quality one-year SOHO MDI Doppler-velocity measurements, the spectral model suggested in this article allows for an accurate description of the observational power spectra over a wide (continuous) range of frequency and in the degree range up to l of about 200. The oscillation frequencies and frequency splittings measured by the new technique appear to have smaller systematic errors when compared with published results.

A distinctive feature of the technique is its potential ability to reduce the random errors in the frequency and frequency-splitting measurements of the p modes that penetrate into the deep solar interior. This benefit comes from the global approximations of the spectral parameters of the model (excitation amplitude A, acoustic reflectivity R, line asymmetry S, and composite background B). At both low- and intermediate-degree l these parameters can be represented by slowly varying functions of frequency only (a possible exception is the solar background B, which requires further study). When inferred from the large volume of high-quality, intermediate-degree data, these global approximations can be used in measurements at lower degree l. The subsequent reduction in the number of free parameters in the mode-fitting algorithm is expected to bring significant improvement to the accuracy and precision of the frequency and frequency-splitting measurements in this important region of the solar oscillation power spectrum (which contains information on the solar core).

A problem that is not yet resolved is that the performance of the technique degrades when analyzing the MDI data of a degree higher than about 200. This is apparently due to an inaccuracy in the predicted amplitudes of the spatial leaks. This leads to poor stability (and unavoidable systematic errors) of the frequencies and frequency splittings at a higher degree where the spatial leaks start to blend into a single ridge. The origin of this problem may be related, at least in part, with the asymmetric distortion of the point-spread function of the MDI instrument (see, e.g., Korzennik et al. 2004, 2013), which is not accounted for in our model.

Our future work is targeted at processing the entire volume of the SOHO MDI data accumulated over the 15 years of the MDI operational lifetime, as well as at the implementation of our technique to the analysis of Solar Dynamics Observatory (SDO) HMI data. On the theoretical side, the ultimate goal of the project is a streamlined helioseismic inversion. Here frequency measurements will be eliminated from the analysis, and the parameters of the rotating and aspherical solar model will be matched directly with p-mode power spectra. This approach will provide the benefits of streamlined regularization by eliminating the problems associated with error correlation and possible mode misidentification in the frequency measurements.

We thank Jesper Schou and Tim Larson for providing Doppler-velocity power spectra from their analysis of SOHO MDI data and for illuminating discussions of possible instrumental and observational distortions in these measurements. We also thank an anonymous referee for suggestions that helped to improve the presentation. SOHO is a project of international cooperation between ESA and NASA. This work was supported in part by the UK STFC under grant PP/E001459/1.

Using the identity

$$\begin{eqnarray} e^{\pm i\theta }&=&R^{\mp 1/2}\left[\cos \left(\theta -{i\over 2}\ln R\right)\pm i\sin \left(\theta -{i\over 2}\ln R\right)\right], \nonumber\\ && \end{eqnarray} \tag{ A1 }$$

Equation (1), which describes the amplitude spectra, can be written as

$$\begin{eqnarray} {\cal A}_{\rm obs}&=&{i\over 2}\left({1\over R}{\cal A}_1+{\cal A}_2\right)\cot \left(\theta -{i\over 2}\ln R\right)\nonumber\\ &&-{1\over 2}\left({1\over R}{\cal A}_1-{\cal A}_2\right)+{\cal B}_c. \end{eqnarray} \tag{ A2 }$$

Introducing the independent variable

$$\begin{equation} z=\theta -{i\over 2}\ln R, \end{equation} \tag{ A3 }$$

which is now a complex quantity, and expanding $\cot z$ in partial fractions, we have

$$\begin{eqnarray} {\cal A}_{\rm obs}&=&{1\over 2}\left({1\over R}{\cal A}_1+{\cal A}_2\right) \left[{i\over z}+\sum \limits _{k=1}^\infty \left({i\over z-\pi k}+{i\over z+\pi k}\right)\right]\nonumber\\ &&-{1\over 2}\left({1\over R}{\cal A}_1-{\cal A}_2\right)+{\cal B}_c. \end{eqnarray} \tag{ A4 }$$

The right-hand side of this equation has simple poles at z = ±πk,   k = n + 1, and describes ${\cal A}_{\rm obs}$ as being produced by a superposition of damped harmonic oscillations of consecutive radial orders. Since the complex amplitudes ${\cal A}_1$ and ${\cal A}_2$ are common for all the modes, the oscillations are excited in a coherent way, and the resultant power spectrum cannot be represented, in general, by a simple sum of power spectra of individual modes (modes described by the same spherical harmonic cannot be considered as uncorrelated). It can also be seen from Equation (A4) that when ${\cal A}_1, {\cal A}_2$, and R are constants, θ changes linearly with ω, ${\cal B}=0$ and the contribution of neighboring resonances is discarded, and an individual resonant line in the power spectrum $|{\cal A}_{\rm obs}|^2$ will have a symmetric Lorentzian profile only if $\arg ({\cal A}_1/R+{\cal A}_2)=\arg ({\cal A}_1/R-{\cal A}_2)$, which is equivalent to $\arg ({\cal A}_1)=\arg ({\cal A}_2)$.

To address the profile of the resultant power spectrum, we perform a similar decomposition but for the observational power $|{\cal A}_{\rm obs}|^2$ (Equation (19)). Using Equation (15), we have

$$\begin{eqnarray} |{\cal A}_{\rm obs}|^2&=&A^2\left[{\cos 2S\over 8R^2}\, {1+{2R\over 1-R^2}\tan 2S\sin 2\theta \over {1+R^2\over 2R}-\cos 2\theta }\right.\nonumber\\ &&\left.-{1\over (1-R^2)^2}\left({1+R^2\over 4R}\cos 2S-{1\over 2}\right)\right]+B^2.\nonumber\\ && \end{eqnarray} \tag{ A5 }$$

Using (1 + R²)/(2R) = cosh ln R,   (1 − R²)/(2R) = sinh ln R, and the standard expression for wrapped Cauchy distribution

$$\begin{equation} \sum \limits _{k=-\infty }^\infty {\beta \over \pi [\beta ^2+(\alpha -\mu +2\pi k)^2]}= {1\over 2\pi }{\sinh \beta \over \cosh \beta -\cos (\alpha -\mu)}, \end{equation} \tag{ A6 }$$

we arrive to the desired expansion

$$\begin{eqnarray} |{\cal A}_{\rm obs}|^2&=&A^2\left[{-}{\ln R\cos 2S\over 8R(1-R^2)} \left(1+{2R\over 1-R^2}\tan 2S\sin 2\theta \right)\right.\nonumber\\ &&\left.\times \sum \limits _{k=-\infty }^\infty {1\over (\theta -\pi k)^2+{1\over 4}\ln ^2R}\right. \nonumber\\ &&\left.-{1\over (1-R^2)^2}\left({1+R^2\over 4R}\cos 2S-{1\over 2}\right)\right]+B^2. \end{eqnarray} \tag{ A7 }$$

Consider a simple scenario when A, R, S, and B do not depend on frequency ω, and θ is proportional to ω, with θ = πω/Δω, so that the resonant frequencies are equidistant. Discarding normalization and additive constants, the power spectrum is proportional to

$$\begin{equation} \left[1+{2R\over 1-R^2}\tan 2S\sin \left({2\pi \omega \over \Delta \omega }\right)\right] \sum \limits _{k=-\infty }^\infty {1\over (\omega -k\Delta \omega)^2+\gamma ^2} \end{equation} \tag{ A8 }$$

with

$$\begin{equation} \gamma ^2=\left({\Delta \omega \over 2\pi }\ln R\right)^2. \end{equation} \tag{ A9 }$$

When S = 0, the power spectrum is a simple sum of Lorentzian profiles, an unexpected conclusion. When the asymmetry parameter differs from zero, the spectrum is modulated with a harmonic function of frequency with period equal to the frequency spacing Δω.

In current work in solar seismology, a popular functional form for the asymmetric line profiles of solar p modes is suggested by Nigam & Kosovichev (1998):

$$\begin{equation} {(1+Bx)^2+B^2\over 1+x^2}, \end{equation} \tag{ A10 }$$

where x = (ω − ω_n)/γ; this expression was obtained by using a power-series expansion at small ω − ω_n and hence is only applicable in the vicinity of the resonant frequency. One obvious problem with this profile it that when x → ±∞, the power density tends to a non-zero value of B², which corresponds to an infinite mode energy. Matching the line profiles predicted by expressions (A8) and (A10) in the vicinity of the resonance, we have S ≃ B. Korzennik et al. (2004) suggested modification to the Nigam & Kosovichev's profile of the form

$$\begin{equation} {1+\alpha \left(x-\alpha /2\right)\over 1+x^2}. \end{equation} \tag{ A11 }$$

This profile brings negative power density when x → ∞ or x → −∞, depending on the sign of the asymmetry parameter α. Matching our result with this profile in the vicinity of x = 0, we have S ≃ α/2.

Another approach to modeling observational signal over a continuous frequency range (all the way between the consecutive acoustic resonances) was suggested recently by Hindman (2011). The expression suggested for the observational amplitude differs markedly from our Equation (A4): instead of a single frequency-dependent function ${\cal A}_1/R+{\cal A}_2$ that defines the normalization of all the partial fractions i/(z − πk) (Equation (A4)), each partial fraction enters the summation with individual complex amplitude. Line asymmetry in the power spectrum is interpreted in Hindman (2011) as a result of mode interference (cross-terms between partial fractions describing individual damped harmonic oscillations). In our view, the origin of the line asymmetry is very different: when ${\cal B}_c=0$, the asymmetry comes from the emission, by the same excitation event, of both an upward- and downward propagating waves; the sign of the asymmetry and its strength are thus governed by the parity and depth of the excitation source (see, e.g., Roxburgh & Vorontsov 1995, 1997 and references therein). When observations are interpreted in terms of damped harmonic oscillations, this vision can be supported by a simple illustrative example. Consider an instantaneous excitation of a single mode of frequency ω₀ and damping rate γ, followed by another similar excitation after a short time interval Δt specified by the temporal behavior of the excitation event. The combined power spectrum is proportional to

$$\begin{equation} \left|{1\over \omega -\omega _0-i\gamma }+{e^{-i\omega \Delta t}\over \omega -\omega _0-i\gamma }\right|^2= {2\left(1+\cos \omega \Delta t\right)\over (\omega -\omega _0)^2+\gamma ^2}, \end{equation} \tag{ A12 }$$

and the spectral line is asymmetric when ω₀Δt ≠ πk, where k is an integer, i.e., when the two excitations are not exactly in phase (or in counter phase). The pair of consecutive excitations can be interpreted as being produced by the upward and the downward wave. It is seen from this example that line asymmetry is governed directly by the spatio-temporal properties of a single excitation event and does not require interference with modes of neighboring frequencies. Interpretation of line asymmetry it terms of mode interference also faces difficulty when the frequency separation Δω is much bigger than the line width, i.e., when the amplitude of a neighboring mode is negligibly small in the vicinity of a sharp resonance.

This appendix should be read in conjunction with Paper I; it contains the replacement to Section 2.2 of Paper I, to allow an explicit account of the leakage-matrix computation for non zero solar B angle. The rotation of the coordinate system is now performed with β = π/2 + B. Equations (25) and (26) of Paper I are replaced with more general transformation relations

$$\begin{equation} Y_{lm}(\theta ^{\prime },\phi ^{\prime })=\sum \limits _{m^{\prime }}d_{mm^{\prime }}^{(l)}(\beta)Y_{lm^{\prime }}(\theta,\phi), \end{equation} \tag{ B1 }$$

$$\begin{equation} Y_{lm}(\theta,\phi)=\sum \limits _{m^{\prime }}d_{m^{\prime }m}^{(l)}(\beta)Y_{lm^{\prime }}(\theta ^{\prime },\phi ^{\prime }), \end{equation} \tag{ B2 }$$

with transformation coefficients

$$\begin{eqnarray} d_{m^{\prime }m}^{(l)}(\beta)&=&\sum \limits _t{(-1)^t\sqrt{(l+m)!(l-m)!(l+m^{\prime })!(l-m^{\prime })!}\over (l-m^{\prime }-t)!(l+m-t)!t!(t+m^{\prime }-m)!}\nonumber\\ &&\times \left(\cos {\beta \over 2}\right)^{2l+m-m^{\prime }-2t}\left(\sin {\beta \over 2}\right)^{2t+m^{\prime }-m}, \end{eqnarray} \tag{ B3 }$$

where t is run over those integer values for which all the factorials are those of positive integers or zero (see Equations (A1) and (A2) of Paper I); the transformation coefficients satisfy symmetry relations

$$\begin{equation} d_{m^{\prime }m}^{(l)}(\beta)=(-1)^{m+m^{\prime }}d_{mm^{\prime }}^{(l)}(\beta)=d_{-m,-m^{\prime }}^{(l)}(\beta), \end{equation} \tag{ B4 }$$

$$\begin{equation} d_{m^{\prime }m}^{(l)}(\pi -\beta)=(-1)^{l+m^{\prime }}d_{m^{\prime },-m}^{(l)}(\beta). \end{equation} \tag{ B5 }$$

The recurrence relations, convenient for calculating the transformation coefficients, are available in Varshalovich et al. (1988). After proper sign corrections (due to a transposed definition of the expansion coefficients, see Equation (5.5.1) of Varshalovich et al. (1988), these relations are Equations (4.8.16) and (4.8.17) of Varshalovich et al. (1988)

$$\begin{eqnarray} &&\sin \beta \left[(l+m)(l-m+1)\right]^{1/2}d_{m^{\prime },m-1}^{(l)}(\beta)\nonumber\\ &&\quad- 2(m^{\prime }-m\cos \beta)d_{m^{\prime }m}^{(l)}(\beta)\nonumber\\ &&\quad+\sin \beta \left[(l+m+1)(l-m)\right]^{1/2}d_{m^{\prime },m+1}^{(l)}(\beta)=0, \end{eqnarray} \tag{ B6 }$$

$$\begin{eqnarray} &&\sin \beta [(l+m)(l-m+1)]^{1/2}d_{m-1,m^{\prime }}^{(l)}(\beta)\nonumber\\ &&\quad+ 2(m^{\prime }-m\cos \beta)d_{m,m^{\prime }}^{(l)}(\beta)\nonumber\\ &&\quad+\sin \beta \left[(l+m+1)(l-m)\right]^{1/2}d_{m+1,m^{\prime }}^{(l)}(\beta)=0. \end{eqnarray} \tag{ B7 }$$

Due to symmetry relations (Equation (B4)), only a quarter of the matrix of transformation coefficients at each l needs to be evaluated using the recurrence relations; when Equation (B7) is in use the recurrence starts at m' = l using the explicit expressions

$$\begin{equation} d_{lm}^{(l)}(\beta)=\left[{(2l)!\over (l+m)!(l-m)!}\right]^{1/2} \left(\cos {\beta \over 2}\right)^{l+m}\left(\sin {\beta \over 2}\right)^{l-m}. \end{equation} \tag{ B8 }$$

The static part of the leakage matrix is

$$\begin{equation} S_{l^{\prime }l}^{m^{\prime }m}(\beta)= \sum \limits _{\mu =-l}^l\sum \limits _{\mu ^{\prime }=-l^{\prime }}^{l^{\prime }} d_{\mu m}^{(l)}(\beta)d_{\mu ^{\prime } m^{\prime }}^{(l^{\prime })}(\beta)C_{l^{\prime }l}^{\mu ^{\prime }\mu }. \end{equation} \tag{ B9 }$$

This expression replaces Equation (35) of Paper I. Deviation of β from π/2 (non-zero B angle) brings "prohibited" leaks (with l + l' + m + m' odd). The static part of the leakage matrix obeys symmetry properties

$$\begin{equation} S_{l^{\prime }l}^{-m^{\prime },-m}(\beta)=(-1)^{m+m^{\prime }}S_{l^{\prime }l}^{m^{\prime }m}(\beta), \end{equation} \tag{ B10 }$$

$$\begin{equation} S_{l^{\prime }l}^{m^{\prime }m}(\pi -\beta)=(-1)^{l+l^{\prime }+m+m^{\prime }}S_{l^{\prime }l}^{m^{\prime }m}(\beta). \end{equation} \tag{ B11 }$$

It can be seen from the last relation that the "prohibited" leaks are zero when β = π/2. At small B = β − π/2, the amplitudes of "prohibited" leaks are proportional to B, and variations of the amplitudes of "unprohibited" leaks is quadratic in B. When analyzing the one-year power spectra, we thus calculate the leakage matrix with B = 511, which gives a proper annual average for B².

Using orthogonality and normalization properties of spherical harmonics, it can be seen that

$$\begin{equation} \sum \limits _{m=-l}^ld_{pm}^{(l)}(\beta)d_{qm}^{(l)}(\beta)=\delta _{pq}; \end{equation} \tag{ B12 }$$

this relation gives

$$\begin{equation} \sum \limits _{m=-l}^l\sum \limits _{m^{\prime }=-l^{\prime }}^{l^{\prime }}\left(S_{l^{\prime }l}^{m^{\prime }m}(\beta)\right)^2= \sum \limits _{m=-l}^l\sum \limits _{m^{\prime }=-l^{\prime }}^{l^{\prime }}\left(C_{l^{\prime }l}^{m^{\prime }m}\right)^2, \end{equation} \tag{ B13 }$$

which means that the total power of leaks with given l' to target l does not depend on the inclination angle β.