{'fraction_non_alphanumeric': 0.039813960947836895, 'fraction_numerical': 0.021330339601459445, 'mean_word_length': 4.522307096202812, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Measurement of low signal-to-noise-ratio solar p modes in spatially-resolved helioseismic data

We present an adaptation of the rotation-corrected, m-averaged spectrum technique designed to observe low signal-to-noise-ratio, low-frequency solar p modes. The frequency shift of each of the 2l+1 m spectra of a given (n,l) multiplet is chosen that maximizes the likelihood of the m-averaged spectrum. A high signal-to-noise ratio can result from combining individual low signal-to-noise-ratio, individual-m spectra, none of which would yield a strong enough peak to measure. We apply the technique to GONG and MDI data and show that it allows us to measure modes with lower frequencies than those obtained with classic peak-fitting analysis of the individual-m spectra. We measure their central frequencies, splittings, asymmetries, lifetimes, and amplitudes. The low-frequency, low- and intermediate-angular degrees rendered accessible by this new method correspond to modes that are sensitive to the deep solar interior down to the core and to the radiative interior. Moreover, the low-frequency modes have deeper upper turning points, and are thus less sensitive to the turbulence and magnetic fields of the outer layers, as well as uncertainties in the nature of the external boundary condition. As a result of their longer lifetimes (narrower linewidths) at the same signal-to-noise ratio the determination of the frequencies of lower-frequency modes is more accurate, and the resulting inversions should be more precise.


INTRODUCTION
Our knowledge of the structure and dynamics of the solar interior has been considerably improved by the use of measurements of the properties of the normal modes of oscillation of the Sun. However, the Sun's interior is far from being fully understood, and better measurements of the mode parameters will also help to better understand the mode excitation and damping mechanisms as well as the physical properties of the outer layers by better constraining the turbulence models. A large number of predicted acoustic oscillation modes, defined by their radial orders (n) and their angular degrees (l), are not yet observed in the lowfrequency range (i.e., approximately below 1800 µHz), because the amplitude of the acoustic modes decreases as the mode inertia increases as the frequency decreases, while the solar noise from incoherent, convective motions increases: thus the signal-to-noise ratio (SNR) of those modes is progressively reduced. Moreover, these lowfrequency p modes have very long lifetimes, as much as several years, which results in very narrow linewidths, hence precise frequency measurements.
Thanks to the long-duration helioseismic observations collected by the spaced-based instruments Michelson Doppler Imager (MDI) (Scherrer et al. 1995) and Global Oscillations at Low Frequencies (GOLF) (Gabriel et al. 1995) onboard the Solar and Heliospheric Observatory (SOHO) spacecraft, and by the ground-based, multi-site Global Oscillation Network Group (GONG) (Harvey et al. 1996) (Chaplin et al. 1996), the frequency resolution is continuously improving and the observation of lower radial-order solar p modes is becoming possible. Their precise mode parameter determination is of great interest for improving our resolution throughout the solar interior because they cover a broad range of horizontal phase velocity, and thus a broad range of depths of penetration. Moreover, these low-frequency modes have lower reflection points in the outer part of the Sun, which make them less sensitive to the turbulence and the magnetic fields in the outer layers, where the physics is poorly understood.
The usual mode-fitting analysis consists of fitting the 2l + 1 individual-m spectra of a given multiplet (n, l), either individually or simultaneously. Such fitting methods fail to obtain reliable estimates of the mode parameters when the SNR of the individual-m spectra is low. Instead, various pattern-recognition techniques have been developed in an effort to reveal the presence of modes in the low-frequency range (see e.g., Schou 1998;Appourchaux et al. 2000;Chaplin et al. 2002;Broomhall et al. 2007, and references therein). In the case of spatially-resolved helioseismic data (such as GONG and MDI observations), m-averaged spectra appeared to be a powerful tool, since for a given multiplet (n, l), there exist 2l + 1 individual-m spectra, which can result in an average spectrum with a SNR ≫ 1 once the individual-m spectra are corrected for the rotationand structure-induced frequency shifts. The m-averaged spectra were employed early in the development of helioseismology by Brown (1985), but were replaced by fitting the m spectra individually as the quality and the SNR of the data improved. However, years later, in order to fully take advantage of the long-duration helioseismic GONG and MDI instruments and reach lower frequencies in the solar oscillation spectrum, Schou (1998Schou ( , 2002Schou ( , 2004 and  Appourchaux et al. (2000) used the m-averaged spectra corrected by the modeled solar rotation to detect new low radial-order p modes and to set upper limits on the detectability of the g modes. These authors demonstrated the potential advantage of such rotation-corrected, maveraged spectra.
We present here an adaptation of the m-averaged spectrum technique in which the m-dependent shift parameters are determined by maximizing the quality of the resulting average spectrum. The analysis is performed on long-duration time series of the spatially-resolved helioseismic GONG and MDI observations of the low-and medium-angular degrees (1 ≤ l ≤ 35). This range of oscillation multiplets samples the radiative interior down to the solar core. In Sec. 2, we introduce the different datasets used in this analysis. In Sec. 3, we describe this new technique in order to observe low signal-to-noiseratio, low-frequency p modes, explaining the different steps of the analysis from the mode detection to peakfitting. In Sec. 4, we demonstrate that this method allows us to successfully measure lower-frequency modes than those obtained from classic peak-fitting analysis of the individual-m spectra by comparing with other measurements obtained from coeval datasets. In Sec. 5, we present the mode parameters of these long-lived, lowfrequency acoustic modes down to ≈ 850 µHz extracted from the analysis of 3960 days of GONG observations using the m-averaged spectrum technique. Finally, Sec. 6 summarizes our conclusions.

OBSERVATIONS
Details of the spatially-resolved helioseismic observations collected by both GONG and MDI used for this work (the starting and ending dates, and their corresponding duty cycles) are given in Table 1. Coeval 2088day observations of GONG and MDI were analyzed for oscillation multiplets with angular degrees from l = 1 to l = 35, and are then directly compared to those of Korzennik (2005) for l ≤ 25 measurements of the same datasets. We also applied the analysis to 3960 days of GONG data (1 ≤ l ≤ 35), which constitutes so far the longest time series (≈ 11 years, spanning most of solar cycle 23) of spatially-resolved observations analyzed.

METHOD
An m-averaged spectrum corresponds to the average of the 2l + 1 individual-m components of an oscillation multiplet (n, l), thus reducing the non-coherent noise. Before averaging, each m spectrum of a given mode (n, l) is shifted by a frequency that compensates for the effect of differential rotation and structural effects on the frequencies. The m-averaged spectrum concentrates, for a given multiplet (n, l), all of the 2l+1 m components, as it would be if the Sun were a purely-spherical, non-rotating object. Thus, the average of the 2l + 1 individual-m spectra considerably improves the SNR of the resulting m-averaged spectrum.

Determination of the shifts
The m-averaged spectrum is obtained by finding the estimates of the splitting coefficients, commonly called a-coefficients, which maximizes the likelihood of the maveraged spectrum. The a-coefficients are individually estimated through an iterative process, with the initial values taken from a model. Thus, for a given mode (n, l), the frequency shift δν nlm is parameterized by a set of coefficients, as: where a i (n, l) are the splitting coefficients, and P i l,m corresponds to the Clebsch-Gordan polynomial expansion as defined by Ritzwoller & Lavely (1991). In this definition, the odd orders of the a-coefficients describe the effects of solar rotation, while the even orders correspond to departures from spherical symmetry in the solar structure as well as to quadratic effects of rotation. Each a i is chosen to maximize the likelihood of the m-averaged spectrum. This is performed through an iterative procedure. For a particular order i of the coefficients a i , a range of values is scanned around its initial value, while the other a i ′ =i s are kept fixed to their previously estimated values.
For each scanned value of a i , the individual-m spectra are shifted by the corresponding Clebsch-Gordan polynomials, and the mean of these 2l + 1 shifted spectra is taken. The mean power spectrum is then fitted using a Maximum-Likelihood Estimator (MLE) minimization as described in Sec 3.2 and its likelihood determined.
For a Monte-Carlo simulation, the left panel of Fig. 1 shows the variation of the likelihood from the MLE minimization as a function of the first splitting coefficient a 1 , showing a well defined minimum which represents the best value of a 1 . The artificial power spectra were simulated following the methodology described in Fierry Fraillon et al. (1998). We also examined the sensitivity of the mode linewidth and the entropy as criteria for determining the best shifts. In our case, the entropy (Shannon 1948) can be seen as a measure of randomness in the m-averaged spectrum, S, and is defined as Both linewidth and entropy show well defined minima around the input value of a 1 (middle and right panels of Fig. 1 respectively). Indeed, the m-averaged spectrum gets narrower as a 1 converges to its input values of 0.4 µHz and a 1 = 400 nHz. Similar variations are obtained for all the a i s. As detailed in Appendix A, the use of these different criteria to determine the best estimates of the a-coeficients returns consistent results.
The iteration is performed until the difference between two iterations in each of the computed a i coefficients falls below a given threshold (such as 0.25σ in the case of a 1 ). Also, in order to remove any outliers, some quality checks are performed after each measure of an a i which needs to fall within a constrained range of values. For example, a ±15% window around its theoretical expectation is used for a 1 . Here, we fitted only the six first a i in the Clebsch-Gordan expansion, even though the quality of the data supports the determination of higher-order coefficients.
Finally, low SNR peaks in the m-averaged spectrum (after adjustment) are tested against the H0 hypothesis. In the framework of that hypothesis, the resulting spectra are tested against a statistics pertaining to pure noise (χ 2 with 2(2l+1) d.o.f). This test has been widely applied to helioseismic observations in the search for longlived, low radial-order p modes and g modes (see e.g., Appourchaux et al. 2000). In the present analysis, we rejected peaks that have a greater than 10% chance of being due to noise in the 238 analyzed windows, each containing 288 frequency bins. Here the fixed number of bins was chosen because we know that the range of theoretical frequency lie within 1.5 µHz or so. Figures 2 and 3 illustrate the advantage of using the m-averaged spectrum technique in the case of two oscillation multiplets for 2088 days of GONG data, where the m-averaged spectra before and after the correction for the splitting coefficients are shown. These examples show the m-averaged spectra of the modes l = 3, n = 5 at ≈ 1015.0 µHz (Fig. 2), and l = 16, n = 4 at ≈ 1293.8 µHz (Fig. 3), as well as the corresponding m − ν diagrams. These two examples were chosen to demonstrate the performance at different SNR levels. The corresponding 10% probability levels are given. The m − ν diagrams in the case of the mode l = 3, n = 5 (right panels on Fig. 2) do not show any high SNR structure before or after correction. However, the m-averaged spectrum after correction clearly shows the target mode (lower left-panel on Fig. 2), with an unambiguous detection level. The mode l = 16, n = 4 presents a higher SNR (Fig. 3) and its m − ν diagram shows that the individual-m spectra line up after correction (lower-right panel on Fig. 3). The estimated splitting coefficients of the low-frequency modes with 1 ≤ l ≤ 35 measured in the 3960-day GONG dataset are shown in Fig. 4 as a function of frequency and ν/L (with L = l(l + 1)), which is approximately propor-tional to the sound speed at the mode's inner turning point. Modes with selected ranges of radial orders are represented with different colors and symbols.

Extraction of the mode parameters
For a given mode (n, l), the best estimates of the splitting coefficients determined as discussed in Sec 3.1 are used to calculate its m-averaged spectrum. When N independent power spectra are averaged together, the statistics of the mean power spectrum corresponds to a χ 2 with 2 × N degree-of-freedom (d.o.f.) statistics. Appourchaux (2003) demonstrated that the mean of 2l + 1 independent power densities, which has a χ 2 with more than 2 d.o.f. statistics, can be correctly fitted with a Maximum-Likelihood Estimator (MLE) minimization code developed for spectra following a χ 2 with 2 d.o.f. statistics. The asymmetric Lorentzian model of Nigam & Kosovichev (1998) was used to describe the m-averaged spectrum, as: where x n,l = 2(ν − ν n,l ) Γ n,l .
Then, for a given mode (n, l), the central frequency, the Full-Width-at-Half-Maximum (fwhm), and the power height of the spectral density are respectively ν n,l , Γ n,l , and H n,l . The peak asymmetry is described by the parameter α n,l , while B n,l represents an additive, constant background level in the fitted window. The first spatial leaks (δl = 0, δm = ±2), commonly called m-leaks, are also included in the fitting model and added to Eq. 2. The frequencies of the m-leaks are set from the central frequency of the target mode using the previously measured splitting coefficients (Sec. 3.1). Their peak asymmetries are assumed to be the same as that of the target  -Dalsgaard et al. 1996). The dot-dashed lines on the left hand-side panels give the 10% probability limit that a peak is due to noise in the 238 windows, 1.5 µHz wide. The illustrated spectral window contains the 2l + 1 components of the represented multiplet.  -Dalsgaard et al. 1996). The dot-dashed lines on the left hand-side panels give the 10% probability limit that a peak is due to noise in the 238 windows, 1.5 µHz wide. The illustrated spectral window contains the 2l + 1 components of the represented multiplet.
mode, while their fwhms are a free parameter of the fit and different from the target mode. The amplitude of the m-leaks is specified to be a fixed fraction of the central peak, which was estimated from the leakage matrix developed especially for the GONG (Hill & Howe 1998) and MDI (J. Schou, private communication) data. The first spatial leaks in the m-averaged spectrum were determined by averaging for a given multiplet (n, l) the δm±2 leaks over the entire 2l + 1 spectra.
The size of the fitting window, Ω ν , is proportional to the first estimates of the mode width, Γ n,l , and centered around the frequency of the target mode. It is defined as: where ∆ν r is the frequency resolution of the power spectrum. The first spatial leaks are always included in the fitting range by adding the offset ∆ δm = 800 nHz. -First six a i splitting coefficients of the low-frequency p modes with 1 ≤ l ≤ 35 measured in 3960 days of GONG data. The a-coefficients are represented as a function of the mode frequency ν (left column), and ν/L, with L = p l(l + 1) (right column). The different colors and symbols correspond to selected ranges of radial orders n: green triangles, modes with n = 1, 2; purple stars, n = 3, 4; blue upside-down triangles, n = 5, 6; orange squares, n = 7, 8; and black dots, n ≥ 9.
The multiplicative factor 20 ensures a good sampling of the mode profile in the low-frequency range. A comparable definition of the fitting window was adopted by Korzennik (2005). Bad fits were removed based on a set of quality criteria based on the fitted mode parameters and associated uncertainties, such as, (1) the error of the mode frequency must be less than its mode width; (2) the SNR must be larger than 1; and (3) the mode width must be larger than the frequency resolution. A discussion on the impact of the fitting model (asymmetry, spatial leaks) on the extracted mode parameters used to describe the m-averaged spectrum can also be found in the Appendix B. Figure 5 shows examples of the m-averaged power spec-tra for four different radial orders n of the multiplet l = 17, and the corresponding best MLE fits, which included the mode asymmetry and the δm±2 spatial leaks. The blending of the first m leaks is particularly clear as the linewidths increase with increasing frequency.

Mode parameter and a-coefficient uncertainties
The mode parameter uncertainties are established in the usual manner by the inverse of the covariance matrix. However, because the m-averaged spectrum is fitted using a MLE minimization and, as explained in Appourchaux (2003), the formal uncertainties must be normalized by the square root of the number of averaged spectra, i.e., in our case, by √ 2l + 1. But this aposteriori error normalization is correct only if the 2l + 1 spectra of a given (n, l) mode have the same variance (or SNR). Since the condition of equal SNR among the m spectra within a multiplet is not satisfied in our case, the uncertainties of the mode parameters have to be taken as a first approximation only. However, Monte-Carlo simulations show that this error normalization holds even in the case of m-dependent SNR (see Sec. 3.3.2).
It can also be derived that the errors on the acoefficients can be estimated as follows: where i is the a-coefficient order and P i l,m the associated Clebsch-Gordan polynomials. The derivation of Eq. 5 is detailed in Appendix C.
3.3.1. m/l dependence of the signal-to-noise ratio Figure 6 shows the dependence in m/l of the SNR in the GONG data. This was obtained with modes observed in the 3960-day GONG dataset below 2000 µHz and of angular degree up to l = 35. Both mode amplitude and background noise depends on the azimuthal order m and can be described with polynomials with only even terms, the polynomials being different for both parameters. Note that any frequency dependence of the m/l dependence is averaged out in Fig. 6.
The m/l dependence of the SNR implies that the acoefficients are not exactly orthogonal and that their errors are correlated (see Appendix C). However, as a first order approximation, the errors on the a-coefficients can be estimated by using Eq. 5 (see Sec. 3.3.2).

Validation of the error estimates: Monte-Carlo simulations
The formal uncertainties of the mode parameters and of the a-coefficients were verified through Monte-Carlo simulations. The artificial power spectra were simulated following the methodology described in Fierry Fraillon et al. (1998). In a first series of simulated spectra, the m dependence in amplitude within a given multiplet (n, l) was introduced, the SNR being symmetric in |m| around the m = 0 spectrum. In a second series, no m dependence was introduced, i.e., a constant SNR over m. The mean values of the formal errors returned by the MLE minimization were compared to the RMS value of the corresponding fitted parameter. The Monte-Carlo simulations showed that in both cases the formal uncertainties of the m-averaged spectra determined as in Sec 3.3 using a MLE minimization are a very good approximation of the errors.

COMPARISON WITH OTHER MEASUREMENTS
4.1. Comparison with spatially-resolved observations (l ≤ 25) GONG and MDI use two independent peak-finding approaches to extract the mode parameters. Developed in the early 1990s, and mostly unchanged since, they provide mode parameters on a routine basis. Time series of 108 days are used by the GONG project (Anderson et al. 1990), while the MDI project uses 72-day time series (Schou 1999). Recently, Korzennik (2005) developed a new and independent peak-finding method of the individual-m spectra, optimized to take advantage of the long, spatially-resolved, helioseismic time series available today from both projects. Korzennik (2005) applied his peak-fitting to extract the low-and medium-degree (l ≤ 25) mode parameters from both GONG and MDI observations using one 2088-day long time series, as well as using five overlapping segments of 728 days. In order to compare our results obtained with the m-averaged spectrum technique, we applied the procedure described in Sec. 3 to the same 2088 days of GONG and MDI observations (Table 1). Figure 7 shows the l−ν diagrams of the low-frequency modes measured with the two different analyses in the case of the 2088-day GONG (left panel) and MDI (right panel) datasets. The modes measured by Korzennik (2005) with a classic peak-fitting method of the individual-m spectra are represented by the open circles. We considered that a given mode (n, l) from Korzennik (2005) was detected when at least two of the 2l + 1 m spectra were successfully fitted, which is enough to obtain estimates of the corresponding central frequency and first splitting coefficient a 1 . The red dots represent modes measured with the m-averaged spectrum technique which were not observed by Korzennik (2005). A significantly larger number of low-frequency modes in the 2088-day GONG and MDI datasets (respectively 45 and 14 new modes) down to ≈ 900 µHz can be measured using the m-averaged spectrum technique.

Mode parameter and uncertainty comparisons
In order to check the accuracy of the technique and to identify any potential bias in our analysis, we compare the central frequencies and splitting coefficients obtained by the two methods. The individual-m frequencies of Korzennik (2005) were fitted using a Clebsch-Gordan polynomial expansion (Ritzwoller & Lavely 1991) in order to estimate the corresponding central frequencies and a-coefficients of each (n, l) multiplet. The formal uncertainties of the individual-m frequencies were used as fitting weights. The left panel on Fig. 8 shows the distribution of the differences in central frequencies below ≈ 1800 µHz of the common modes between the 2088day GONG estimates measured using the m-averaged spectrum technique and from Korzennik (2005) (as represented on Fig. 7), demonstrating that there is no frequency dependence over the analyzed low-frequency  range. The distribution was fitted by a Gaussian function, and its associated parameters (mean, standard deviation) are indicated on Fig. 8. While, on average, the GONG central frequencies obtained using the maveraged spectrum technique are less than 1 nHz smaller than Korzennik (2005)'s estimates, this offset is not significant -the corresponding standard deviation being about 5 times larger. The MDI frequencies estimated with the m-averaged spectrum technique give comparable, insignificant mean differences with Korzennik (2005). Similar results are obtained with the splitting coefficients.
We also compared the low frequencies (ν 1800µHz, see Fig. 7) estimated in both the 2088-day GONG and MDI datasets using the m-averaged spectrum technique. The right panel of Fig. 8 represents the distribution of the frequency differences of the common modes, in the sense GONG minus MDI. The mean difference value is of -0.17 ± 1.99 nHz, i.e. the GONG and MDI lowfrequency modes are essentially the same. The mean difference in Korzennik (2005)'s central frequencies between the 2088-day GONG and MDI datasets for modes below 1800 µHz is of 0.35 ± 5.40 nHz. The splitting coefficient estimates are also consistent between the two datasets with in the case of the a 1 coefficient a mean difference of -0.04 ± 0.31 nHz.
The left panel of Fig. 9 shows the formal 1σ uncertainties of the central frequencies of the measured common modes between the 2088-day GONG dataset analyzed in the present analysis and the coeval 2088-day GONG dataset from Korzennik (2005) up to ≈ 1800 µHz. Our estimates of the frequency uncertainties are much smaller than those quoted by Korzennik (2005). However, Fig. 10 in Korzennik (2005) suggests that the errors are overestimated, and that his results "might be too conservative". Korzennik (2008) reported that in the case of a 2088-day long time series, as a first estimate, a multiplicative factor of 0.75 needs to be applied to the frequency uncertainties reported in Korzennik (2005). However, despite these uncertainty scaling issues, while the Korzennik (2005)'s uncertainties show an increase with decreasing frequency from ≈ 1500 µHz, the uncertainties returned from the maveraged spectrum technique do not show this increase, thanks to the higher SNR of the m-averaged spectrum than for the individual-m spectra.
The uncertainties on the a-coefficients returned by the m-averaged spectrum technique are also smaller than the ones obtained by fitting the individual-m spectra, as shown on the right panel of Fig. 9 in the case of the a 1 coefficients. As for the frequencies, the a-coefficients of the modes below ≈ 1500 µHz are better constrained using the m-averaged spectrum technique.

4.2.
Comparison with Sun-as-a-star observations (l ≤ 3) The spatially-resolved GONG and MDI instruments are not optimized to observe low-degree solar p modes below l ≤ 3, unlike the Sun-as-a-star, integrated-light instruments such as the space-based instrument GOLF onboard SOHO and the ground-based, multi-site BiSON network. The low-degree modes are of particular interest as they reach the very deep interior of the Sun. However, the spatially-resolved observations are still able to observe such low-degree oscillations.
Low-degree (l ≤ 3) modes down to ≈ 1000 µHz are observed in both GONG and MDI data with the maveraged spectrum technique as illustrated in Figs. 7 and  (2005) using a classic peak-fitting method of the individual-m spectra, while the red dots correspond to the additional modes measured using the m-averaged spectrum technique that were not observed by Korzennik (2005). The ridges of same radial order are also indicated from n = 2 to n = 12. 12. In order to test the capability and the precision of the m-averaged spectrum technique to observe lowdegree, low-frequency modes in spatially-resolved data, measurements obtained for ≈ 11 years of the Sun-as-astar GOLF and BiSON instruments were compared with the 3960-day GONG dataset and the 2088-day GONG and MDI datasets. The GOLF data were independently analyzed by two mode-fitting algorithms (R. A. García, private communication 4 , and P. Boumier, private communication 5 ). The BiSON observations come from a combination of frequencies obtained with different long time series analyzed in order to measure low-frequency modes (W. J. Chaplin, private communication. See also Chaplin et al. 2002;Broomhall et al. 2007). However, while imaging instruments give us access to all of the 2l + 1 individual-m components, only l + 1 components can be clearly observed with integrated-sunlight technique. These various components have different spatial structure over the solar sur- face, which can lead to differences in the extracted central frequencies of the multiplet. Chaplin et al. (2004) and  derived expressions to predict the differences between the low-degree frequencies extracted from spatially-resolved (as MDI and GONG) and Sun-as-a-star (as GOLF and BiSON) observations. However, in the following, as a first approximation, we compared directly the extracted mode parameters.
The comparisons of the estimated mode frequencies and a 1 rotational splittings between the common lowdegree (1 ≤ l ≤ 3), low-frequency modes in the two types of observation are shown on the left and right panels respectively of Fig. 10. The three different datasets and analysis methods give consistent results, for both the frequency and the splitting coefficient a 1 . Of course, this is only assuming that the different subsets of observed multiplets from both types of observational technique "see" the same central frequencies.
Thanks to decade-long available datasets, the lowdegree, low-frequency modes are today measured lower than 1200 µHz with high precision, demonstrated by  Korzennik (2005) fitting the individual-m spectra (red dots) and by using the m-averaged spectrum technique (black plus signs).  ≤ 3), low-frequency p modes observed in spatially-resolved GONG and MDI data using the m-averaged spectrum technique and in two Sun-as-a-star, full-disk observations, GOLF and BiSON, over comparable periods of time. The GOLF data were analyzed by two independent peak-fitting algorithms (R. A García, filled symbols; P. Boumier, crosses). The comparison with BiSON is illustrated with the open symbols. The frequency differences were normalized by the combined errors. Right panel -Same as on the left panel, but for the splitting coefficient a 1 . the consistency in the extracted parameters from different instruments using distinct and independent analysis. Figure 10 also demonstrates that spatially-resolved observations can provide as accurate measurements of the low-degree modes as the Sun-as-a-star instruments do. Moreover, the m-averaged spectrum technique allows the observation of lower radial-order l = 3 modes than the integrated-light GOLF and BiSON observations, for commensurate observation lengths, thanks to the observations of the 2l + 1 components (Fig. 11).

MODE PARAMETERS OF THE LOW-FREQUENCY OSCILLATIONS
The m-averaged spectrum technique has been applied to 3960 days of GONG observations (see Sec. 2), spanning most of the 11 years of solar cycle 23. The analysis covered low-frequency modes with angular degrees from l = 1 to l = 35. Oscillation multiplets well below 1000 µHz were detected with good precision, such as the modes l = 4, n = 4 at ≈ 913.5 µHz; l = 9, n = 3 at ≈ 930.5 µHz; l = 16, n = 2 at ≈ 912.1 µHz; or l = 31, n = 1 at ≈ 907.5 µHz. Some examples are illustrated in Figs. 2 and 3. These low horizontal-phase-velocity modes do not penetrate deeply into the Sun, but their very high inertias afford higher precision frequencies for the inversions. It is clear from Sec. 4 that this method allows us to observe modes that are otherwise lost in the background of each individual-m spectrum of a given multiplet (n, l), and thus unobservable with a classic peak-fitting analysis. The l − ν diagram of the observed low-frequency modes (1 ≤ l ≤ 35) down to n = 1 and ≈ 850 µHz in the 3960-day GONG dataset and 2088-day GONG and MDI datasets with the m-averaged spectrum technique is shown on Fig. 12. Figure 13 shows the fitted mode fwhms Γ n,l (upperleft panel) and mode heights H n,l (upper-right panel) of the measured low-frequency oscillations. The fitted background level is also represented on the right panel. The fwhms and heights are extremely valuable tests of models of the physical processes responsible for the mode damping and excitation by the turbulent convective motions in the outer layers of the Sun: the mode damping is inversely related to the fwhm of the mode, and the mode excitation is proportional to the mode height times the mode fwhm squared (for a detailed description, see, e.g., Salabert & Jiménez-Reyes 2006). The leveling off of the mode widths observed below ≈ 1100 µHz, despite their dispersion becoming larger, could be a resolution effect, the peaks being then so narrow that the limiting resolution of the spectrum becomes an issue. Moreover, Schou (2004) did not observe such behaviour at low frequency in MDI data with a 2952-day time series.

Mode linewidths, heights, and background levels
As indicated by different colors and symbols on Fig. 13, the fitted mode widths follow ridges for equal radial orders n. This dependence on angular degree (l) is directly related to the mode inertia (I) in terms of a power law, as illustrated on the lower-left panel of Fig. 13. The l dependence in the mode fwhms is removed when represented as a function of the mode inertia I.

Mode asymmetry
The mode parameters extracted through the routine GONG and MDI peak-fitting pipelines are obtained by use of symmetric Lorenztian profiles (Anderson et al. (1990) and Schou (1992) respectively). However, it was demonstrated that ignoring the peak asymmetry in the description of the acoustic modes leads to biais in the estimated mode parameters (see Appendix B and Thiery et al. 2000). Today, most of the estimates of the mode asymmetries have been restricted to low degrees (l ≤ 3) only, from Sun-as-a-star, integrated-sunlight observations. However, Korzennik (2005) used asymmetric profiles and presented estimates of the peak asymmetry for modes with angular degrees 1 ≤ l ≤ 25, obtained with GONG and MDI observations. Recently, Larson & Schou (2008) are planning to reprocess all the MDI medium-l data including a set of corrections and improvements (such as the mode asymmetry) in the MDI pipeline algorithm itself.
The asymmetry parameter (α n,l ) in the low-frequency range, obtained by fitting the 3960-day GONG maveraged spectrum (1 ≤ l ≤ 35), is shown on the lowerright panel of Fig. 13. The extracted peak asymmetry is well constrained down to ≈ 1400 µHz, with a mean value of about −0.044 ± 0.002, and no discernable l dependence. The average asymmetry observed in the m-averaged spectrum is consistent with other measurements. For instance, the mean value observed by Korzennik (2005) was about −0.04 for modes below 2000 µHz and l ≤ 25, once his estimates are transformed back into the Nigam & Kosovichev (1998)'s definition of the peak asymmetry. A comparable mean value is also observed at the lowest frequencies for which asymmetries were reported in Sun-as-a-star, integrated-sunlight observations (e.g., Thiery et al. 2000). Figure 14 shows the frequency differences (in µHz) between the fitted low-frequency modes observed in the 3960-day GONG dataset using the m-averaged spectrum technique and the corresponding theoretical values calculated from Christensen-Dalsgaard's model S (Christensen-Dalsgaard et al. 1996). The corresponding frequency uncertainties were multiplied by 20 to render them visible. These comparisons are represented as a function of the angular degree (left panel), of the frequency (middle panel), and of the inner turning point (right panel). Modes of equal radial orders are connected. As these differences between observed and theoretical frequencies show, there is still room to improve the model of solar internal structure. Note that the right panel on Fig. 14 illustrates also the wide range of depths of penetration that these low-frequency modes cover.

CONCLUSION AND DISCUSSION
We presented here an adaptation of the rotationcorrected, m-averaged spectrum technique to observe low signal-to-noise-ratio, low-frequency solar p modes in spatially-resolved helioseismic data. For a given multiplet (n, l), the shift coefficients describing the differential rotation-and structural-induced effects are chosen to maximize the likelihood of the m-averaged spectra. The average of the 2l + 1 individual-m spectra can result in a high signal-to-noise ratio when the individual-m spectra have a too low signal-to-noise ratio to be successfully fitted. This technique was applied to long time series of the spatially-revolved GONG and MDI observations for low-frequency modes (i.e., approximately below 1800 µHz) with low-and intermediate-angular degrees (1 ≤ l ≤ 35). We demonstrated that it allows us to measure lower frequency modes than with classic peak-fitting analysis of the individual-m spectra. Figure 15 shows the new low-frequency solar p modes observed in spatiallyresolved data using the m-averaged spectrum technique in long time series of both GONG and MDI observations. Their central frequencies and splitting a 1 coefficients, as well as their associated uncertainties are indicated in Table 2. These normal modes of oscillation were predicted but were not measured previously. The potential of the m-averaged spectrum technique for increasing our knowledge of the solar interior is clearly illustrated on Fig. 15, p modes well below 1000 µHz being measured with a high accuracy thanks to their longer lifetimes. We also demonstrated that the m-averaged spectrum technique returns unbiased results with no systematic differences with other long-duration measurements, which also include the asymmetry in the mode profile description.  and height H n,l (upper-right panel) as a function of frequency for the low-frequency p modes l ≤ 35 extracted from 3960 days of GONG data with the m-averaged spectrum technique. The different colors and symbols correspond to selected ranges of radial orders n: green triangles, modes with n = 1, 2; purple stars, n = 3, 4; blue upside-down triangles, n = 5, 6; orange squares, n = 7, 8; and black dots, n ≥ 9. The frequency resolution of the analyzed dataset is indicated on the upper-left hand-side plot, and the fitted background level (B n,l ) is also represented on the upper-right hand-side plot (small black dots). Lower-left panel -Same as above, but for the mode fwhm Γ n,l as a function of the mode inertia (I), calculated from Christensen-Daslgaard model S (Christensen-Dalsgaard et al. 1996). The frequency resolution is also indicated. Lower-right panel -Same as above, but for the mode asymmetry (α n,l ).
The oscillation parameters of these low signal-to-noise-ratio, low-frequency modes, such as their central fre-  quencies, splittings, asymmetries, lifetimes, and heights were measured. These low-frequency p modes contribute to improve our resolution throughout the solar interior since they sample a large range of penetration depths. Moreover, because these modes have lower upper turning points in the outer part of the Sun, they are less sensitive to the turbulence and magnetic fields in the outer layers, which should make them extremely valuable for the study of the physical processes responsible for the oscillation excitation and damping by the turbulent convection.
We would like to recall that Schou (1992)'s peakfinding approach consists in fitting the individual-m spectra simultaneously by using a model in which the shift coefficients are introduced, while in the present technique, the best shifts are determined first, based on the calculation of figure-of-merits (Sec. 3.1 and Appendix A), and then the rotation-corrected, m-averaged spectrum is fitted (Sec. 3.2).
The development of the m-averaged spectrum technique towards both higher frequencies and larger angular degrees is one of the next step to be addressed, as also the analysis of shorter datasets, such as the canonical 108-and 72-day time series.
This work utilizes data obtained by the Global Oscillation Network Group (GONG) program, managed by the National Solar Observatory, which is operated by AURA, Inc. under a cooperative agreement with the National Science Foundation. The data were acquired by instruments operated by the Big Bear Solar Observatory, High Altitude Observatory, Learmonth Solar Observatory, Udaipur Solar Observatory, Instituto de Astrofísica de Canarias, and Cerro Tololo Interamerican Observatory. The GOLF and MDI instruments onboard SOHO are cooperative efforts to whom we are indebted. SOHO is a project of international collaboration between ESA and NASA. BiSON is funded by the Science Technology and Facilities Council (STFC). We thank the members of the BiSON team, and colleagues at the host institutes at each of the BiSON sites. The authors are particularly grateful to S. G. Korzennik for providing us with the 2088-day MDI dataset, and to J. Schou for the MDI leakage matrix. The authors thank R. A García, P. Boumier, and W. J. Chaplin for providing estimates of GOLF and BiSON mode frequencies and a 1 rotational splittings observed in decade-long time series. The authors also thank S. J. Jiménez-Reyes and J. Schou for their useful comments on the manuscript, and S. G. Korzennik for helpful discussions during the various stages of this work. D. S. acknowledges the support of the NASA SEC GIP grant NAG5-11703. This work has been partially funded by the grant PNAyA2007-62650 of the Spanish National Research Plan.

FIGURES-OF-MERIT AND DETERMINATION OF THE A-COEFFICIENTS
The best estimates of the splitting a-coefficients are obtained by maximizing the likelihood of the m-averaged spectrum (see Sec. 3.1) through the calculation of a figure-of-merit (FOM). However, as shown on Fig. 1, other criteria to define a FOM can be used, such as the narrowest peak (i.e., the minimum mode linewidth), or the minimum entropy of the resulting m-averaged spectrum. In order to compare the actual mode parameters and associated uncertainties obtained with two different definitions of the FOM, we applied the m-averaged spectrum technique to the 3960-day GONG dataset by using both the maximum likelihood and the narrowest peak in the m-averaged spectrum as FOM. Figure A16 shows the corresponding central frequencies, and the odd a 1 , a 3 , and a 5 splitting coefficients of the common, measured low-frequency p modes. The associated formal uncertainties are also represented. The two FOMs return consistent mode parameters within the error bars, the difference between the two being within the 3σ limit for all of the mode parameters.

IMPACT OF THE FITTING MODEL USED
As a test of the dependence of the measured frequencies on the fitting model used to describe the m-averaged spectra, we fitted the m-averaged spectra using three different models: an asymmetric Lorenztian profile (Eq. 2) including the closest δm ± 2 spatial leaks (hereafter called A2 and used as the reference model); a symmetric Lorenztian profile including the δm ± 2 spatial leaks (hereafter S2); and an asymmetric Lorenztian profile (Eq. 2) but omitting the neighbouring δm ± 2 spatial leaks (hereafter A). Figure B17 shows the differences as a function of frequency in the 2088-day GONG low frequencies estimated using the m-averaged technique between S2 and A2 (red dots), and between A and A2 (black plus signs), in both cases A2 being the reference model. Ignoring the peak asymmetry in the fitting model leads to a systematic underestimation of the mode frequency as the frequency increases, the effect becoming particularly large above ≈ 1400 µHz (red dots). The differences become much larger than 3σ, for example, at ≈ 1800 µHz, the fitted frequencies between S2 and A2 are about 20σ apart. These results obtained for modes below 2000 µHz confirm previous observations, e.g. Thiery et al. (2000) who analyzed low-degree modes above 2000 µHz in 805 days of GOLF data.
On the other hand, omitting the spatial leaks has no effect below ≈ 1600 µHz, as they become well separated from the main peak because the corresponding mode linewidths are much smaller than their frequency separation. As the frequency increases, the mode linewidths increase, and ignoring the spatial leaks in the fitting model of the m-averaged spectrum between about 1600 and 2000 µHz leads to an underestimation of the target mode frequency, the maximum difference occuring around 1800 µHz. The frequency separation between the target mode and the mleaks then becomes comparable to their linewidths and the lines blend together in the m-averaged spectrum. Above 2000 µHz, this underestimation seems to vanish. Indeed, at that frequency range, the mode linewidths are much larger than the frequency separation, and the first spatial leaks (at least) are totally blended into the target mode in the m-averaged spectrum, having a much lower impact on the frequency determination. However, the effect of ignoring the peak asymmetry is much larger than that from ignoring the m leaks even in the frequency range where the m leaks have the strongest impact. For instance, at 1800 µHz, the effect on the frequency underestimation by ignoring the mode asymmetry is about seven times larger than by ignoring the m leaks.
As an example of the other mode parameters, the right panel of Fig. B17 shows the differences in the extracted mode linewidths between the different fitting models. The color code is the same as for the differences in frequency represented on the left panel of Fig. B17. Ignoring the presence of the m leaks in the fitting model leads to a 35% overestimation at most of the extracted linewidths in the low-frequency range showing a maximum mismatch around 1900 µHz. Interestingly, if the m leaks are omitted, the linewidths are underestimated below ≈ 1600 µHz, showing a maximum 10% underestimation around 1500 µHz. On the other hand, ignoring the peak asymmetry has a very small influence on the fitted linewidths in the low-frequency range. However, above ≈ 1800 µHz, the linewidths extracted using an asymmetric profile are systematically larger than the ones returned using a symmetric profile. -Left columnl − ν diagram, and a 1 , a 3 , and a 5 splitting coefficients for the common, low-frequency p modes measured in the 3960-day GONG dataset obtained by using the maximum likelihood (circles) and the narrowest peak (dots) of the m-averaged spectrum as figure-of-merit to determine the best estimates of the a-coefficients (see Sec. 3.1 and Fig. 1). Right column -Associated 1σ formal uncertainties (nHz) of the mode central frequencies, and of the a 1 , a 3 , and a 5 splitting coefficients.

DERIVATION OF ERRORS FROM AN M -AVERAGED SPECTRUM
The derivation of the errors of the mode central frequencies and of the a-coefficients measured from the m-averaged spectrum technique is detailed here.
Approximation of the statistics of the m-averaged spectrum The m-averaged spectrum is obtained from the summation of 2l + 1 spectra assumed to be with χ 2 with 2 d.o.f statistics each having a different mean or signal-to-noise ratio. All of the individual m-spectra are independent from each other. The solar background noise is assumed to depend on m with a polynomial with only even terms (0, 2, etc...). The amplitude of the modes is assumed to depend on m with a different polynomial also with even terms (0, 2, etc...). In a first step, the a i -coefficients are calculated to maximize the likelihood of the resulting m-averaged spectrum.
Using Appourchaux (2003), we can derive an approximation of the statistics of the summation of the 2l + 1 spectra. The statistics of the m-averaged spectrum S can be approximated by a Gamma law given by: B17.-Effect of asymmetry and spatial leaks on the fits. Left panel -Differences (in nHz) in the 2088-day GONG frequencies estimated using different models to describe the m-averaged spectrum: symmetric Lorentzian profile minus asymmetric Lorentzian profile (Eq. 2) both models including the first δm ± 2 spatial leaks, i.e., S2 minus A2 (red dots); both asymmetric Lorentzian profiles (Eq. 2) but ignoring the closest δm ± 2 spatial leaks for one of them, i.e., A minus A2 (black plus signs). The differences in the extracted frequencies using different fitting profiles are represented as a function of frequency. Right panel -Same as the left panel, but for the differences (in µHz) in the fitted mode linewidths.
which shows that using the MLE applied to a χ 2 with 2 d.o.f as prescribed by Appourchaux (2003) is in the case of the m-averaged spectrum a good approximation. It is not an approximation when averaging several power spectra of identical mean (or variance), i.e. when α = β = ρ = 0. Note that what we minimize is the sum over a range of frequency that can be approximated as: Error bars on the central frequencies Error bars for frequency are derived from the inverse of the Hessian (second derivative of L) as: Toutain & Appourchaux (1994) showed that we could express the error bars as a function of the mode profile P (=A + B) as: σ −2 ν0 = T (2l + 1) 1 P 2 (ν) where T is the observation time. The 2l + 1 factor is due to the fact that the likelihood is 2l + 1 times larger than the likelihood of Toutain & Appourchaux (1994) (cf Eq. 23). Eq. (C26) shows that the error bars on the frequencies in the m-averaged spectrum will be √ 2l + 1 smaller than for the mean of the individual modes. In deriving Eq. (C26), we assumed that u = P. This is an approximation good enough for getting the error bars on the frequency but not on the a i .

Error bars on the a-coefficients
The error bars on the a-coefficients are derived from the inverse of the Hessian (second derivative of L) as: As shown by Toutain & Appourchaux (1994), these coefficients can be related to the mode profile as using Eq. (C23): using the following property: where P i l,m (m/l) are the Ritzwoller-Lavely polynomials. And we finally get: h ij = T l 2 m P i l,m (m/l)P j l,m (m/l) We recognize the error bars of the frequency for the m spectrum depending on the inverse of the signal-to-noise ratio β m (as in Libbrecht 1992): Finally Eq. (C30) can be written as: The errors σ a1 scale like l − 3 2 (as in Veitzer et al. 1993). If the SNR is the same for all m, then we have σ −2 ν0 = (2l+1)σ −2 m . Thus, by simply using the orthogonality property of the P i l,m polynomials, and as given in Sec. 3.3 (Eq. 5), we obtain the following expression to calculate the error bars of the a-coefficients in the m-averaged spectrum: σ −2 ai = l 2 2l + 1 m P i l,m (m/l) 2 σ −2 ν0 .
All terms off of the diagonal are zero. Of course, when the SNR varies with m, the off-diagonal terms are non-zero and correlations appear.