Asteroseismology of $\beta$ Cephei stars: The stellar inferences tested in hare and hound exercises

As defined in Eq.(1), the $\chi^{2}$ function would follow a $\chi^{2}_{\mu}$ statistic with $\mu$ degrees of freedom, provided it respects two conditions: first, the frequencies $\nu_{\mathrm{obs,i}}$ have to be independent and secondly, the theoretical model has to be expressed by a linear combination²²2or at least does not depart too strongly from linearity function of the parameters, $a_{j}$, which are adjusted to obtain the best fit to the observations. These conditions are not met in the case of $\beta$ Cephei pulsators. Indeed, the frequencies $\nu_{\mathrm{th,i}}\equiv\nu_{\mathrm{th,i}}(a_{1},...,a_{j},...a_{n})$, rely on global quantities symbolised by the $n$ parameters $a_{j}$, with $n=5$. These five parameters are $M$, $X$, $Z$, $\alpha_{\mathrm{ov}}$, and the age, with $X$ the initial hydrogen mass fraction, and $Z$ the initial metallicity. For instance, in the simple case of radial oscillation modes, the relation of their frequencies $\nu_{\mathrm{th}}\propto 1/\tau_{\mathrm{dyn}}$, with the dynamical timescale $\tau_{\mathrm{dyn}}=\sqrt{R^{3}/GM}$ ($R$ being the radius, and $G$ the gravitational constant), already reveals a non-linear dependence on $M$.

More generally, due to the mixed character of the low radial order modes (using $k$ for the radial order hereafter) in $\beta$ Cephei stars, the frequency spectra present features known as avoided crossings. The mixed modes present both an acoustic (p-) and gravity (g-) mode behaviour and appear when the propagation cavities of p- and g-modes are only separated by a small evanescent region. As the $\beta$ Cephei stars develop a chemical composition gradient as their convective core recedes, a sharp and localised peak in the Brunt-Väisälä frequency appears, favouring the presence of mixed modes. The avoided-crossing phenomenon is shown in Fig. 1 in light blue in the main-sequence evolution of the dimensionless ³³3$\tau_{\textrm{dyn}}$ is used for the nondimensionalisation frequencies ($\omega$) of low-radial order modes of a typical $\beta$ Cephei model. The sharp variations of the frequencies induce a non-linear relation between the frequencies and the parameters of the model (e.g. Scuflaire, 1974; Shibahashi, 1979).

Thus we use a re-sampling method to explore the accuracy and precision on the stellar parameters derived from the selected best-fit model. We perform Monte-Carlo simulations by randomly drawing frequencies that we call pseudo-observed frequencies. Assuming $N_{obs}$ frequencies have been observed for a given star, we build $M_{\mathrm{sim}}$ new sets, each composed of $N_{obs}$ pseudo-observed frequencies. These latter are drawn one-by-one from Gaussian distributions of standard deviation, $\sigma_{\mathrm{\nu,i}}$ centred on the $\nu_{\mathrm{obs,i}}$. The values of $\sigma_{\mathrm{\nu,i}}$ are tailored not only to account for the observational errors, but also the uncertainties in the physics of the models. A natural choice is hence to set $\sigma_{\mathrm{\nu,i}}=\max(\sigma_{\nu_{\mathrm{obs,i}}},\sigma_{\nu_{\mathrm{theor,i}}})$, where $\sigma_{\nu_{\mathrm{obs,i}}}$ and $\sigma_{\nu_{\mathrm{theor,i}}}$ are the observational and theoretical errors, respectively (see details in Sect. 2.2).

A best-fit stellar model (based on the minimum in $\chi^{2}$) is then derived for each of the $M_{\mathrm{sim}}$ simulated sets of pseudo-observed frequencies. The stellar parameters of the best-fit model of each simulation are collected. The distributions of stellar parameters obtained in this way are then considered as the distributions of the solution in the space of parameters, accounting for the observational errors and theoretical uncertainties of the models (see Sect. 2.2).

Indeed, we read these distributions built from the simulations as probability density distributions. The medians of the marginal distribution of a given stellar parameter are picked up as the inferred results from the Monte-Carlo (MC hereafter) simulations. We also derive confidence intervals on the estimated stellar parameters. The lower and upper limits of the confidence intervals at the 1-$\sigma$ level, conventionally set to 68.3%, correspond to the $(1-\sigma)/2\times M_{\mathrm{sim}}$ and $(1+\sigma)/2\times M_{\mathrm{sim}}$ values of the ordered marginal distribution.

We can then compare these confidence intervals to the single set of stellar parameters obtained from the minimum in $\chi^{2}_{\mathrm{min,true}}$, which in this case is based on the actual observed frequencies.

By carrying out tests with assumed different numbers of observed frequencies (keeping it representative of typical $\beta$ Cephei observations) and numbers of simulations, we found that increasing the number of simulations above $M_{\mathrm{sim}}$=1000 did not change the estimated errors nor did they induce changes in the distributions of the inferred stellar parameters. We consequently adopted that value in our different modelling exercises.

To determine the value of $\sigma_{\mathrm{\nu,i}}$ to run the simulations, we review the expected observational and theoretical errors affecting $\beta$ Cephei frequencies. The observational errors are typically of 10^-4-10^-6 cycle per day (c/d) for well-studied $\beta$ Cephei stars, see for instance Handler et al. (2004) and Aerts et al. (2004) in the case of ground-based spectroscopy and photometry, respectively; or Degroote et al. (2012) for the case of space-based photometry.

Different theoretical uncertainties contribute to the theoretical error. We first tested the uncertainty in the micro-physics of the model by comparing the adiabatic frequencies between stellar models with identical stellar parameters, but differing in the chemical mixture and opacity. Their theoretical pulsation frequencies varied by $\sim 10^{-3}$–10^-2 c/d. On the scale of macro-physics, stellar rotation, not included in the models of our grid, is likely one of the processes potentially impacting at most the stellar structure, but also impacting directly the physics of oscillations. Including the effect of rotation on oscillations and on the stellar structure (i.e. distortion), Briquet et al. (2012b) found a theoretical difference up to $\sim$ 5$\times$ 10^-2 c/d.

The computed frequencies of oscillation also depend whether the energy equation is treated under the adiabatic assumption. However, given the non-adiabaticity of superficial stellar layers, a non-adiabatic approach can reveal shifts in the frequencies. We hence computed adiabatic and non-adiabatic frequencies for a sample of models representative of the theoretical grid. We used for this purpose the LOSC ( Liège OScillation Code) adiabatic and MAD non-adiabatic codes, in their default setting mode, as detailed respectively in Scuflaire et al. (2008a) and Dupret (2001). We found differences between adiabatic and non-adiabatic frequencies from $\sim$ 10^-4 to 10^-2 c/d, depending on the $\ell$ degree and evolutionary stage. We finally took the maximum of these differences for $\sigma_{\nu_{\mathrm{theor,i}}}$, that is 10^-2 c/d, as the value representative of the theoretical uncertainties. The observational errors on frequencies are typically two to three orders smaller than the theoretical ones, so that the limits on asteroseismic precision shall be dominated by current model uncertainties.

To test the potential of classical constraints, we need to estimate what are the typical values of their observational errors in the case of $\beta$ Cephei stars. The fundamental parameters are usually derived from spectroscopic observations. In particular, the effective temperature of early B-type stars can be determined from the ionisation balance of spectroscopic lines of He, C, O, Ne, Si or Fe. When a maximum of these indicators is detected, the error on $T_{\mathrm{e}}$ can be as low as 200-300 K (e.g. Nieva & Przybilla, 2012). However all these lines are not always available: instead strong lines of Ne and Si are then used, leading to an error of $\sim$1000 K (T. Morel, private communication). We adopt as the error on $T_{\mathrm{e}}$ an intermediate value of 700 K for the hare and hound exercises. The surface gravity $\log g$, is estimated with the help of a fit to the wings of H Balmer lines. In that case, a minimum 0.15 dex reasonably accounts for the possible sources of error (T. Morel, private communication). We note that depending of the cases, for examples in reason of fast rotation or radial velocity variability in a binary system, the error can be larger, around 0.30 dex.

The spectra of early-type stars present many lines of N, O and Fe, and to a lower extent of C, enabling to determine their photospheric abundances. Those of Ne, Mg, Si, and S can also be detected under favourable conditions, then giving an insight on the most abundant elements. The present-day ratio $Z/X$ can then be calculated by assuming the abundances of minor elements follow the solar mixture. Different studies have revealed that abundances of B stars in the solar neighbourhood differ and are metal poorer than the solar abundances (Przybilla et al., 2008; Morel, 2009; Niemczura et al., 2009; Nieva & Simón-Díaz, 2011). The difference was particularly marked with past solar abundance determinations such as the Grevesse & Noels (1993), hereafter GN93, one. Revision of solar abundances led to an important decrease of the solar metallicity in 2005 (Asplund et al., 2005, hereafter AGS05). A subsequent revision (Asplund et al., 2009) has moderated this downward revision of the solar metallicity, but the AGS05 appears today the more representative of the chemical mixture of neighbour B stars. Once the solar distribution assumed, obtaining the metallicity from $Z/X$ still requires to determine $X$, for instance from a solar calibrated model. Morel et al. (2006) derived the metallicity in this way for a sample of $\beta$ Cephei stars, finding an average error on $Z$ of $\sim$ 0.002.

If the number of individual abundances is insufficient, metallicity is derived from the relation [Fe/H]=[$Z/X$] that assumes the stellar abundances are distributed following the solar mixture. In the literature we find errors on the Fe abundance of early-type stars from $\sim$0.1 to 0.2 dex (e.g. Morel et al., 2006; Nieva & Przybilla, 2012). This translates to typical errors on the metallicity $Z$ of $\sim$0.003 and 0.005. As we see in the next sections, the uncertainty on the metallicity derived with help of the MC simulations are typically of $\sim$0.002–0.004, except for of a few cases.

The adiabatic frequencies of oscillation are determined for all the models from the Zero Age Main Sequence (ZAMS) to the Terminal Age Main Sequence (TAMS), with the help of the LOSC tool, adopting the standard boundary conditions of the code (see exact details in Scuflaire et al., 2008a). For each model on the main sequence, the radial ($\ell$=0) and non-radial ($\ell$=1,2,3) low-order p- and g-modes are computed. For each model and each $\ell$, the modes are computed in the same fixed interval of adimensional frequencies, chosen as representative of $\beta$ Cephei stars.

The density of models in a restricted region of the grid is shown in Fig. 2 (the whole grid ranges in $T_{\textrm{e}}$ from $\sim$ 14,400 to $\sim$ 36,000 K and surface gravities from 260 to 22,400 cm/s², or $\log g\sim$ 2.40 – 4.35). There are typically several dozens of models per 50 K and per 20 cm/s² boxes. The denser parts correspond to models in the end of the main sequence, and which overlap with less massive models presenting higher values of $\alpha_{\textrm{ov}}$. Overshooting not only extends the MS lifetime but increases as well the luminosity of evolutionary tracks at the end of MS so that the models with high overshooting cross the tracks of more massive ones.

We selected for this series of tests a model computed with the same physics than in the grid, with the goal of isolating the only influence of the seismic indicators. Contrary to the grid, the frequencies of this t1 model are computed with the MAD non-adiabatic code, as explained in the previous section.

We first analysed the impact of misidentifying a $\ell$ degree. In practice, although the methods⁵⁵5such as the moment method based on spectroscopic line variations (Balona, 1986; Aerts, 1996) or the analysis of photometric bandpass ratio Cugier et al. (1994). developed for identifying the modes in $\beta$ Cephei stars are generally able to constrain the angular degree of detected pulsations, in some cases the identification can be ambiguous. That was for example the case of $\nu$ Eri (see De Ridder et al., 2004, in particular their Figs. 5 and 6), whatever the method used. The same difficulty occurred for the photometric bandpass ratios analysis of the 12 Lac star (Handler et al., 2006). In a second time, we have focused on cases where no identification of the mode is possible, which is typically the case of stars observed in a single photometric bandpass.

We changed the physics of the model used as a target so that it differed from that of the grid. We examined the consequences on the quality of the global asteroseismic parameters, as well as on the constraints on the internal structure. We have seen in Sect.5.1.1 that there is a significant degeneracy for the determination of overshooting and metallicity. This degeneracy will be particularly accentuated if the detailed chemical composition of the star differs from that used in theoretical models. For instance, we recalled in Sect. 2.2 that between past determination of solar abundances like GN93 and the revision of AGS05, the metallicity decreased by $\sim$30%. Hence we tailored the t2-gn93 exercise to test the impact of the adopted composition. We used for the target t2 model the solar mixture of GN93, while it is the AGS05 mixture that is adopted in the grid. The other properties of the t2 models are given in Table 2. The set of input frequencies for the exercise are reported in Table 3 and is composed of six frequencies, with pairs of modes of degree $\ell_{0}$,$\ell_{1}$, and $\ell_{2}$.

The output parameters slightly overestimate the parameters of the t2 models with $M$=11.4 $\textrm{M}_{\odot}$ and $R$=6.07 $\textrm{R}_{\odot}$. The best-fit model reproduces $\tau_{\textrm{dyn}}$ within 1% the one of the t2 model, confirming the influence of radial modes on the modelling. We see in the top panel of Fig. 11 that the lowest values of the merit function are indeed located in ridges, which actually correspond to places of equal $\tau_{\textrm{dyn}}$. The parameters $\alpha_{\textrm{ov}}$ and $Z$ are overestimated and underestimated, respectively (see column 3 of Table 7). This comes as expected from the degeneracy in $\alpha_{\textrm{ov}}$ and $Z$, accentuated by the difference in the chemical mixture. This effect was already observed in the modelling of HD129929, then from an ad hoc change in the mixture made by the authors (Thoul et al., 2004). Actually, in reason of the predominance of nuclear energy production by CNO cycle in B stars, for a given X and Z, a GN93 model will be more luminous than an AGS05 one. The reason is that C, N, and O are more abundant constituents of the metal mixture in the first case. Since we try to model a GN93 model with AGS05 models, the Z and $\alpha_{\textrm{ov}}$ likely adapt and do not correspond to those of the GN93 model.

Including the classical constraints. With predicted $Z$ = 0.014 and $\alpha_{\textrm{ov}}$ = 0.25 (using 1-$\sigma$ error box, see column 5 of Table 7), the discrepancy between the solution and the t2 model is reduced, but $X$ is now clearly overestimated, illustrating again the degeneracy on the overshooting and chemical composition determinations. The results of the re-sampling method are not very different (column 6 of Table 7), although the errors on the parameters are larger and so enclose the actual values of the t2 model. In Fig. 11, the global minimum and the lowest values of $\chi^{2}$ surrounding it, lie outside the 1-$\sigma$ error box on the classical parameters (top panel). In the bottom panel, the local minimum found in the 1-$\sigma$ error box is not well defined and its surrounding $\chi^{2}$ values are very similar: the solution is not well constrained, resulting in larger uncertainties on the derived parameters. The 3-$\sigma$ constraints do not bring any information as the global minimum lies within it. This global minimum is clearly brighter and hotter than the target t2 model.

Overall, the inaccuracy induced by the difference in the chemical mixture remains limited to a few % on global parameters such as $M$ and $R$. However, it is difficult to estimate the accuracy of parameters such as metallicity and overshooting given the degeneracy and compensatory role they play. Verifying how the limits of the convective core and $\nabla\mu$ are reproduced would better hint at the robustness of the asteroseismic modelling on this question. We present in Fig. 12 the internal profile of $X$ of the t2 and its seismic best-fit models, showing that the size in mass of the fully mixed regions and $\nabla\mu$ are not retrieved by the modelling (accuracy errors of 16 and 12%, respectively). Moreover, one mode ($k_{1}$, $\ell_{1}$) of the input dataset is mixed with a dominant g-character, in principle an optimal information on these stellar layers. Yet, we see from Table 5 that the accuracy on these limits is much better in terms of $r/R$, $\sim$5%. It highlights that the asteroseismic information is defined in terms of variables related to oscillations, so essentially sensitive to $r$. While the acoustic structure of the star can be well reproduced, this is not necessarily the case of the mass distribution in the layers. Here it is clearly the difference in the chemical mixture that hinders the recovery of mass of the convective core and its overshooting region.

We did two other tests where the target models were computed with the GN93 mixture and with OPAL opacities (versus OP in the grid), but in one case, no overshooting was included. In both cases, the additional change in the physics did not alter precedent conclusion: we were still able to infer $M$ and $R$ with a good accuracy. But in the case with overshooting, we faced the same degeneracy on the chemical composition and overshooting determination. It resulted in an important discrepancy of $\sim$20% on the $m/M$ limit of the fully-mixed central region. Nevertheless, in the case without overshooting, the solution correctly retrieved the absence of overshooting (but failed at the chemical composition). The limits on the convective core were still of $\sim$5% in terms of $r/R$, but remained of the same order when expressed in $m/M$.

We have considered insofar models where the extra mixing was treated with a classical instantaneous prescription (see e.g. Maeder, 1975). Other processes can be responsible for extra mixing near the convective core, as for example turbulent mixing induced by rotation. Although this latter process has almost the same impact on the stellar evolutionary tracks as the overshooting (see e.g. Talon et al., 1997), the process acts like a diffusive process leading to smoother chemical composition gradients at the boundary of the convective core. The effect is in principle noticeable in comparison to the sharp profile generated in that same region by a very efficient mixing.

We explored this question in the t6-asp-diff exercise, in which turbulent mixing is implemented as a diffusive process parametrised by a turbulent diffusion coefficient, $D_{\textrm{t}}$. We computed the t6 target model (see Table 2) with a value of $D_{\textrm{t}}$ = 70,000 cm.s^-2, calibrated to correspond to a model computed with the Geneva code (Eggenberger et al., 2008) with an initial equatorial rotational velocity on the ZAMS ($V_{\textrm{i}}$) of 50 km/s. Doing so, the t6 model is equivalent (in terms of the extent of the central fully-mixed region) to a model without diffusive mixing but $\alpha_{\textrm{ov}}$ = 0.05.

The input frequency set is composed of two radial modes, three $\ell_{1}$, and two $\ell_{2}$ modes. The modelling with these seismic constraints recovers with good accuracy $M$, $R$, and $\alpha_{\textrm{ov}}$ (see column 3 of Table 7), but $X$ and $Z$ are overestimated.

The re-sampling method does not predict different results and it fails to determine consistent errors on the parameters (see column 4 of Table 7). The reason is that the global minimum is deeply marked in the $\chi^{2}$ map in Fig. 13, so that other solutions hardly emerge from the MC simulations. The global minimum lies within the 1-$\sigma$ error box on $T_{\textrm{e}}$–$\log g$, and would do the same considering typical 1-$\sigma$ observational error bars on $Z/X$. The classical constraints appear as no help in this case.

Of prime importance, we finally look at the limits of the convective core and $\nabla\mu$ in Fig. 14. The limit of the convective core is determined with a precision of  5% (see also Table 5) in $m/M$, as it could be expected, since the t6 target model and the grid share the same chemical mixture (AGS05). However, the modelling fails at determining the location of $\nabla\mu$, and so to be sensitive to the nature of the extra mixing. This is somehow expected as the grid did not include any diffusive mixing, confirming that getting constraints on the internal processes depends on the stellar models used. It suggests that a diffusive mixing should be as well considered in the grid for assessing its presence and efficiency in $\beta$ Cephei stars. This would add a degree of freedom, under the form of a diffusive coefficient taken as an additional parameter to adjust. Lovekin & Goupil (2010) tested the addition of a parameter in the fit, although applied to retrieve the rotation velocity of the star. Their results were encouraging as the convergence towards a reliable modelling of $\beta$ Cephei stars was not hindered by the additional parameter to fit.

We also did another test (not detailed here) considering a target model with diffusive mixing and the GN93 mixture, still with the grid only including instantaneous mixing. As expected, the inaccuracy increased, in particular on the determination of the convective core. Its inferred mass was then underestimated by 15%. The solution also exceeded by 15% the limit in $m/M$ on $\nabla\mu$, confirming the insensitivity to the nature of extra mixing with an inappropriate model grid. Despite the internal structure was less-well characterised, the global stellar parameters were recovered with the same accuracy and precision than it the tests where the target stars included no diffusive mixing and the GN93 mixture in Sect. 5.3.