Determination of the rotation of the solar core requires very accurate data on splittings for the low-degree modes which penetrate to the core, as well as for modes of higher degree to suppress the contributions from the rest of the Sun to the splittings of the low-degree modes. Here we combine low-degree data based on 32 months of observations with the BiSON network and data from the LOWL instrument. The data are analysed with a technique that specifically aims at obtaining an inference of rotation that is localized to the core. Our analysis provides what we believe is the most stringent constraint to date on the rotation of the deep solar interior.
The angular velocity of the solar core has been the subject of considerable debate. Early helioseismic investigations (Claverie et al. 1981; Duvall & Harvey 1984; Duvall et al. 1984) provided some indications for rotation of the core at a speed somewhat exceeding the surface rotation rate. More recently, measurements by the IRIS network (Lazrek et al. 1996) have indicated that the core rotates at, or perhaps slightly above, the surface rate. Somewhat similar results have been obtained through inversion of combined data from the GOLF and MDI instruments on the SOHO spacecraft (Corbard et al. 1998). On the other hand, inferences from BiSON network data (Elsworth et al. 1995) and from LOWL data (Tomczyk, Schou & Thompson 1995a) have shown some evidence for core rotation below the mean rate of the solar envelope. The present state of disagreement is discussed further by Eff-Darwich & Korzennik (1998).
It is observed that the surface rotation rates of solar-like stars generally decrease with the inferred age of the stars. To account for this, it is assumed that the stars arrive on the main sequence in a state of rapid rotation and subsequently lose angular momentum through magnetically controlled stellar winds. This slows down the outer convection zone, in which the magnetic fields are anchored and where angular momentum is redistributed very efficiently by convection. The deceleration of the radiative interior depends on the transport of angular momentum between this region and the convection zone. Models that rely on transport through turbulence induced by the rotational shear (e.g. Chaboyer, Demarque & Pinsonneault 1995) suggest that the transport is relatively inefficient, possibly leaving behind a fairly rapidly rotating core. On the other hand, even a weak large-scale magnetic field would be sufficient to couple very efficiently the interior and the convection zone, leading to essentially solid-body rotation (Rosner & Weiss 1985; Mestel & Weiss 1987; cf. also Charbonneau & MacGregor 1993). There have also been speculations that the rotation of the core may be variable, perhaps on the time-scale of the solar cycle (e.g. Gough 1985). Evidently, constraints on the rotation of the core are of great value to help distinguish between these diverse theoretical models.
The observed global p-mode oscillations of the Sun are labelled with three quantum numbers: the radial order n, degree l and azimuthal order m. Such a mode has horizontal structure described by the spherical harmonic The leading-order effect of the Sun's rotation is to split the frequencies of modes that have the same n and l but different m: to this order, the difference between the angular frequency ωnlm of the (n, l, m) mode and the frequency ωnl0 of the corresponding m=0 mode is m times a mode-weighted measure of the Sun's interior rotation rate. We define the splitting to beequation (2) need only be taken from zero to π/2.
Inferences about the core rotation are complicated by the fact that the dominant contributions to the frequency splittings come from the outer parts of the Sun, where the kernels Knlm are largest. Thus, in the analysis, data for the low-degree modes which penetrate to the core must be combined with data with higher-degree modes, to eliminate the contributions from the rest of the Sun.
Except at the lowest degrees l, we do not determine observationally all the splittings Δωnlm individually. Rather, the data analysis is made by representing the m variation of the splitting (for fixed n and l) with a low-order polynomial in m:Schou, Christensen-Dalsgaard & Thompson (1994) are used, and kmax is taken to be 3, so that we have so-called a coefficients a1, a3 and a5. For the BiSON data we have individual sectoral cyclic frequency splittings Δνnll≡(2π)−1Δωnll.
Here there are M data, and the index i numbers the data di. These may be either individual splittings m−1Δωnlm, or expansion coefficients 2πa2k−1(n,l) resulting from expansion (3) of Δωnlm. Obviously the kernels Ki depend on the type of data considered.
In practice, of course, equation (4) does not hold exactly because the data contain errors. We assume that the errors in the di are independent, and normally distributed with zero means and with standard deviations σi that are estimated as part of the analysis of the observations.equation (4). For the scaled data, the error covariance matrix is the identity matrix.
Quite a number of techniques have been developed and applied to the problem of inverting the integral constraints (4) to infer the solar rotation Ω(r, θ) (see e.g. Schou et al. 1998). Here we use three methods. The first, regularized least-squares (RLS), has the virtues that it is quite fast and that it attempts to find a solution that fits the data. The other two methods are both of the class of optimally localized averages (OLA) methods, which are based upon an original idea by Backus & Gilbert (1968). The two variants we use are described further in Section 3. Although they are computationally more expensive and slower, the OLA methods have the advantage that they make a better job of obtaining localized information about the rotation in the interior. In this paper we are particularly concerned to obtain a localized estimate of the rotation rate in the solar core, and so the use of the OLA methods is particularly appropriate.
2 Observations and data analysis
The BiSON rotational splittings come from the analysis of a 32-month power spectrum, generated from velocity residuals collected between 1994 May 16 and 1997 Jan 10, i.e. at or near the solar activity cycle 22/23 minimum. Adjacent l=2 and 0, and l=3 and 1 modes were fitted in pairs with maximum likelihood estimators (χ2, two degrees-of-freedom statistics assumed) as rotationally split, Lorentzian multiplets with an associated flat background offset. The first temporal sidebands were also included in the model. Powell's multidimensional direction-set minimization algorithm (Press et al. 1986) was used to perform the fitting. Of relevance here, a single parameter describing a symmetric rotational splitting pattern was fitted to each non-radial mode multiplet. Further, the associated uncertainties are those derived from the inverse of the Hessian matrix of each fit. As the sectoral (i.e. m=±l) mode components are the strongest at all l in the ‘unresolved’ BiSON observations, the fitted values are essentially sectoral splittings.
The values of the fitted splittings for l=1 and 2 are given in Table 1. The quoted uncertainties are generally an increasing function of frequency, which is a consequence of the greater linewidth of higher frequency modes. As the very small uncertainties on the lowest frequency splittings increase the importance of these data to the inversions, it behoves us to look with special care at those fits. As examples we have carefully reviewed the fitting of the l=1, n=9,10 modes. The resonant peaks of these modes appear with good signal-to-noise ratio in the BiSON power spectrum. The maximum-likelihood fits to the peaks also appear sound, i.e. the residuals generated by dividing the raw spectrum through by the fit show no evidence of any structure that might be indicative of poor convergence. In order to investigate further the reliability of these splittings, we performed a series of Monte Carlo simulations. We generated and then fitted artificial p-mode-like data that were designed to mimic the observed characteristics of the low-l pairs that include the l=1, n=9 and l=1, n=10 modes (see Chaplin et al. 1997). We repeated this procedure many times for each low-l pair in order to build up fitted parameter distributions from the independent realizations processed at each l and n. The results indicate that the ‘limit’ estimates of the splittings (i.e. the average of the fitted artificial splittings at each l and n) are unbiased. We also studied the distribution of the fitted formal uncertainties returned by each fit. These would appear to indicate that the errors on the real BiSON splittings at n=9 and 10 are perhaps near the ‘low’ end of the distribution one might expect, given the characteristics of the real modes, and the quality and length of the data set. Further, the ratio of the external and internal errors in the distribution of artificial splittings would appear to suggest that — on average — any formal error bars may underestimate the actual precision in the splittings by roughly 15 to 20 per cent or so. Nevertheless, we do not see this as reason to reject these two splittings in the analysis.
At frequencies above 3500 μHz, the decreasing modal lifetimes give rise to substantial line blending — this makes the task of obtaining reliable estimates of the parameters of a low-l pair increasingly difficult. This is particularly true for the rotational splittings, because at such frequencies their magnitude is only a small fraction of the modal line widths. We therefore include splittings for analysis up to n=23 (at both l=1 and 2) but no higher. Several studies have indicated a substantial bias in fitted rotational splittings at high n, with a clear tendency to overestimate their true magnitude (Appourchaux et al. 1997; Chang, Gough & Sekii 1997; Chaplin et al. 1998b). Chaplin et al. also demonstrated that the formal uncertainties in the splittings returned by ‘disc-integrated’ mode-fitting procedures are anti-correlated with the magnitude of the extracted splittings, an effect that increases at higher n. Further, they showed that changes in the correlation as a function of frequency mimic the variation of the mode linewidths, i.e. the larger anti-correlations displayed at high n result from more severe mode blending in the power spectrum. Estimates for l=1 exhibit the largest bias and correlation, owing to the closer separation in frequency of the sectoral mode components (to which full-disc observations are most sensitive). In light of the above, we therefore add a note of caution regarding the high-n splittings included in Table 1.
For l3 we used splittings from the LOWL instrument (Tomczyk et al. 1995b; LOWL is an abbreviation for low-degree, with degree denoted by L) which has been operating continuously at Mauna Loa, Hawaii since 1994 February 26. The LOWL employs a magneto-optical filter to measure the Doppler shift of the solar surface in the potassium line at 770 nm, with a spatial resolution of approximately 20 arcsec, sufficient to observe modes with 0l99.
For the results presented here we have used two years of observations covering the period 1994 February 26 to 1996 February 25. For computational convenience the data were analysed as two separate 1-yr long time series and averaged. The fitting for each year was performed as described by Schou (1992). Owing to the rapid increase in linewidth with frequency, only modes with ν3500 μHz have been used. The analysis of the observations to produce splitting data is described in detail elsewhere (Schou & Tomczyk, in preparation). To stabilize the fits, individual splittings Δωnlm were not sought: rather, the fits were performed in terms of a coefficients, described above (see equation 3).
2.3 Compatibility of the data sets
The combination of the BiSON and LOWL data sets should produce reasonably homogeneous results, given that both systems employ the same basic physical principles, i.e. using an atomic reference standard to observe the oscillations via Doppler velocity shifts of the same Fraunhofer line. To investigate the compatibility further, we examine the splittings from the two instruments. As stated above, the BiSON measurement is essentially a sectoral mode splitting. Fig. 1(a) illustrates the BiSON l=1 and l=2 splittings together with the LOWL data combinations (a1+a3+a5) equivalent to sectoral splittings Δωnll for modes with l=3–5, plotted against ν/(l+1/2) (closely related to turning-point position; see Section 3.1 below). There is no obvious inconsistency between the BiSON and LOWL data. This also suggests that the lack of exact contemporaneity of these data sets is not a significant concern in the inversion of the frequency splittings.
3 Inversion techniques
3.1 Properties of kernels
Our ability to infer the rotation of the deep solar interior depends on the availability of data for deeply penetrating modes. In general, the behaviour of the rotational kernels reflects the fact that the modes have a lower turning point at rrt, such that c(rt)/rt≃ω/(l+1/2), c being the adiabatic sound speed. Thus modes of low degree and/or high frequency penetrate most deeply. For low-degree modes, the behaviour at radii smaller than the turning-point radius is determined by the expansion of the eigenfunctions around the regular singular point at r=0. This shows that, near the centre, the kernels behave like Ki∝r2l; however, for l=1 the leading order terms cancel, so that Ki∝r4. Hence splittings for l=1 are less sensitive to the e core rotation than might otherwise have been expected. These properties are illustrated in the kernels shown in Fig. 2, for an angular velocity Ω=Ω(r) that does not depend on latitude. Further properties of the core contributions to the splitting are discussed in Section 4, below.
3.2 General aspects of inversion
It follows from equations (4) and (6) that
The errors in two solution points (be they e.g. at two different locations, or from two different solution methods, or from one method with two different choices of parameters in an inversion method) are correlated in general, because the solutions use the same data. If (1) and (2) are two solution points with corresponding inversion coefficients and then the correlation between (1) and (2) is given by
3.3 Regularized least-squares inversion
In the two-dimensional regularized least-squares (RLS) inversion method, we minimize a sum of terms, one of which is essentially the χ2 of the fit to the data; the remainder of the expression serves to regularize the solution, by penalizing rapid variations that might otherwise be introduced in an attempt to fit the data or their errors. Specifically, we minimizeSchou et al. 1994 for full details). The relative importance of the regularization terms in equation (11) depends upon the choice of values for μr and μ: larger values reduce the effect of data errors and hence give a less rapidly varying solution, but if they are chosen to be too large then the solution will not give a good fit to the data because there is no reason a priori why the true rotation rate should have very small second derivatives with respect to r and θ.
To minimize expression (11) numerically, we discretize the rotation rate as a piecewise bilinear function on a rectangular grid in (r, θ) space, namely(11) then becomes a problem in linear algebra: the coefficients and hence, by equation (12), the solution itself depend linearly on the data, so that the solution is of the form (6) with a set of inversion coefficients i for each point on the solution. Thus averaging kernels (equation 7) exist which relate the solution to the true rotation rate as in equation (8); and the error properties of the solution are as stated in equations (9) and (10).
3.4 Ola techniques
We use two optimally localized averages (OLA) techniques, the so-called subtractive (SOLA) and multiplicative (MOLA), cf. Pijpers & Thompson (1992). In the case of SOLA, the inversion coefficients are found so as to make the averaging kernel an approximation to a given target T (r0, θ0, r, θ), by minimizing
It may be noted that the matrix multiplying the vector of inversion coefficients on the left-hand side of equation (14) is independent of (r0, θ0); thus it can be inverted (or, more efficiently, factorized) once and for all.Equation (18) leads to the following set of linear equations for i(r0, θ0): equation (15); here
Note that the matrix Jij(r0, θ0) must be set up, and equation (19) solved, separately for each target location (r0, θ0). Thus the MOLA technique is generally much more demanding on computational resources than is the SOLA technique. However, both have been implemented much more efficiently than hitherto. In the present implementation, the factorization of equations (14) (SOLA) and (19) (MOLA) is computed using an iterative method based on the Lanczos bidiagonalization algorithm (e.g. Golub & Van Loan 1989). This makes it possible to take advantage of the algebraic structure of i(r, θ) in the calculation, and thereby reduces the computing time by almost two orders of magnitude — see Larsen & Hansen (1997) and Larsen (1998).
The SOLA target functions were taken to be Gaussians, symmetrized around the equator (cf. Schou et al. 1994). In accordance with the resolution properties of the kernels, the target widths in radius were chosen to be proportional to the sound speed at the target location (e.g. Thompson 1993), whereas the target widths in latitude were taken to be constant. The choice of the weight function J in the MOLA method was motivated by the desire to concentrate the averaging kernel as near as possible to the core. We have used
As no localization in latitude was attempted, Δ was chosen to be 3R; the localization in radius was controlled by varying Δr.
Having found inversion coefficients by either the SOLA or the MOLA method, the inferred value (r0, θ0) is then obtained from equation (6).
4.1 Properties of the inferred rotation rate
Our inferences regarding the solar internal rotation are illustrated in Fig. 3. Panel (a) shows radial cuts through the RLS solution at latitudes 0° (equator), 30° and 60°. Panel (b) shows SOLA inversion results at the same three target latitudes. Also shown in panel (b) are four MOLA results targeted in the core, achieving different levels of confinement to the deep interior. In the convection zone (r≳0.7R) we see surface-like latitudinally dependent rotation, as have previous investigations (e.g. Thompson et al. 1996). Beneath the base of the convection zone the inversions reveal no significant latitudinal dependence in the rotation rate, although both RLS and SOLA indicate that the equatorial rate is higher than those at other latitudes at r≃0.55R, and the latitude 60° rate is higher at r≃0.4R. The SOLA solution in particular indicates a local maximum in the rate at all latitudes at r≃0.4R, though this is barely a 1σ result. In the core (r≲0.3R), where we believe the MOLA technique has produced better localized results than ever before, there is an apparent trend to slower rotation: we discuss the inferences in the core in more detail below. However, we already note that our results are consistent with the rotation profile of the outer core being a flat continuation of the rotation rate in the rest of the radiative interior.
The confluence of the RLS solution at all latitudes below r=0.25R is simply an indication that we have no differentiation in latitude in the core: regardless of the nominal latitude, the averaging kernel looks the same. One might call this an artefact of the inversion, but principally it is a consequence of having so few m values at low degree, which severely limits the latitudinal resolution. The linear trend with radius in this region is an artefact: the second-derivative smoothing causes such a linear extrapolation in regions where there is insufficient data to influence the solution. Reducing the regularization would weaken the linear trend, but at the expense of increasing substantially the error bars on the solution.
The apparent discrepancy between SOLA and RLS in the convection zone is caused by the different latitudinal smoothing created by the averaging kernels from the two methods. The SOLA averaging kernels are essentially everywhere positive, and so they to some extent average out the latitude variation in the convection zone. This pulls down the inferred rotation rate at the equator, where the rotation is fastest, and pushes it up at high latitudes where the rotation is slower. The features in the inversion profiles around r=0.9R are likely artefacts of the LOWL data for modes with turning points around 0.9R, where diurnal sidelobes in the Fourier spectra cause systematic errors in the fitted mode parameters.
Fig. 4 illustrates some averaging kernels. Fig. 4(a) shows an RLS averaging kernel at the equator and r=0.17R. The inversion at that point has a formal error (cf. Fig. 3) of about 21 nHz. Fig. 4(b) shows the averaging kernel for the deepest of the three filled-circle MOLA inversions in Fig. 3(b), which corresponds to a similar formal error (20 nHz). The advantage of the MOLA kernel is readily apparent: the kernel has none of the near-surface structure, nor the ringing (positive and negative structure) outside the core, that is possessed by the RLS kernel. The MOLA kernel therefore gives a much better localized measure of the core rotation.
The RLS inversion explicitly attempts to fit the data, and large residuals (i.e. data minus fit) would be an indication of inconsistent data sets or underestimated data uncertainties. We note from Fig. 5 that the RLS residuals look entirely consistent with the quoted uncertainties on the data, and reveal no inconsistency between the BiSON and LOWL data sets (cf. also Fig. 1).
4.2 Relative weights given to different data
The MOLA inversion coefficients ci (Fig. 6) show the weight assigned by the inversion to each datum di (cf. equation 6). Panel (a) corresponds to the deepest of the filled circles in Fig. 3, but the relative weights are similar for the other filled-circle cases. Note in particular that the two very precise low-frequency l=1 splittings, of radial order n=9 and 10 (two open triangles in Fig. 6a) contribute strongly in these MOLA inversions. Indeed, we have found that the slight jump in the inferred SOLA/MOLA rotation rate (Fig. 3) would be eliminated if these two modes were not included. This is fairly evident from the fact that their splitting values are about 20 nHz greater than the ‘ambient’ 430-nHz rotation rate of much of the radiative interior, coupled with the size of their inversion coefficients. It follows, however, from the Monte Carlo simulations described in Section 2.1 that there is no apparent reason to throw out those two splittings. We also note (Fig. 5) that the residuals from an RLS fit look entirely consistent with the quoted uncertainties on the data, and provide no argument for weeding out these or any others of our low-degree data.
When the MOLA inversion is pushed even harder (deepest open circle in Fig. 3), the inversion coefficients (Fig. 6b) show a very interesting qualitative difference. We see now that the precisely determined low-frequency l=2 modes have most weight. (Furthermore, by comparing the largest l=1,2 coefficients in Fig. 6b with the values in Table 1, one sees that there is a very strong correlation between the precision of the splitting measured and the weight assigned in the inversions.) Fig. 2 shows that the kernels of the low-frequency modes have their deepest maximum at or outside r=0.2R. Beneath this, the l=1 kernel quickly drops to zero (cf. the discussion in Section 3.1). Thus these modes are useful for forming an averaging kernel around r=0.2R, but not significantly deeper than that. The modes of degree l=2 and higher, although at comparable frequency they have turning points further out from the core, also have quite significant ‘tails’ extending quite far beyond their lower turning points. If the errors on the corresponding splittings are sufficiently small, these tails can be used by the inversion to construct an averaging kernel at greater depth than can be attained by the l=1 modes. To understand the importance, particularly of the l=2 data, to the MOLA construction of a kernel at around r=0.15R, we look at the partial integrals of mode kernels from r=0 to r=0.15R (upper panel of Fig. 7). At most frequencies it is the l=1 kernels that have the largest partial integrals. However, at the lowest frequencies the tails of the l=2 kernels cause them to have larger partial integrals than the l=1 kernels, even though the nominal turning points of the l=1 modes are deeper than those of the l=2 modes. The MOLA inversion constructs an averaging kernel from these tails; but in doing so it must also take account of the data uncertainties (cf. equation 5). Thus in the lower panel of Fig. 7 we illustrate the partial integrals divided by the (scaled) data uncertainty σ— see the Fig. caption. The l=2 data have smaller errors than the l=1 data, and so this reinforces the importance of the l=2 data over the l=1 data. There is a close correlation between the modes with the largest values of σ−1∫Knl dr and those with largest coefficients (Fig. 6b).
4.3 Does the Sun have a slowly rotating core?
Finally, we return to the evidence for a slowly rotating core presented by the downward trend of the four MOLA points in Fig. 3. One should beware that the errors on the four inversion values are correlated, because they are of course constructed from the same data. In fact, the correlation is strong: using as reference the least deeply penetrating of the large filled circles in Fig. 3, and in order of increasing depth, the correlations are 0.958, 0.850, and 0.409 respectively. To investigate this further, we have made inversions of artificial data constructed using a solar-like rotation in the convection zone and a rigidly rotating radiative interior. To these artificial splittings we have added 20 realizations of random errors, with independent Gaussian distributions corresponding to the estimated uncertainties in the real solar data. Fig. 8 shows the results of inverting (with identical SOLA and MOLA procedures to those used on the real data) the error-free artificial splittings and two of the 20 realizations with added noise. It can be seen that one of our 20 realizations produces a strong downward trend in the core very much like what we find for the Sun, even though the actual underlying rotation is uniform. Another case generally gives a similarly strong upward trend. This is not surprising: we would expect a couple of cases as extreme as one-and-a-half standard deviations away from the mean in 20 realizations, and this is roughly what we have. Because the first three MOLA points are so strongly correlated (see above) they all show essentially the same extreme behaviour at once, leading the eye to see an organized trend. On the other hand, the innermost MOLA point is not so strongly correlated, as noted above, and follows the trend of the other MOLA points in only one of the two cases illustrated. The downturn in the core inferred from real solar data (Fig. 3) is actually less extreme, roughly a 1σ result, so is even more likely to arise by chance as a result of random noise in the data, even if the core rotation were flat. We do not, therefore, regard the downward trend in Fig. 3 as significant.
Using a modification of the MOLA method which penalizes the averaging kernels for having any structure far from the centre of the Sun, we have succeeded in making measures of the Sun's rotation that are strongly localized in the solar core. We believe that these results are the best to date in terms of clean localization in the core. Formally, our most likely solution is that the core is rotating more slowly than the rest of the radiative interior. However, our solution in r=0.15–0.2R is consistent with a flat rotation curve, because the slower rotation is no more than a 1σ result. We obtained a qualitatively similar downturn in core rotation in our preliminary study (Chaplin et al. 1998a), but with the present revised estimates of splittings and their associated errors, the downturn appears to be less significant. The localization we have now achieved is much superior to that attained in earlier studies with BiSON and LOWL data separately (Elsworth et al. 1995; Tomczyk et al. 1995a), which also indicated a slower rotation in the solar core. Charbonneau et al. (1998) inverted the LOWL data alone and produced a result consistent with a flat rotation profile in the core: although they used a genetic algorithm rather than a linear inversion, we believe that our localization is also superior to that implicitly achieved by them. We emphasize that our new measures of the rotation are truly localized at the radii we claim for them to a very great extent (shown precisely by the averaging kernels in Fig. 4): by contrast, earlier quoted values for the rotation at e.g. r=0.1R were less well-localized extrapolations.
To obtain localized measures of the rotation any deeper than those presented in this paper, it would be highly desirable to have splittings for g modes. (Data for g modes are highly desirable for structural inversions also, but at least in that case we have data about deeper regions from p modes of degree l=0.) However, some further progress in constructing localized measures of the core rotation might yet be possible with p modes, thanks to the very high precision with which splittings of the low-frequency modes can be determined. In particular, for l=2 modes this compensates for the lower turning points of a mode moving further from the centre as frequency diminishes (Fig. 7). We are studying the potential to push the rotation inversions even deeper with future observations of these modes (Christensen-Dalsgaard et al., in preparation).
We are grateful to T. Corbard and an anonymous referee for comments on earlier versions of the text. This work was supported in part by the UK Particle Physics and Astronomy Research Council, by the Danish National Research Foundation through its establishment of the Theoretical Astrophysics Center, and by the US National Science Foundation through base funding of HAO/NCAR. WJC acknowledges the support of an ESA Internal Fellowship. We thank all those who are — or have been — associated with the BiSON global network and all our hosts at the network sites.