Low-degree solar rotational splitting from 45 yr of BiSON observations

ABSTRACT

We present solar low-degree rotational splitting values based on a new analysis of Sun-as-a-star observations from the Birmingham Solar Oscillations Network, covering a 16 425-d period from 1976 December 31 to 2021 December 20 with a duty cycle of 57 per cent. The splitting values are estimated from the power spectrum using a Markov chain Monte Carlo sampling method, and we also present for comparison the results from an analysis of 100 realizations of synthetic data with the same resolution and gap structure. Comparison of the scatter in the results from the synthetic realizations with their estimated uncertainties suggests that for this data set the formal uncertainty estimates are about 30 per cent too small. An upward bias in the splittings at frequencies above 2200 μHz, where the components are not fully resolved, is seen in both the observed and synthetic data. When this bias is taken into account, our results are consistent with a frequency-independent synodic rotational splitting value of 400 nHz.

1 INTRODUCTION

The low-degree solar p modes are one of the very few tools available to probe the structure and rotation of the deepest parts of the solar interior. In order to use the information they provide, we need both to determine their properties as accurately as possible and to have realistic estimates of the uncertainties on these measurements. The Birmingham Solar Oscillations Network (BiSON) has monitored these oscillations since the mid-1970s, originally in short observing campaigns with one or more ground-based instruments and later from an automated six-site network, giving us the longest time series of such observations available to date. In this article, we present solar rotational splitting values obtained by fitting a single Fourier power spectrum covering almost the entire lifetime of the project to date.

Solar rotation lifts the degeneracy between spherical harmonics with the same degree l > 0 and different azimuthal order m, giving rise to a multiplet in which the variation of the frequency with m contains information about the rotation and asphericity. The solar rotation period of approximately 27 d corresponds to an angular frequency of about 430 nHz, but because the Earth orbits the Sun once a year, which translates to an angular frequency of 31.7 nHz, the ‘synodic’ rate that is measured from Earth-based observations is approximately 400 nHz. This means that the difference in frequency between the m = −l and m = l components at a given degree l and radial order n is about l × 800 nHz. Fig. 1, which is based on the first-guess table discussed in Section 3.5, shows how the m = ±l splittings for l = 1, 2, and 3 relate to the mode width; because peaks need to be separated by at least twice their full width at half-maximum Γ to avoid a biased or unstable fit when using unconstrained maximum likelihood estimation fitting (see e.g. Howe & Thompson 1998), we show both the Γ and 2Γ values as a function of frequency. The l = 1 components are separated by 2Γ only below about 1800 μHz, while the outer components of l = 3 are 2Γ apart up to 3200 μHz. The figure also shows the separation between frequencies for the (l + 2, n − 1)/(l, n) mode pairs for l = 0 and l = 1; we can see that at frequencies above about 3600 μHz for l = 2/0 and 3800 μHz for l = 3/1 the mode pairs are not resolved, so it would be challenging if not impossible to estimate the splittings of individual modes in this regime to any meaningful precision.

BiSON data have been used to estimate low-degree rotational splittings since shortly after the deployment of the worldwide network. Elsworth et al. (1995) used BiSON data from 1992 January to 1994 August, divided into 16-month segments, to measure the rotational splittings of l = 1 modes between 1473 and 2559 μHz (radial order n from 9 to 17), of l = 2 modes between 1536 and 2486 μHz (9 ≤ n ≤ 16), and of l = 3 modes from 2138 to 2676 μHz (13 ≤ n ≤ 17). Chaplin et al. (2001) analysed an 8-yr set of BiSON data from 1992 January to 1999 December; they give splittings down to l = 1, n = 9 (1473 μHz), with the lowest frequency l = 1 splitting having an uncertainty of 5.3 nHz.

Davies et al. (2014a) obtained splittings and other mode parameters for modes at frequencies down to 1185 µHz (n = 7), using a Bayesian algorithm to fit to a spectrum based on 22 yr of BiSON data from 1991 January to 2012 December. They cite a ‘duty cycle efficiency’ for this data set of 68.3 per cent in a time series prepared to optimize for low-frequency noise.

García et al. (2008) analysed a time series of 4182 d from 1996 April 11 to 2007 September 22 from the Global Oscillations at Low Frequencies (GOLF) instrument onboard the Solar and Heliospheric Observatory, with 94 per cent duty cycle. They were able to measure splittings down to the l = 1, n = 7 (1185 μHz) mode. They cite Chaplin et al. (2006) on the possible bias of splittings by leakage from higher degree modes (e.g. l = 4 impinging on the l = 2 fitting window), and hence they use a 45 μHz window to fit each mode pair. They state that the l = 1 splittings are ‘roughly constant’ for the smaller fitting windows up to 3400 μHz but increase at higher frequencies than that, whatever window is used, so they trust the splittings only up to 3400 μHz.

2 DATA

The main data set that we analyse here is based on nearly 45 yr of data from BiSON, prepared as described by Davies et al. (2014b). The time series is zero-padded to give an integer multiple of 365 d and covers 45 × 365 = 16 425 d from 1976 December 31 to 2021 December 20, with an overall duty cycle of 57 per cent. The performance of the network in its early years has been described by Chaplin et al. (1996), and a more current overview was given by Hale et al. (2016).

3 ANALYSIS

We form the acoustic power spectrum by carrying out a fast Fourier transform of the prepared time series between the selected dates, in which any missing data are replaced by zeros in the detrended time series.

3.1 Peak profile

We fit the spectrum using a model in which each mode or rotationally split component is an asymmetric Lorentzian function of the frequency ν, described by the formula

where

ν₀ is the frequency of the Lorentzian component, Γ its width, h its height, and b is a fractional parameter characterizing the asymmetry. The expression simplifies to the normal Lorentzian for b = 0. This is the formulation proposed by Nigam & Kosovichev (1998), but the quadratic term in b is suppressed, as proposed by Fletcher et al. (2009), in order to ensure that the value of the expression is small far from the central frequency.

3.2 Rotational multiplets

Because rotational splitting lifts the degeneracy between modes of the same l and different m, for each (l, n) there are potentially 2l + 1 components of different m. As the inclination of the Sun’s rotation axis to the observer is close to zero, we assume that components with l − m odd have zero amplitude: Davies et al. (2014a), who also neglect these components, estimate their size at less than 3 per cent of the |m| = l power. Therefore, in practice, we only deal with l + 1 components for a mode of degree l: m = ±1 for l = 1, m = 0, m = ±2 for l = 2, and m = ±1, m = ±3 for l = 3. Furthermore, we assume that the frequency and amplitude of components within a multiplet are symmetric about m = 0, and that components m, m′ within a multiplet are separated by (m − m′)δΩ, where δΩ is our ‘splitting’ parameter; following Davies et al. (2014a), as we are concerned only with low-degree modes we do not take into account differential rotation or asphericity. The asphericity term for l = 2 is just detectable in BiSON data at epochs of high solar activity (Chaplin et al. 2003), but while it might bias the frequency measurement it should not bias the splitting. For the l = 3 multiplet, because of differential rotation, measuring the splitting in this way is not quite equivalent to measuring the first-order term of the polynomial expansion of frequency as a function of m, but the difference (which we estimate at around 5  nHz) will be within the uncertainties for all but the lowest frequency modes, and for those the signal-to-noise ratio of the inner components is low, so it is not practical to fit the differential rotation term. The heights of the m = 0 component of the l = 2 multiplet and the m = ±1 components for l = 3 modes are scaled by a visibility factor V_l relative to the m = ±l component; we base the prior distributions for these factors on values of 0.54 for V₂ and 0.38 for V₃. These values correspond to the ones used to construct synthetic BiSON-like data as described by Chaplin et al. (2006), based on the fitting results of Chaplin et al. (2001). The asymmetry and width parameters are also assumed to be the same for all components within a multiplet. The full model is built up by adding the asymmetric Lorentzian functions for each component in an l = 2/0 or l = 3/1 pair to a flat background offset c.

3.3 Window-function convolution

Ground-based observations, even with a network, tend to have a daily periodicity in their observing window; the solar spectrum is convolved with the power spectrum of the window function (Lazrek & Hill 1993). This results, among other effects, in ‘sidelobes’ at 1 d⁻¹ (11.57 μHz) on either side of each solar mode, which for some orders coincide with the other mode in an (l, n)/(l − 2, n + 1) pair. To mitigate this, as the last step of evaluating the model function we convolve the model with the power spectrum of the window function (a sequence of ones and zeros of the same length as the data, where a non-zero value corresponds to data being present), truncated at ±100 μHz. The convolution is implemented such that the model outside the fitting window is assumed to take the background offset value c, to avoid step-function effects in the convolved model that would occur if it were set to zero. To be specific, we use the astropy.convolution.convolve_fft function and set the fill_value keyword to c.

3.4 Fitting window

We fit the modes in pairs, (l + 2, n − 1)/(l, n), where l is 0 or 1. To specify the fitting window, we first take the mid-point of the mode pair and then choose the upper and lower limits to be at least 22 μHz above and below this point. We also require that the central frequency of each multiplet in the range should be no less than 15 μHz from the end of the range, extending the range if necessary to ensure this. In practice, this means that the l = 2/0 pairs are fitted in a 44 μHz window and the l = 3/1 pairs, which are more widely separated, in a window of about 50–55 μHz. By using this relatively narrow window, we avoid interference from the weak l = 4 mode, which typically appears about 25 μHz below the l = 2 peak, as shown in Fig. 2.

3.5 Optimization

We use the affine-invariant sampler from the python emcee package (Foreman-Mackey et al. 2013) to perform a Markov chain Monte Carlo (MCMC) sampling of the model posterior distribution given the observed spectrum. The spectral density is |$\chi ^2_2$| distributed, and so the log-likelihood function is given by (e.g. Anderson, Duvall & Jefferies 1990)

where the sum is over the N frequency bins in the range around each mode pair, M is the model, and S is the observed spectrum. Note that strictly this equation only applies for data without gaps; introducing gaps lowers the effective resolution of the power spectrum and means that the bins are not independent.

The priors we apply to each of the model parameters are presented in Table 1. We use prior distributions based on the Gaussian function for all of the parameters of our model, specified by a centroid value μ and a width σ. Hard constraints – where the log-probability goes to −∞ if the limit is exceeded – are used to keep the splitting and the visibility factor positive, modifying the underlying Gaussian prior for those parameters. The centroids of the prior distributions for the frequency (ν) and width (Γ) parameters are taken from a first-guess table in which the frequencies are based on the results of Broomhall et al. (2009) and the mode widths are from a smoothed version of earlier historical fits to low-degree modes. For amplitude, line width and background offset terms we vary the logarithm of the parameter, so it is the logarithmic value that is drawn from the distribution. Some additional prior constraints are included in the calculation of the prior probability function in order to keep the widths within a mode pair reasonably close, and to ensure that the asymmetry values within a pair are strongly correlated. A further hard constraint on the difference between frequencies in a mode pair is applied to prevent the modes in the pair from swapping places, which can otherwise occasionally happen when the mode width is approaching the separation between modes.

We use 100 walkers over 2000 steps and discard any walkers with an acceptance fraction below 0.16. While there is no formal criterion for the convergence of an MCMC fit, from visual inspection of the evolution of the parameters during the fit we found that the parameters usually settle around their final value within the first 500 steps.

The full information about the fit result is contained in the posterior distributions of the parameters, and we illustrate some of these below in Section 4. We use the median and half of the difference between the 16th and 84th percentile of the marginalized posterior distributions as summary statistics on each of the model parameters when a single uncertainty is needed, while the positive and negative uncertainties shown in Table 2 come from the difference between the median and the 84th and 16th percentiles, respectively. In most cases, the posteriors for the splitting are Gaussian in form and close to symmetrical, so these are equivalent to 1σ errors.

3.6 Synthetic data

To test the fitting procedure, we used synthetic BiSON-like ‘SolarFLAG’ data prepared as described by Chaplin et al. (2006) and Howe et al. (2015). The SolarFLAG data are constructed with |$\delta \Omega ={400\, {\mathrm{nHz}}}$| and visibility ratios of 0.54 for l = 2, m = 0: m = 2 and 0.38 for l = 3, m = 1: m = 3, the same values used for the prior distributions in our fitting. We generated 500 independent 11-yr realizations, each covering a solar cycle with activity variation based on Cycle 23. We then concatenated these realizations in sets of five, applying the BiSON gap pattern, to simulate the full BiSON history, yielding 100 realizations of a 45-yr spectrum. Each of the synthetic spectra was fitted in exactly the same way as the real BiSON one. By averaging the fitting results from many realizations, we can both check that the uncertainties in our results are appropriate and uncover any systematic bias in the results.

4 RESULTS

4.1 Synthetic data test with multiple realizations

We show in Fig. 3 the superimposed posterior distributions of the splitting parameter for fits to our 100 synthetic 45-yr spectra. We can see that there is some spread in the centroids of the distributions, increasing with frequency. In Figs 4(a) and (b), we show the mean value μ_σ of the 1σ error estimated from the distributions and the standard deviation σ_μ of the mean value, and in Fig. 4(c) we plot the ratio of μ_σ to σ_μ. We can see that the ratio values cluster below the y = 1 line, which suggests that the uncertainties from the fitting are underestimated, by about 30 per cent on average. To be more precise, the mean and standard errors of the μ_σ:σ_μ ratio are 0.64 ± 0.05 for l = 1, 0.74 ± 0.04 for l = 2, and 0.77 ± 0.03 for l = 3. This level of underestimation seems reasonable given that, because of the duty cycle of the data, the frequency bins will not be completely independent and hence equation (3) is not strictly appropriate. In a separate test where we fitted synthetic data with 100 per cent duty cycle, there was no such underestimate of the errors.

Fig. 5 shows the mean centroid μ_μ of the splitting estimates over our 100 realizations, plotted as a function of frequency with errors taken from σ_μ. Here, we can see that there is an upward bias, increasing with frequency, in the splitting estimates where the mode components are not fully resolved (above n = 14), but this is generally within the uncertainties and decreases with l. In the test with 100 per cent duty cycle, this bias is slightly reduced for l = 2 modes but still present.

4.2 Synthetic data test for sensitivity to the visibility scale factor

It has been shown (Chaplin et al. 2004, 2006) that choosing the wrong visibility ratios for the l > 1 components will result in a systematic bias of the splittings. In order to verify that our chosen method of using a constrained free parameter for the ratio will ameliorate this bias, we carried out three sets of tests in which we fitted a single realization of the SolarFLAG spectrum using a central value of the visibility factor prior that was scaled relative to our usual value by factors ranging from 0.7 to 1.3 in steps of 0.15. In the first set of tests, the σ of the prior distribution was set at the central value multiplied by 0.001; in the second set, it was 0.2 times the central value; and in the final set it was 0.4 times the central value. For each mode within each of the three sets, we then fitted a linear variation to the derived splitting values as a function of the prior central value for the visibility ratio, thus obtaining a value for the sensitivity of the splitting estimate to the central value of the visibility ratio factor. The results are shown in Fig. 6. This is similar to the test described by Chaplin et al. (2004), and the results in the first panel of Fig. 6 are very similar to theirs; above about 2500 μHz the sensitivity is such that a 100 per cent overestimate of the value of the visibility factor would shift the inferred splitting of the l = 2 and l = 3 modes upward by about 50 nHz. For the modes where the components are fully resolved, there is no sensitivity to the visibility factor. In the other two panels, we see that the sensitivity of the l = 2 and l = 3 splittings to the visibility factor is substantially reduced by using the wider prior distributions. Above about 3400 μHz, the change of prior has very little effect. Although at the 0.2 prior width the bias is not completely eliminated, we have sufficient confidence in our estimate of the true visibility ratios that we chose to use this prior rather than a wider one. In practice, when fitting the BiSON data, the differences between the splitting results for a 0.4 and 0.2 width of the visibility-factor prior are well within the uncertainties.

4.3 BiSON data

We now turn to fitting the 45-yr BiSON spectrum.

Fig. 7 shows selected mode-pair fits, plotted on a logarithmic scale. The smoothed input power spectrum is plotted, but it is almost obscured by the fitted model, giving a visual indication that the fit is working well. In the appendix, we show the corresponding ‘corner plots’ for these fits.

In Fig. 8, we show the prior and posterior distributions for the visibility ratios for the l = 2 and l = 3 modes as a function of frequency. For most of the modes, these are dominated by the prior; at low frequencies, this is because the signal-to-noise ratio of the inner components is low, while at high frequencies the components are not resolved. For a small number of modes between about 1700 and 2200 μHz, the visibility parameter appears to be somewhat constrained by the fit, and the posterior distributions suggest that the choice of prior was appropriate. We consider it reasonable to use the same prior for all frequencies. In recent resolved-Sun work (e.g. Larson & Schou 2015; Korzennik 2023), a term dependent on the radial order is included in the calculation of the leakage matrix (to which our visibility factors correspond) to account for the horizontal displacement of the modes, but according to Larson & Schou (2015) it is reasonable to neglect this for high-order modes.

Fig. 9 shows the prior and posterior distributions for the mode splitting, and Table 2 gives the median values and ±1σ uncertainties taken from the distributions. We note that the distributions for the two highest orders of l = 1 (which are unresolved) and the lowest order of l = 3 (which is lost in the background noise) are prior dominated and are not included in the table. There is a spurious second peak in the distribution for the l = 2, n = 7 mode at 1250 μHz, reflected in a very asymmetric uncertainty estimate for this mode, which is not included in the table. Again, this is a case where the signal-to-noise ratio of the mode is poor and the fit may be affected by a noise spike. Otherwise, the splittings appear to be well constrained and the posterior distributions are close to Gaussian. In Table 3, we show the splittings and uncertainties after adjusting for the bias and underestimation of uncertainties identified by using the synthetic data: the uncertainties have been scaled by the μ_σ:σ_μ ratio from the Monte Carlo experiment and the difference between the mean splitting estimates from the synthetic data and the true 400 nHz value has been subtracted from the splitting estimates.

We can estimate a mean splitting over several modes by combining samples drawn from the posterior distribution for each. After applying the adjustments used in Table 3, for 8 ≤ n ≤ 15 we obtain values of 399.3|$_{-9.6}^{+6.4}$| nHz for l = 1, 400.3|$_{-7.4}^{+4.2}$| nHz for l = 2, and 403.5|$_{-5.0}^{+5.4}$| nHz for l = 3, while for 16 ≤ n ≤ 23 the values are 402.6|$_{-36.0}^{+30.1}$| nHz for l = 1, 411.6|$_{-12.3}^{+11.5}$| nHz for l = 2, and 407.6|$_{-12.0}^{+9.5}$| nHz for l = 3. All of these values are consistent with the 400 nHz value used to construct the SolarFLAG data.

Finally, in Fig. 10 we show the splitting values as a function of frequency in error bar form, overlaid on the results from Fig. 5. The trends at frequencies above 2000 μHz appear generally consistent, with a few notable outliers. The low observed splitting value for the l = 1, n = 17 mode at 2559 μHz is particularly striking and seems to be robust, as it occurs even in fitting with different choices for the priors, but it is most likely a result of the stochastic nature of the mode excitation rather than an intrinsic feature of the underlying rotation profile. The l = 1 splittings in general show some scatter around the synthetic average values, but there is no obvious systematic trend.

5 DISCUSSION AND CONCLUSIONS

We have derived rotational splitting parameters for modes of l = 1, 2, and 3 from a BiSON spectrum spanning 16 425 d from 1976 December 31 to 2021 December 20. Our algorithm assumes a single splitting parameter per mode, with the visibility ratio for the inner components of l = 2 and l = 3 modes as a variable parameter. The estimated value for the visibility ratio is dominated by its prior except for the few orders where the inner components are well resolved and have good signal-to-noise ratios. We present the posterior distributions from our MCMC fits as well as the centroid estimates and their uncertainties. It should be noted that the l = 3 splitting values do not include a differential rotation term and are not strictly equivalent to the first-order term of a polynomial expansion of frequency as a function of m.

We have presented the results for modes where the splittings are not resolved, for the sake of completeness. However, the Monte Carlo experiments we have performed with simulated data suggest that the splitting estimates for these modes (n > 14) are systematically biased upwards and should not be used for inversions, at least not without appropriate correction. This bias appears in l = 1 modes as well as the higher degree modes, so we believe it is intrinsic to the problem of fitting unresolved peaks and not related to the modelling of the visibility function of the modes nor to the neglect of various higher order and yearly effects that are included neither in our fits nor in the synthetic data. The upward trend in splitting with frequency in the fits to BiSON data appears to be consistent with that from SolarFLAG data where the splitting is constant at 400 nHz, and therefore these results do not suggest that there is any real upward trend in the splitting with frequency.

Our Monte Carlo experiments also reveal that the uncertainties from the fitting are underestimated by about 30 per cent. We believe this is due to the reduction in the effective resolution of the power spectrum because the overall duty cycle is relatively low due to sparse observations early in the time series. The uncertainties should be corrected accordingly when using the splittings together with data from other sources, and we have provided a table of data and uncertainties adjusted based on the results of our Monte Carlo experiments.

We note that the BiSON observations have a variable duty cycle, being sparser in the early years, and hence the measurements from the 45-yr spectrum are not a uniform average over those years but more heavily weighted towards the post-1992 epoch.

In future work, we hope to compare splitting estimates from a shorter subset of the BiSON data with those from space-based Sun-as-a-star instruments such as GOLF and resolved-Sun instruments such as the Global Oscillations Network Group, the Michelson Doppler Imager, and the Helioseismic and Magnetic Imager. We also plan to refine the fitting algorithm to provide definitive values for the mode frequencies.

SOFTWARE

Additional software used in this work, which have not explicitly been mentioned above, are listed below:

• python (Van Rossum & Drake Jr 1995)

• astropy (Astropy Collaboration 2013, 2018)

• matplotlib (Hunter 2007)

• numpy (Harris et al. 2020)

• scipy (Virtanen et al. 2020)

• corner (Foreman-Mackey 2016)

ACKNOWLEDGEMENTS

We acknowledge the support of the UK Science and Technology Facilities Council (STFC) through grant ST/V000500/1. The computations described in this paper were performed using the University of Birmingham’s BlueBEAR HPC service, which provides a High Performance Computing service to the University’s research community. See http://www.birmingham.ac.uk/bear for more details. This research made use of NASA’s Astrophysics Data System. We would like to thank all those who have been associated with BiSON over the years. We thank Dr Annelies Mortier for useful suggestions. We thank the anonymous referee for helpful comments.

DATA AVAILABILITY

The BiSON time series analysed here is available at http://bison.ph.bham.ac.uk/opendata. The SolarFLAG time series is available on reasonable request to the authors.

References

APPENDIX A: CORNER PLOTS

Figs A1–A4 show the ‘corner plots’ for sample fits to the BiSON spectrum. These are the fits that are shown in Fig. 7. We can see that there some correlations among the parameters, especially for the higher frequency case. There tends to be a positive correlation between the visibility scale parameter and the l = 2 or 3 splitting, which is not unexpected. The l = 3 amplitude is somewhat negatively correlated with the visibility scale parameter, and hence it is also anticorrelated with the l = 3 splitting. None of the correlations we see are of particular concern, given that our final results are taken from the marginalized posterior distributions.