© ESO, 2016

1. Introduction

One of the features of many solar dynamo models is the rotational shear layer below the convection zone, the tachocline (Spiegel & Zahn 1992), where the magnetic field is thought to be generated in some models (e.g. Charbonneau 2010). Despite knowing the spatially-resolved internal rotation profile of the Sun (Schou et al. 1998), it is poorly constrained in other Sun-like stars. Asteroseismology is the only tool available to measure the rotation of stellar interiors.

In the power spectrum, modes of oscillation have strong peaks at their oscillation frequency. When the rotation rate is slow the oscillation is perturbed so that the mode becomes a multiplet in frequency, ω_nlm = ω_nℓ + mδω_nℓ, ignoring any latitudinal differential rotation, where n is the radial order, ℓ is the harmonic degree, and m is the azimuthal order of the mode, described by a spherical harmonic. These azimuthal components are called splittings.

From the splittings, an inversion technique can be used to infer the internal rotation profile. This has been performed for stars with large rotation rates (and easily measured splittings) and modes sensitive to different depths, tightly constraining the core rotation rate (e.g. Deheuvels et al. 2012, 2014). However, in Sun-like stars, inferring the radial differential rotation (RDR, hereafter) is more challenging. Lund et al. (2014) found that it is unlikely that asteroseismology of Sun-like stars will result in reliable inferences of the RDR profile, such as can be performed for red giants. This is predominantly due to the degeneracy of the information contained in the observable modes: they are sensitive to roughly the same regions of the stellar interior (see Figs. 3 and 4 in Lund et al. 2014). They posit that, at best, it may be possible to obtain the sign of the RDR gradient. Despite the inherent challenges, rotational splitting as a function of frequency has been measured in Sun-like stars by Nielsen et al. (2014, 2015) using the short-cadence long duration observations from Kepler (Borucki et al. 2010).

Schunker et al. (2016, hereafter Paper I) conclude that it would not be possible to infer a radially resolved rotation profile for a single Sun-like star using regularised inversion techniques. As an alternative, they demonstrate that direct functional fitting of a step-function can constrain the sign of the RDR gradient better and can even be used to retrieve the surface rotation rate of a subgiant star.

Paper I also explores the sensitivity of linear inversions to the stellar models used in the inversion. They conclude that the inversions are insensitive to the uncertainties in the stellar model, compared to the level of noise in the splittings. They therefore propose that it may be possible to exploit this insensitivity to perform ensemble asteroseismic inversions across a broader range of stellar types.

Ensemble asteroseismic inversions have the advantage that the uncertainties will be reduced with the extra information from more stars, from which an ensemble inversion can be solved to infer an average RDR profile. The upcoming PLATO 2.0 mission (Rauer et al. 2014) will observe tens of thousands of solar-like stars, offering us the opportunity to exploit this ensemble fitting method to constrain the RDR gradient in Sun-like stars.

In this paper we explore the possibility of using ensemble asteroseismic fits of Sun-like stars to infer the average RDR. We first solve the forward problem to compute the rotational splittings for a broad range of stellar models with a synthetic rotation profile (Sect. 2). We describe the fitting function and the ensemble fit method (Sect. 3). We show the results for independent fits (Sect. 4) and apply the ensemble fit to all 15 stellar models in Sect. 5. We then extend our results to estimate how well we can constrain the size of the step function along the main sequence using subsets of the Sun-like stars that are expected to be observed by PLATO (Sect. 6). We present our conclusions in Sect. 7.

2. The forward problem: modelling splittings

2.1. Stellar models

We compute a set of 15 stellar models of main sequence solar-like oscillators that covers F to K-type stars, with a range of metallicities and ages using the Modules for Experiments in Stellar Astrophysics code (MESA¹ revision 6022, Paxton et al. 2013). Opacities are taken from OPAL (Iglesias & Rogers 1996) and Ferguson et al. (2005) at high and low temperatures, respectively. The equation of state tables are based on the 2005 version of the OPAL tables (Rogers & Nayfonov 2002), and nuclear reaction rates are drawn from the NACRE compilation (Angulo et al. 1999) or Caughlan & Fowler (1988), if not available in the NACRE tables. We used newer specific rates used for the reactions (Imbriani et al. 2005) and (Kunz et al. 2002). Convection was formulated using standard mixing-length theory (Böhm-Vitense 1958) with mixing-length parameter α = 2. The solar metallicity and mixture was that of Grevesse & Sauval (1998) and atomic diffusion was included, using the prescription of Thoul et al. (1994).

Models were initialised on the pre-main-sequence with a central temperature T_c = 9 × 10⁵ K. We computed tracks for masses 0.8, 1.1 and 1.4 M_⊙ and initial metallicities [ Fe / H ] = −0.5 and 0.5. Along each track, we recorded the models with central hydrogen abundances X_c = 0.60, 0.35, and 0.10 to represent stars towards the beginnings, middles, and ends of their main-sequence lives. We excluded models from the low-metallicity 1.4 M_⊙ track because all had effective temperatures above 7000 K and would therefore either not be solar-like oscillators, or their oscillation frequencies would be difficult to identify and measure. The parameters of the models are presented in Table 1. The rotational kernels for each mode of oscillation in the star can be computed from these models.

2.2. Synthetic rotation and mode splittings

The Sun’s convection zone shows radial and latitudinal differential rotation. The analysis of the acoustic modes which have low sensitivity in the deep interior suggests solid body rotation in the core (Schou et al. 1998). Thus, we chose a step-function profile as the simplest representation of the expected rotation profile of a star. Specifically, we chose (1)where ΔΩ_c = −350 nHz, R is the stellar radius, and r_BCZ is the radius of the base of the convection zone for the particular stellar model (see Table 1).

We ignore any latitudinal differential rotation. The rotational splittings are a linear function of the angular rotation rate, , where K_nℓ are the rotational kernels of the star. For δω_n2 we compute the error-weighted mean of the splittings for the | m | = 1,2 components, where the errors are the uncertainties modelled in Sect. 2.3.

2.3. Noise model for the rotational splittings

Global rotational splittings have been unambiguously measured (e.g. Gizon et al. 2013) for Sun-like stars. The only attempt to measure variations in the splittings for Sun-like solar-like oscillators was done as a function of mode-set (ℓ = 2,0,1, Nielsen et al. 2014), and not for independent modes. Therefore, we need to model the measured uncertainties on the splittings. We scale the uncertainties on the splittings from a function of the uncertainties on the frequencies, which is similar to Paper I.

We begin with a model for the uncertainties on frequencies from Libbrecht (1992, Eq. (2)): (2)where Γ is the FWHM of the mode linewidth, T is the observing time, which we set to the expected three years of future asteroseismology observing campaigns, and f(β) = (1 + β)^{1 / 2} [ (1 + β)^{1 / 2} + β^{1 / 2} ] ³, where β is the inverse signal-to-noise ratio B(ν) /H_ℓm. Here we are specifically employing ν to indicate a function of frequency, to distinguish against our use of σ(ω_ℓm), which indicates the uncertainty in the mode frequency ω_ℓm. Here, B(ν) is the background noise spectrum and H_ℓm is the mode height. We determine Γ relative to the radial order of the mode closest to the frequency of maximum power, ν_max. H_ℓm and B(ν) are scaled from the solar values relative to ν_max.

We first model the linewidths of the Sun by fitting the values for each mode in Table 2.4 of Stahn (2011) with a third order polynomial as a function of radial order relative to the mode closest to ν_max to get Γ(n). The function is then scaled using the relation (Chaplin et al. 2009) for each star.

We then model the envelope of mode heights in the Sun using the function (3)(Stahn 2011) relative to ν_max, and scaled by the surface gravity (Chaplin et al. 2009). The parameters of the amplitude , and ν_c, σ₁, σ₂ of the asymmetric Lorentzian, P, are given in Stahn (Table 2.3, 2011). Here, we choose the inclination of the stellar rotation axis to be i = 90° in the mode visibility term, ℰ_ℓm(i) since, statistically, this inclination is most likely. We discuss the statistical effect of random inclinations later in this paper (Sect. 6). The distribution of multiplet power as a function of stellar inclination for ℓ = 1,2 modes can be found in Gizon & Solanki (2004).

The background term, B(ν), is modelled using two Harvey laws (Harvey 1985) with solar values from Table 2.1 of Stahn (2011), which are scaled by the frequency of maximum power, B/B_⊙ ∝ ν_max/ν_{⊙ max}. We note that including the background term increases the uncertainties in the splittings by up to 20% for the stars with higher T_eff.

Lastly, we compute the uncertainties for the frequency splittings slightly differently to Paper I, which assumed that all modes were equally visible. The uncertainty on the splittings is then σ_m(δω) = σ_m(ω) /m, and the average for each harmonic degree is . Our final model for the uncertainties on the splittings is a quadratic fit to the uncertainties of the splittings for each ℓ relative to ν_max. Figure 1 shows the uncertainties for the ℓ = 1 modes at an inclination of 90° with ten radial orders evenly distributed about the frequency of maximum power. The uncertainties on the ℓ = 2 modes are half that of the ℓ = 1 modes. We add a Gaussian-distributed random realisation of the noise that is centred at zero and with standard deviation given by the uncertainties to the splittings of each stellar model.

Fig. 1 Uncertainty on the mode splittings ℓ = 1, | m | = 1 for each stellar model, as indicated by the temperature. The vertical ticks indicate the frequency of maximum p-mode power, νmax, in the envelope. Hotter stars have modes with higher oscillation frequencies and larger uncertainties on the splittings.

3. Inversion method: functional fitting

We call forcing a specific function to be fit to the data “functional fitting”. We choose to fit a step function as the simplest estimate of the rotation above and below the base of the convection zone, where the size of the step is a proxy for the globally averaged RDR gradient, from the core to the surface. Traditionally one would fit for the rotation rate below and above the convection zone. We reformulate this slightly so that we fit for the rotation rate below the convection zone and for the size of the step. For a single star the function is (4)where we fit the rotation rate below the base of the convection zone, Ω^∗, and the step size for the coefficients, ΔΩ. In this way we can directly get a measure of the overall gradient of the RDR and its uncertainties. The rotational splitting is linearly related to the coefficients , where a₁ = Ω^∗ and a₂ = ΔΩ, and .

For an ensemble fit to N stars, we fit for independent rotation rates below the convection zone for each star, , but for a common step size for all stars, ΔΩ. The splitting is then the sum over all coefficients, , where .

4. Independent asteroseismic fits

Paper I shows that inversions for stars using different stellar models within the spectroscopic uncertainties have similar uncertainties. To determine over which range of stars this property extends and thus over what range of stellar types an ensemble fit would be beneficial, we compute independent asteroseismic fits for each star. First, we compute the fit for each stellar model independently (N = 1) using only ℓ = 1 modes. Table 2 (Cols. 4 and 5) shows the inverted rotation coefficients and uncertainties for each of the models. We impose a different random noise realisation on the splittings for each star. The uncertainties in the inverted rotation coefficients are large and, although the synthetic rotation rate below the BCZ and the step size are identical for each model, the inversion returns a broad range of Ω^∗ and ΔΩ values, where the step can be negative or positive. If we double the number of splittings to use ℓ = 1, | m | = 1 and ℓ = 2, | m | = 2 modes (since we are using i = 90°), the uncertainties are reduced (Table 3, Cols. 4 and 5) but still not enough to successfully constrain the size or sign of the RDR gradient. As shown in Paper I for regularised least squares inversions, we see that the uncertainties for each of the stellar models across the main sequence are of the same order. This suggests that the uncertainty in the inverted coefficients will reduce as the square root of the number of stars, and all stars across the main sequence can be included in an ensemble inversion.

5. Ensemble asteroseismic fits

We now explore the ensemble fit for all of the 15 stars (N) at once. The ensemble fits of the ℓ = 1 modes for all of our modelled stars are shown in Table 2, Cols. 6 and 7. One thing to note is that the weighted mean of the step size for the independent fits is −216 ± 294 nHz, showing that the ensemble inversion provides slightly better uncertainties than simply averaging the independent inversions.

The uncertainties are reduced from the independent cases and the inverted values for vary less, but the uncertainty is not small enough to constrain the sign of ΔΩ. Using only the ℓ = 1 mode splittings, this method cannot differentiate between a solid body rotation and a rotation gradient for this synthetic profile.

In Table 3 Cols. 6 and 7, we show the results from fits including the ℓ = 2, | m | = 2 mode splittings. The uncertainties in the ensemble fit show that the sign of the step can be correctly constrained since . For this noise realisation, the inverted values for and ΔΩ_∗ are quite close to the input values (Ω_c = 1035 nHz and ΔΩ_c = −350 nHz).

5.1. Variations in the rotation profile

This is a linear inversion and we can only invert for the average step size. To demonstrate the stability of the method, we compute the ensemble fit again, but with the synthetic rotation rate in the core perturbed by a Gaussian distributed random realisation of the noise centred at zero and standard deviation of 50% of 1035 nHz (≈250 nHz) (Tables 2 and 3, Cols. 8 and 9). The variation in values of Ω_c is larger than when the synthetic Ω_c are all the same, reflecting the variation of the synthetic rotation profile in the stars. Varying the rotation rate below the convection zone was a choice to simply demonstrate the behaviour of the fitting method in response to an ensemble of different rotation rates. We could have equally chosen to vary the step size, or both parameters. However, the resulting uncertainties in the parameters would have remained the same.

5.2. Lower radial order mode splittings

We explore the consequences for the inversions by adding the next five lower radial orders to all of the stars, so that we use a total of 30 mode splittings. In Table 4, we show that the uncertainties are reduced, which demonstrates that splittings measured from lower frequency modes better constrain the fit (Fig. 1).

6. Measuring radial differential rotation along the main sequence

Our ensemble fit successfully inferred the sign of the rotation rate for a broad range of stellar models along the main sequence. We now select three ensembles of stars based on spectral type to determine the sensitivity of our method to differences in rotation across the main sequence. This was simply a choice to demonstrate the potential of the ensemble fitting method, however it is not known how or if differential rotation varies systematically across the main sequence. This is something we would like to determine in the future using our ensemble fitting method. Analysing ensembles of stars based on other physical characteristics, such as the bulk rotation rate, activity level, or age, are other options.

We estimate the uncertainty in the ensemble fits using 1000 stars. We estimate the proportion of stellar types using the stellar counts from the Geneva-Copenhagen Survey (Nordström et al. 2004) of bright, main-sequence cool dwarfs. Using the temperatures, we estimate that approximately 60% will be F-type stars (T_eff ≥ 6000 K), 35% will be G-type (5000 ≤ T_eff< 6000 K) stars, and 5% K-type (T_eff< 5000 K) stars. We note that PLATO will make two long duration observations of >1000 bright stars (as in the Geneva-Copenhagen Survey) of stellar type F5 to K7 and 20 000 less bright stars, that will be suitable for measuring rotation (Rauer et al. 2013).

As shown by Paper I and in Sect. 4 of this paper, the uncertainties for all the independent inversions for each star will reduce with the square root of the number of stars used in the ensemble fit. Table 5 shows the mean of the uncertainties for all the models of each stellar type in our sample. If we extrapolate these uncertainties to the number of stars of each type expected to be observed by PLATO, the uncertainties can be significantly reduced (right-hand side of Table 5 ). Therefore, given the uncertainties on the step, σ(ΔΩ_∗), in Table 5, it would be possible to distinguish between the average RDR gradient for different stellar types across the main sequence if the differences are larger than ≈200 nHz.

Thus far, we have only used uncertainties for stars with a rotation axis perpendicular to the line of sight (i = 90°). Assuming stars have a rotation axis randomly inclined to the line of sight, the probability distribution of a star having inclination i is a sini distribution. To estimate the effect of stellar inclination on the uncertainties on the fitting coefficients, we computed the functional fit for each star with inclinations from 10° to 90° and computed the sini weighted average. The uncertainty on the ensemble of each type of star in our sample is shown in Table 6 (left). On the right, we have scaled for the number of stars in a sample of 1000 (Table 6, right). The uncertainties on the fits are larger at lower inclinations and so there is a slightly larger uncertainty than at i = 90°.The smallest uncertainties that are estimated for on the step, σ(ΔΩ_∗), in Table 6 for ℓ = 1,2 modes show that it would be possible to distinguish between the average RDR gradient for different stellar types across the main sequence if the differences are larger than ≈100 nHz. The uncertainties on the RDR themselves allow us to place constraints down to the order of solar RDR (≈30 nHz). The trend in the rotation rate could be constrained if the differences are larger than ≈50 nHz.

7. Conclusion

Supported by the fact that the inversions are insensitive to the uncertainties in the stellar models (Paper I), we found that ensemble asteroseismic fits are a feasible method to extract basic rotation parameters. We found that the uncertainties in the independent RDR fits of F, G, and K-type stars are of the same order (Sect. 4). We demonstrated that for 15 stars across the main sequence it is possible to distinguish between solid body rotation and a RDR gradient of ≈200 nHz using splittings of ℓ = 1,2 modes.

However, ensemble fittings and inversions have limitations. The assumption behind ensemble fitting is that the RDR of the stars in the ensemble behaves systematically. This technique is a linear inversion to retrieve the average RDR step function size. Any outliers will have an effect on the result depending on the number of modes used in the inversion, the level of noise on the splittings, and the magnitude of the outlier. For example, if by some coincidence half of the ensemble stars had a RDR with exactly the opposite sign and magnitude of the other half, then the retrieved RDR from the ensemble fitting would indicate solid-body rotation. Although, the variance of the inferred values for the rotation below the convection zone does reflect the variance of the input values. Since the RDR of stars along the main sequence is not known, and this is what we would like to constrain, there is no prior information to help select a useful ensemble. Therefore, doing ensemble fitting for ensembles of stars based on different physical characteristics, such as bulk rotation rate, activity level, or age, will help us to understand for which physical processes RDR is important.

Further constraints to the inversions for RDR could be implemented by using surface rotation constraints from starspot rotation or vsini measurements. Measuring rotational splittings of lower radial order or higher harmonic degree rotational splittings with SONG (Pallé et al. 2013, a new ground-based network to measure the surface Doppler velocities of stars) may also help to reduce the uncertainties, especially at lower frequencies.

With the advent of PLATO, we will be able to observe thousands of main-sequence Sun-like stars, and begin to trace differences in rotation for subgroups of stellar populations: e.g. stellar type, bulk rotation rate, or activity-related measurements. We showed that with a total of 1000 stars, the difference in magnitude of the gradient can be detected for F, G, and K-type stars if the difference in the average RDR gradient is larger than ≈100 nHz. We also showed that this would place constraints on the average RDR of the ensemble of stars down to the order of solar RDR, and we could detect trends in the rotation rate across the main sequence above ≈50 nHz. Ensemble inversions can also help to constrain the importance of RDR for stellar activity such as absolute activity level, stellar cycle period, and spot coverage.

Acknowledgments

The authors acknowledge research funding by Deutsche Forschungsgemeinschaft (DFG) under grant SFB 963/1 “Astrophysical flow instabilities and turbulence” (Project A18).

References

All Tables

All Figures

	Fig. 1 Uncertainty on the mode splittings ℓ = 1, | m | = 1 for each stellar model, as indicated by the temperature. The vertical ticks indicate the frequency of maximum p-mode power, νmax, in the envelope. Hotter stars have modes with higher oscillation frequencies and larger uncertainties on the splittings.
In the text